{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":117682,"databundleVersionId":14443416,"sourceType":"competition"}],"dockerImageVersionId":31192,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"### Input data standardization\n#### Brett Olsen\n#### December 11, 2025\n\nFrom our previous look into the Vesuvius data, it's clear that we'll want to normalize the data to a common standard before processing it.  In particular, we'll want to:\n\n1.  Denoise the data\n2.  Normalize it so that background/foreground signal is matched to a common value\n\nLet's get started.","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"# These are used for efficiently loading the 3D tiff files\n!pip install imagecodecs\n!pip install tifffile\n# This will be used for contrast normalization\n!pip install kmeans1d","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-11T19:06:01.760456Z","iopub.execute_input":"2025-12-11T19:06:01.761123Z","iopub.status.idle":"2025-12-11T19:06:13.959260Z","shell.execute_reply.started":"2025-12-11T19:06:01.761094Z","shell.execute_reply":"2025-12-11T19:06:13.957925Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os\n\nimport scipy\nfrom scipy.special import erf\n\nfrom skimage import img_as_float\nfrom skimage.restoration import denoise_nl_means, estimate_sigma\nimport skimage.exposure\n\nfrom collections import Counter\nfrom tifffile import tifffile\nimport imagecodecs\n\nimport kmeans1d\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nsns.set_context(\"poster\")\n\nINPUT_DIR = \"/kaggle/input/vesuvius-challenge-surface-detection/\"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-11T19:06:13.961752Z","iopub.execute_input":"2025-12-11T19:06:13.962507Z","iopub.status.idle":"2025-12-11T19:06:13.971213Z","shell.execute_reply.started":"2025-12-11T19:06:13.962471Z","shell.execute_reply":"2025-12-11T19:06:13.970088Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def remove_zeropadding(data, labels):\n    \"\"\"Our training data is coming in with some zero padding.\n    We're going to filter down to just the regions that aren't\n    entirely zeros, saving the offsets.\n    \"\"\"\n    offsets = [0, 0, 0]\n    \n    idxs = np.where(np.sum(data, axis=(1, 2)) > 0)[0]\n    zslice = slice(idxs[0], idxs[-1] + 1, None)\n    offsets[0] = idxs[0]\n\n    idxs = np.where(np.sum(data, axis=(0, 2)) > 0)[0]\n    yslice = slice(idxs[0], idxs[-1] + 1, None)\n    offsets[1] = idxs[0]\n\n    idxs = np.where(np.sum(data, axis=(0, 1)) > 0)[0]\n    xslice = slice(idxs[0], idxs[-1] + 1, None)\n    offsets[2] = idxs[0]\n\n    if labels is not None:\n        labels = labels[zslice, yslice, xslice]\n    \n    return data[zslice, yslice, xslice], labels, tuple(offsets)\n\ndef load_data(section_id, data_type=\"train\"):\n    if data_type == \"train\":\n        path = os.path.join(INPUT_DIR, \"train_images\", f\"{section_id}.tif\")\n        label_path = os.path.join(INPUT_DIR, \"train_labels\", f\"{section_id}.tif\")\n    elif data_type == \"test\":\n        path = os.path.join(INPUT_DIR, \"test_images\", f\"{section_id}.tif\")\n        label_path = None\n    else:\n        raise ValueError(f\"Unknown data type {data_type}\")\n    if not os.path.exists(path):\n        raise ValueError(f\"File does not exist: {path}\")\n    print(f\"Loading file {path}\")\n    data = tifffile.imread(path)\n    print(f\"Data shape {data.shape}\")\n    if label_path is not None:\n        labels = tifffile.imread(label_path)\n    else:\n        labels = None\n    data, labels, offsets = remove_zeropadding(data, labels)\n    return data, labels, offsets","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-11T19:06:13.972439Z","iopub.execute_input":"2025-12-11T19:06:13.972813Z","iopub.status.idle":"2025-12-11T19:06:14.001486Z","shell.execute_reply.started":"2025-12-11T19:06:13.972788Z","shell.execute_reply":"2025-12-11T19:06:14.000325Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"df = pd.read_csv(os.path.join(INPUT_DIR, \"train.csv\"))\nsection_ids = df[\"id\"].values\n\n# Two samples I know have different distributions\nsection_1 = 1407735\nsection_2 = 567464717\n\ndata1, _, offsets = load_data(section_1)\nprint(offsets)\ndata2, _, offsets = load_data(section_2)\nprint(offsets)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-11T19:06:14.003510Z","iopub.execute_input":"2025-12-11T19:06:14.003896Z","iopub.status.idle":"2025-12-11T19:06:15.510907Z","shell.execute_reply.started":"2025-12-11T19:06:14.003865Z","shell.execute_reply":"2025-12-11T19:06:15.509725Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Let's look at the distribution of signal across all 256 (uint8 data) values in both datasets.","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(12, 4))\nplt.grid()\nbins = np.arange(256)\n\nplt.hist(data1.flatten(), bins=bins, histtype=\"step\", alpha=0.5);\nplt.hist(data2.flatten(), bins=bins, histtype=\"step\", alpha=0.5);\nplt.xlabel(\"Signal\")\nplt.ylabel(\"Voxel Count\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-11T19:06:18.452410Z","iopub.execute_input":"2025-12-11T19:06:18.452771Z","iopub.status.idle":"2025-12-11T19:06:21.449694Z","shell.execute_reply.started":"2025-12-11T19:06:18.452747Z","shell.execute_reply":"2025-12-11T19:06:21.448651Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Great, we can very clearly see the differences between these two:\n\n1. There's a large pileup at a signal of 0.  That suggests there's already _been_ a thresholding at some point upstream, at least for one of the datasets.\n2. There are loosely two separate peaks, a larger one at lower signal and a smaller one at higher signal.  Likely corresponding to foreground and background.\n3. We can see a small spike at maximum signal, likely from the high-contrast inclusions we've seen elsewhere.\n\n---\n\nLet's start by denoising.  We'll just use scipy tools for this; they'll be straightforward.  I like non-local means for this.  It scales well to 3D data, though you have to be careful with `patch_distance`.","metadata":{}},{"cell_type":"code","source":"def denoise(data):\n    # First we cast the image to a float.  This renorms so the data range covers[0, 1] rather than [0, 255]\n    data = img_as_float(data)\n    sigma = estimate_sigma(data)\n\n    # This gets _much_ slower with longer patch distances; we can get most of the benefits with a very short distance.\n    # May want to look into a faster implementation.\n    denoised = denoise_nl_means(data, patch_size=7, patch_distance=1, sigma=sigma)\n\n    return denoised","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-11T19:06:25.427068Z","iopub.execute_input":"2025-12-11T19:06:25.427382Z","iopub.status.idle":"2025-12-11T19:06:25.433429Z","shell.execute_reply.started":"2025-12-11T19:06:25.427359Z","shell.execute_reply":"2025-12-11T19:06:25.431920Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\ndenoised1 = denoise(data1)\ndenoised2 = denoise(data2)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-11T19:06:27.198462Z","iopub.execute_input":"2025-12-11T19:06:27.199660Z","iopub.status.idle":"2025-12-11T19:07:09.798553Z","shell.execute_reply.started":"2025-12-11T19:06:27.199623Z","shell.execute_reply":"2025-12-11T19:07:09.797502Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Cast the raw data to floats for comparisons\ndata1 = img_as_float(data1)\ndata2 = img_as_float(data2)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-11T19:07:09.800258Z","iopub.execute_input":"2025-12-11T19:07:09.800521Z","iopub.status.idle":"2025-12-11T19:07:10.053884Z","shell.execute_reply.started":"2025-12-11T19:07:09.800499Z","shell.execute_reply":"2025-12-11T19:07:10.053003Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Just pick an arbitrary slice.  In general z gives more useful slices for visualization due to the structure of the scrolls.\nz = 150\n\nfor data, denoised in zip([data1, data2], [denoised1, denoised2]):\n    plt.subplots(1, 3, figsize=(15, 5))\n    cmap = \"gray\"\n    \n    plt.subplot(1, 3, 1)\n    plt.imshow(data[z], cmap=cmap)\n    plt.title(\"Original\")\n    plt.subplot(1, 3, 2)\n    plt.imshow(denoised[z], cmap=cmap)\n    plt.title(\"Denoised\")\n    plt.subplot(1, 3, 3)\n    plt.imshow(denoised[z] - data[z], cmap=cmap)\n    plt.title(\"Diff\")\n    \n    plt.tight_layout()\n\n    plt.figure(figsize=(12, 4))\n    plt.grid()\n    \n    plt.hist(data.flatten(), bins=256, histtype=\"step\", label=\"Raw\");\n    plt.hist(denoised.flatten(), bins=256, histtype=\"step\", label=\"Denoised\");\n    plt.xlabel(\"Signal\")\n    plt.ylabel(\"Voxel Count\")\n    plt.legend(loc=0)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-11T19:07:10.054850Z","iopub.execute_input":"2025-12-11T19:07:10.055142Z","iopub.status.idle":"2025-12-11T19:07:15.709022Z","shell.execute_reply.started":"2025-12-11T19:07:10.055112Z","shell.execute_reply":"2025-12-11T19:07:15.707997Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Great.  The denoising clearly helped clean up the signal.  Next I want to handle contrast adjustment.  There are two clear peaks in both datasets:  one for background signal and one for foreground.  But each dataset has background/foreground set to different values.  So let's figure out how to normalize them to the same range.\n\nWe'll see how well getting the centroids of a 2-means clustering works.","metadata":{"execution":{"iopub.status.busy":"2025-12-08T21:55:55.092680Z","iopub.execute_input":"2025-12-08T21:55:55.092991Z","iopub.status.idle":"2025-12-08T21:55:55.097130Z","shell.execute_reply.started":"2025-12-08T21:55:55.092967Z","shell.execute_reply":"2025-12-08T21:55:55.096333Z"}}},{"cell_type":"code","source":"%%time\ndraw = denoised1.flatten()\nnp.random.shuffle(draw)\n_, centroids = kmeans1d.cluster(draw[:len(draw) // 1000], 2)\nprint(centroids)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-11T19:13:40.567433Z","iopub.execute_input":"2025-12-11T19:13:40.567798Z","iopub.status.idle":"2025-12-11T19:13:42.668779Z","shell.execute_reply.started":"2025-12-11T19:13:40.567772Z","shell.execute_reply":"2025-12-11T19:13:42.667766Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"for data, denoised in zip([data1, data2], [denoised1, denoised2]):\n\n    draw = denoised.flatten()\n    np.random.shuffle(draw)\n    _, centroids = kmeans1d.cluster(draw[:len(draw) // 1000], 2)\n    print(centroids)\n    \n    plt.figure(figsize=(12, 4))\n    plt.grid()\n    \n    plt.hist(data.flatten(), bins=256, histtype=\"step\", label=\"Raw\");\n    plt.hist(denoised.flatten(), bins=256, histtype=\"step\", label=\"Denoised\");\n    for x in centroids:\n        plt.axvline(x, color='k')\n    \n    \n    plt.xlabel(\"Signal\")\n    plt.ylabel(\"Voxel Count\")\n    plt.legend(loc=0)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-11T19:20:01.737496Z","iopub.execute_input":"2025-12-11T19:20:01.738476Z","iopub.status.idle":"2025-12-11T19:20:09.548950Z","shell.execute_reply.started":"2025-12-11T19:20:01.738445Z","shell.execute_reply":"2025-12-11T19:20:09.547839Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Cool.  That's doing a reasonable estimate of the two populations.  Let's go ahead and write a renorming function to adjust contrast for these.","metadata":{}},{"cell_type":"code","source":"BACKGROUND_LEVEL = 0.25\nFOREGROUND_LEVEL = 0.50\n\ndef renorm_image_contrast(data):\n    # Do kmeans fitting on the flattened data with k=2 to get centroids of background/foreground\n    draw = data.flatten()\n    np.random.shuffle(draw)\n    _, (current_background, current_foreground) = kmeans1d.cluster(draw[:len(draw) // 1000], 2)\n\n    scaled = (data - current_background) * (FOREGROUND_LEVEL - BACKGROUND_LEVEL) / (current_foreground - current_background) + BACKGROUND_LEVEL\n    scaled = np.clip(scaled, 0, 1)\n    return scaled\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-11T19:27:22.253033Z","iopub.execute_input":"2025-12-11T19:27:22.253442Z","iopub.status.idle":"2025-12-11T19:27:22.260247Z","shell.execute_reply.started":"2025-12-11T19:27:22.253390Z","shell.execute_reply":"2025-12-11T19:27:22.259066Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\nrenormed1 = renorm_image_contrast(denoised1)\nrenormed2 = renorm_image_contrast(denoised2)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-11T19:27:46.381995Z","iopub.execute_input":"2025-12-11T19:27:46.382298Z","iopub.status.idle":"2025-12-11T19:27:51.584997Z","shell.execute_reply.started":"2025-12-11T19:27:46.382276Z","shell.execute_reply":"2025-12-11T19:27:51.584054Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"z = 150\n\nfor data, denoised, renormed in zip([data1, data2], [denoised1, denoised2], [renormed1, renormed2]):\n\n    plt.subplots(1, 3, figsize=(15, 5))\n    cmap = \"gray\"\n    \n    plt.subplot(1, 3, 1)\n    plt.imshow(data[z], cmap=cmap, vmin=0, vmax=1)\n    plt.title(\"Original\")\n    plt.subplot(1, 3, 2)\n    plt.imshow(denoised[z], cmap=cmap, vmin=0, vmax=1)\n    plt.title(\"Denoised\")\n    plt.subplot(1, 3, 3)\n    plt.imshow(renormed[z], cmap=cmap, vmin=0, vmax=1)\n    plt.title(\"Contrast-Adjusted\")\n    \n    plt.tight_layout()\n    \n    plt.figure(figsize=(12, 4))\n    plt.grid()\n    \n    plt.hist(data.flatten(), bins=256, histtype=\"step\", label=\"Raw\");\n    plt.hist(denoised.flatten(), bins=256, histtype=\"step\", label=\"Denoised\");\n    plt.hist(renormed.flatten(), bins=256, histtype=\"step\", label=\"Contrast-Adjusted\")\n    \n    plt.xlabel(\"Signal\")\n    plt.ylabel(\"Voxel Count\")\n    plt.yticks([])\n    plt.legend(loc=0)\n    plt.xlim((0, 1))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-11T19:31:44.578775Z","iopub.execute_input":"2025-12-11T19:31:44.579127Z","iopub.status.idle":"2025-12-11T19:31:51.171426Z","shell.execute_reply.started":"2025-12-11T19:31:44.579100Z","shell.execute_reply":"2025-12-11T19:31:51.170582Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Great.  You can see that the first sample was lighter than the target contrast and so signal got pulled down, while the second one was darker and so got pulled up.\nThis means that now we can use those constants effectively across all our samples for identifying background and foreground regions without worrying about sample-to-sample variation.","metadata":{}},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}