From 14ccb7eca2aa9e00d582d0751ccc385cfa128a46 Mon Sep 17 00:00:00 2001 From: Bradley Lowekamp Date: Wed, 9 Jun 2021 12:55:30 -0400 Subject: [PATCH 1/4] Add Python Jupyter notebook with prototype of visualization range --- notebooks/AutoMinMaxStats.ipynb | 218 ++++++++++++++++++++++++++++++++ 1 file changed, 218 insertions(+) create mode 100644 notebooks/AutoMinMaxStats.ipynb diff --git a/notebooks/AutoMinMaxStats.ipynb b/notebooks/AutoMinMaxStats.ipynb new file mode 100644 index 0000000..cf30b90 --- /dev/null +++ b/notebooks/AutoMinMaxStats.ipynb @@ -0,0 +1,218 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Histogram of EM Data with Automatic Min/Max Visualization \n", + "\n", + "This notebook is for exploration of automatic computation of Min/Max parameters for adjust with image intestity range used to visualized electron microscopy image data." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "input_filename = \"/Users/blowekamp/Desktop/20130412-Vn-poxtomo_8/20130412-Vn-poxtomo_8_full_rec.mrc\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import SimpleITK as sitk\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from scipy.stats import norm\n", + "%matplotlib notebook" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "img = sitk.ReadImage(input_filename)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "stats = sitk.StatisticsImageFilter()\n", + "stats.Execute(img)\n", + "print(stats)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "a_img = sitk.GetArrayViewFromImage(img).flatten()\n", + "bins = 2**np.iinfo(a_img.dtype).bits\n", + "h, b = np.histogram(a_img, bins=bins, range=(np.iinfo(a_img.dtype).min, np.iinfo(a_img.dtype).max + 1))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def rescale_min_max(img, vmin, vmax):\n", + " ss_filter = sitk.ShiftScaleImageFilter()\n", + " ss_filter.SetShift( -float(vmin) )\n", + " ss_filter.SetScale( 255.0/(vmax-vmin) )\n", + " ss_filter.SetOutputPixelType(sitk.sitkUInt8)\n", + " return ss_filter.Execute(img)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "def plot_range_options(img, range_mult=3, percentile_crop=10, bins=1024, fig_filename=None):\n", + "\n", + "\n", + " fa_img = sitk.GetArrayViewFromImage(img).flatten()\n", + " median = np.median(fa_img)\n", + " mad = np.median(np.abs(fa_img-median))\n", + " mean = np.mean(fa_img)\n", + " sigma = np.sqrt(np.var(fa_img))\n", + "\n", + "\n", + " mad_range = (median-range_mult*mad, median+range_mult*mad)\n", + " percentile_range = (np.percentile(a_img,percentile_crop*.5), np.percentile(a_img,100-percentile_crop*.5))\n", + "\n", + " \n", + " sigma_range = (mean-range_mult*sigma, \n", + " mean+range_mult*sigma)\n", + " \n", + " min_max_range = (np.min(fa_img), np.max(fa_img))\n", + "\n", + " fig = plt.figure(figsize=(8, 8), dpi=240)\n", + " ax = fig.add_subplot(2, 4, (1, 4))\n", + "\n", + " x = np.linspace(min(fa_img), max(fa_img), bins)\n", + " plt.plot(x, norm.pdf(x, mean, sigma), color='k', linestyle='--', alpha=0.5, label=f'Fit Gaussian')\n", + "\n", + " plt.hist( fa_img, bins=bins, density=True, color='b')\n", + "\n", + " plt.axvline(x=mad_range[0], color='r', alpha=0.5, label=f'{range_mult} Median Absolute Deviation')\n", + " plt.axvline(x=mad_range[1], color='r', alpha=0.5)\n", + " plt.axvline(x=median, color='r', alpha=0.5, linestyle='--', label=f'Median')\n", + " print(f\"Median: {median} {range_mult} Median Absolute Deviation: {mad_range}\")\n", + "\n", + "\n", + " plt.axvline(x=percentile_range[0], color='g', alpha=0.5, label=f'Middle {100-percentile_crop} Percentile')\n", + " plt.axvline(x=percentile_range[1], color='g', alpha=0.5)\n", + " print(f\"Middle {100-percentile_crop} Percentile: {percentile_range}\" )\n", + "\n", + " plt.axvline(x=sigma_range[0], color='y', alpha=0.5, label=f'{range_mult} Sigma')\n", + " plt.axvline(x=sigma_range[1], color='y', alpha=0.5)\n", + " print(f\"{range_mult} Sigma: {sigma_range}\" )\n", + " \n", + " plt.axvline(x=min_max_range[0], color='k', alpha=0.5, label=f'Min/Max')\n", + " plt.axvline(x=min_max_range[1], color='k', alpha=0.5)\n", + " print(f\"Min/Max: {min_max_range}\" )\n", + "\n", + "\n", + " \n", + " \n", + "\n", + " plt.legend(fontsize=5)\n", + " plt.title(\"Data Ranges and Histogram\")\n", + "\n", + " \n", + " img_slice = img[:,img.GetSize()[1]//2,:]\n", + " \n", + " ax = fig.add_subplot(2, 4, 5)\n", + " plt.imshow(sitk.GetArrayViewFromImage(rescale_min_max(img_slice,*mad_range)), cmap='gray')\n", + " ax.axes.xaxis.set_visible(False)\n", + " ax.axes.yaxis.set_visible(False)\n", + " plt.title(f'{range_mult} Median Absolute Deviation', color='r', fontsize=5)\n", + " \n", + " ax = fig.add_subplot(2, 4, 6)\n", + " plt.imshow(sitk.GetArrayViewFromImage(rescale_min_max(img_slice,*percentile_range)), cmap='gray')\n", + " ax.axes.xaxis.set_visible(False)\n", + " ax.axes.yaxis.set_visible(False)\n", + " plt.title(f'Middle {100-percentile_crop} Percentile', color='g', fontsize=5)\n", + " \n", + " ax = fig.add_subplot(2, 4, 7)\n", + " plt.imshow(sitk.GetArrayViewFromImage(rescale_min_max(img_slice,*sigma_range)), cmap='gray')\n", + " ax.axes.xaxis.set_visible(False)\n", + " ax.axes.yaxis.set_visible(False)\n", + " plt.title(f'{range_mult} Sigma', color='y', fontsize=5)\n", + " \n", + " \n", + " ax = fig.add_subplot(2, 4, 8)\n", + " plt.imshow(sitk.GetArrayViewFromImage(rescale_min_max(img_slice, *min_max_range)), cmap='gray')\n", + " ax.axes.xaxis.set_visible(False)\n", + " ax.axes.yaxis.set_visible(False)\n", + " plt.title(f'Data Min/Max', color='k', fontsize=5)\n", + " \n", + " plt.show()\n", + " if fig_filename:\n", + " plt.savefig(fig_filename)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_range_options(img, range_mult=3, percentile_crop=10, bins=1024, fig_filename=\"figure.pdf\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%history\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From a362b09ff66bf1b05c275be2cf9f8161bc185673 Mon Sep 17 00:00:00 2001 From: Bradley Lowekamp Date: Mon, 14 Jun 2021 10:22:47 -0400 Subject: [PATCH 2/4] Update notebook with streaming image to histogram prototype --- notebooks/AutoMinMaxStats.ipynb | 77 +++++++++++++++++++++++++++++++-- pytools/ng/build_histogram.py | 4 ++ 2 files changed, 77 insertions(+), 4 deletions(-) create mode 100644 pytools/ng/build_histogram.py diff --git a/notebooks/AutoMinMaxStats.ipynb b/notebooks/AutoMinMaxStats.ipynb index cf30b90..2b7d5c8 100644 --- a/notebooks/AutoMinMaxStats.ipynb +++ b/notebooks/AutoMinMaxStats.ipynb @@ -11,7 +11,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -20,7 +20,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -33,7 +33,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -177,13 +177,82 @@ "plot_range_options(img, range_mult=3, percentile_crop=10, bins=1024, fig_filename=\"figure.pdf\")" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "a_img = sitk.GetArrayViewFromImage(img).flatten()\n", + "bins = 2**np.iinfo(a_img.dtype).bits\n", + "\n", + "h, b = np.histogram(list(), bins=bins, range=(np.iinfo(a_img.dtype).min, np.iinfo(a_img.dtype).max + 1))\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def stream_build_histogram(filename:str, histogram_bin_edges, extract_axis = 1, density = False):\n", + " reader = sitk.ImageFileReader()\n", + " reader.SetFileName(filename)\n", + " reader.ReadImageInformation()\n", + " \n", + " extract_index = [0]*reader.GetDimension()\n", + "\n", + " extract_size = list(reader.GetSize())\n", + " extract_size[extract_axis] = 0\n", + " reader.SetExtractSize(extract_size)\n", + "\n", + " h = np.zeros(len(histogram_bin_edges)-1, dtype=np.int64)\n", + " \n", + " for i in range(reader.GetSize()[extract_axis]):\n", + " extract_index[extract_axis] = i\n", + " reader.SetExtractIndex(extract_index)\n", + " img = reader.Execute()\n", + " \n", + " # accumulate histogram density/weights\n", + " h += np.histogram( sitk.GetArrayViewFromImage(img).flatten(), bins=histogram_bin_edges, density=density)[0]\n", + " \n", + " return h, np.array(np.histogram_bin_edges)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bin_edges=range(np.iinfo(a_img.dtype).min, np.iinfo(a_img.dtype).max + 2)\n", + "h, b = stream_build_histogram(input_filename, bin_edges)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "len(b2)" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "%history\n" + "len(bin_edges)" ] }, { diff --git a/pytools/ng/build_histogram.py b/pytools/ng/build_histogram.py new file mode 100644 index 0000000..b4b5828 --- /dev/null +++ b/pytools/ng/build_histogram.py @@ -0,0 +1,4 @@ +import numpy as np + +def build_histogram(file_name: str): + pass \ No newline at end of file From 74371c6462807428235df75f034bd1c9a53b1653 Mon Sep 17 00:00:00 2001 From: Bradley Lowekamp Date: Fri, 18 Jun 2021 09:27:48 -0400 Subject: [PATCH 3/4] Updates to batch process files and improvements in figure Fix how data min/max image was displayed to explicily set visualization range. --- notebooks/AutoMinMaxStats.ipynb | 124 ++++++++++++-------------------- 1 file changed, 45 insertions(+), 79 deletions(-) diff --git a/notebooks/AutoMinMaxStats.ipynb b/notebooks/AutoMinMaxStats.ipynb index 2b7d5c8..2d35550 100644 --- a/notebooks/AutoMinMaxStats.ipynb +++ b/notebooks/AutoMinMaxStats.ipynb @@ -11,16 +11,7 @@ }, { "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [], - "source": [ - "input_filename = \"/Users/blowekamp/Desktop/20130412-Vn-poxtomo_8/20130412-Vn-poxtomo_8_full_rec.mrc\"" - ] - }, - { - "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -28,38 +19,19 @@ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.stats import norm\n", + "from pathlib import Path\n", "%matplotlib notebook" ] }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [], - "source": [ - "img = sitk.ReadImage(input_filename)" - ] - }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "stats = sitk.StatisticsImageFilter()\n", - "stats.Execute(img)\n", - "print(stats)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "a_img = sitk.GetArrayViewFromImage(img).flatten()\n", - "bins = 2**np.iinfo(a_img.dtype).bits\n", - "h, b = np.histogram(a_img, bins=bins, range=(np.iinfo(a_img.dtype).min, np.iinfo(a_img.dtype).max + 1))" + "#a_img = sitk.GetArrayViewFromImage(img).flatten()\n", + "#bins = 2**np.iinfo(a_img.dtype).bits\n", + "#h, b = np.histogram(a_img, bins=bins, range=(np.iinfo(a_img.dtype).min, np.iinfo(a_img.dtype).max + 1))" ] }, { @@ -83,7 +55,7 @@ "outputs": [], "source": [ "\n", - "def plot_range_options(img, range_mult=3, percentile_crop=10, bins=1024, fig_filename=None):\n", + "def plot_range_options(img, sigma_mult=3, mad_mult=4, percentile_crop=10, bins=1024, fig_filename=None, title=None):\n", "\n", "\n", " fa_img = sitk.GetArrayViewFromImage(img).flatten()\n", @@ -93,16 +65,18 @@ " sigma = np.sqrt(np.var(fa_img))\n", "\n", "\n", - " mad_range = (median-range_mult*mad, median+range_mult*mad)\n", - " percentile_range = (np.percentile(a_img,percentile_crop*.5), np.percentile(a_img,100-percentile_crop*.5))\n", + " mad_range = (median-mad_mult*mad, median+mad_mult*mad)\n", + " percentile_range = (np.percentile(fa_img,percentile_crop*.5), np.percentile(fa_img,100-percentile_crop*.5))\n", "\n", " \n", - " sigma_range = (mean-range_mult*sigma, \n", - " mean+range_mult*sigma)\n", + " sigma_range = (mean-sigma_mult*sigma, \n", + " mean+sigma_mult*sigma)\n", " \n", " min_max_range = (np.min(fa_img), np.max(fa_img))\n", "\n", " fig = plt.figure(figsize=(8, 8), dpi=240)\n", + " if title:\n", + " fig.suptitle(title, fontsize=8)\n", " ax = fig.add_subplot(2, 4, (1, 4))\n", "\n", " x = np.linspace(min(fa_img), max(fa_img), bins)\n", @@ -110,19 +84,19 @@ "\n", " plt.hist( fa_img, bins=bins, density=True, color='b')\n", "\n", - " plt.axvline(x=mad_range[0], color='r', alpha=0.5, label=f'{range_mult} Median Absolute Deviation')\n", + " plt.axvline(x=mad_range[0], color='r', alpha=0.5, label=f'{mad_mult} Median Absolute Deviation')\n", " plt.axvline(x=mad_range[1], color='r', alpha=0.5)\n", " plt.axvline(x=median, color='r', alpha=0.5, linestyle='--', label=f'Median')\n", - " print(f\"Median: {median} {range_mult} Median Absolute Deviation: {mad_range}\")\n", + " print(f\"Median: {median} {mad_mult} Median Absolute Deviation: {mad_range}\")\n", "\n", "\n", " plt.axvline(x=percentile_range[0], color='g', alpha=0.5, label=f'Middle {100-percentile_crop} Percentile')\n", " plt.axvline(x=percentile_range[1], color='g', alpha=0.5)\n", " print(f\"Middle {100-percentile_crop} Percentile: {percentile_range}\" )\n", "\n", - " plt.axvline(x=sigma_range[0], color='y', alpha=0.5, label=f'{range_mult} Sigma')\n", + " plt.axvline(x=sigma_range[0], color='y', alpha=0.5, label=f'{sigma_mult} Sigma')\n", " plt.axvline(x=sigma_range[1], color='y', alpha=0.5)\n", - " print(f\"{range_mult} Sigma: {sigma_range}\" )\n", + " print(f\"{sigma_mult} Sigma: {sigma_range}\" )\n", " \n", " plt.axvline(x=min_max_range[0], color='k', alpha=0.5, label=f'Min/Max')\n", " plt.axvline(x=min_max_range[1], color='k', alpha=0.5)\n", @@ -135,30 +109,40 @@ " plt.legend(fontsize=5)\n", " plt.title(\"Data Ranges and Histogram\")\n", "\n", - " \n", - " img_slice = img[:,img.GetSize()[1]//2,:]\n", + " if img.GetSize()[0] >= img.GetSize()[2]:\n", + " if img.GetSize()[1] >= img.GetSize()[2]:\n", + " img_slice = img[:,:, img.GetSize()[2]//2]\n", + " else:\n", + " img_slice = img[:,img.GetSize()[1]//2,:]\n", + " else:\n", + " if img.GetSize()[1] >= img.GetSize()[0]:\n", + " img_slice = img[img.GetSize()[0]//2,:,:]\n", + " else:\n", + " img_slice = img[:,img.GetSize()[1]//2,:]\n", + " \n", " \n", " ax = fig.add_subplot(2, 4, 5)\n", - " plt.imshow(sitk.GetArrayViewFromImage(rescale_min_max(img_slice,*mad_range)), cmap='gray')\n", + " plt.imshow(sitk.GetArrayViewFromImage(img_slice), vmin=mad_range[0], vmax=mad_range[1], cmap='gray')\n", " ax.axes.xaxis.set_visible(False)\n", " ax.axes.yaxis.set_visible(False)\n", - " plt.title(f'{range_mult} Median Absolute Deviation', color='r', fontsize=5)\n", + " plt.title(f'{mad_mult} Median Absolute Deviation', color='r', fontsize=5)\n", " \n", " ax = fig.add_subplot(2, 4, 6)\n", - " plt.imshow(sitk.GetArrayViewFromImage(rescale_min_max(img_slice,*percentile_range)), cmap='gray')\n", + " plt.imshow(sitk.GetArrayViewFromImage(img_slice), vmin=percentile_range[0], vmax=percentile_range[1], cmap='gray')\n", " ax.axes.xaxis.set_visible(False)\n", " ax.axes.yaxis.set_visible(False)\n", " plt.title(f'Middle {100-percentile_crop} Percentile', color='g', fontsize=5)\n", " \n", " ax = fig.add_subplot(2, 4, 7)\n", - " plt.imshow(sitk.GetArrayViewFromImage(rescale_min_max(img_slice,*sigma_range)), cmap='gray')\n", + " plt.imshow(sitk.GetArrayViewFromImage(img_slice), vmin=sigma_range[0], vmax=sigma_range[1], cmap='gray')\n", " ax.axes.xaxis.set_visible(False)\n", " ax.axes.yaxis.set_visible(False)\n", - " plt.title(f'{range_mult} Sigma', color='y', fontsize=5)\n", + " plt.title(f'{sigma_mult} Sigma', color='y', fontsize=5)\n", " \n", " \n", " ax = fig.add_subplot(2, 4, 8)\n", - " plt.imshow(sitk.GetArrayViewFromImage(rescale_min_max(img_slice, *min_max_range)), cmap='gray')\n", + " #plt.imshow(sitk.GetArrayViewFromImage(rescale_min_max(img_slice, *min_max_range)), cmap='gray')\n", + " plt.imshow(sitk.GetArrayViewFromImage(img_slice), vmin=min_max_range[0], vmax=min_max_range[1], cmap='gray')\n", " ax.axes.xaxis.set_visible(False)\n", " ax.axes.yaxis.set_visible(False)\n", " plt.title(f'Data Min/Max', color='k', fontsize=5)\n", @@ -171,10 +155,15 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": false + }, "outputs": [], "source": [ - "plot_range_options(img, range_mult=3, percentile_crop=10, bins=1024, fig_filename=\"figure.pdf\")" + "for path in Path.cwd().glob(\"data/*.mrc\"):\n", + " img = sitk.ReadImage(str(path))\n", + " print(path)\n", + " plot_range_options(img, mad_mult=5, percentile_crop=4, bins=1024, fig_filename=f\"{path.with_suffix('')}_fig.pdf\", title=str(path.name))" ] }, { @@ -182,11 +171,13 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "plot_range_options?" + ] }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -236,31 +227,6 @@ "bin_edges=range(np.iinfo(a_img.dtype).min, np.iinfo(a_img.dtype).max + 2)\n", "h, b = stream_build_histogram(input_filename, bin_edges)" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "len(b2)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "len(bin_edges)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { From 13ab315cedc1f8974fa9680b4405a59c00d5503d Mon Sep 17 00:00:00 2001 From: Bradley Lowekamp Date: Tue, 22 Jun 2021 15:58:47 -0400 Subject: [PATCH 4/4] More experimentation with robust statistics --- notebooks/AutoMinMaxStats.ipynb | 47 ++++++++++++++++++++++++++++++--- pytools/utils/__init__.py | 0 2 files changed, 44 insertions(+), 3 deletions(-) create mode 100644 pytools/utils/__init__.py diff --git a/notebooks/AutoMinMaxStats.ipynb b/notebooks/AutoMinMaxStats.ipynb index 2d35550..180b4c3 100644 --- a/notebooks/AutoMinMaxStats.ipynb +++ b/notebooks/AutoMinMaxStats.ipynb @@ -29,9 +29,50 @@ "metadata": {}, "outputs": [], "source": [ - "#a_img = sitk.GetArrayViewFromImage(img).flatten()\n", - "#bins = 2**np.iinfo(a_img.dtype).bits\n", - "#h, b = np.histogram(a_img, bins=bins, range=(np.iinfo(a_img.dtype).min, np.iinfo(a_img.dtype).max + 1))" + "fname = \"data/ave_2013-1220-dA30_5-BSC-1_7.mrc\"\n", + "img = sitk.ReadImage(fname)\n", + "a_img = sitk.GetArrayViewFromImage(img).flatten()\n", + "bins = 2**np.iinfo(a_img.dtype).bits\n", + "h, b = np.histogram(a_img, bins=bins, range=(np.iinfo(a_img.dtype).min, np.iinfo(a_img.dtype).max + 1))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def histogram_robust_stats(hist, bin_edges):\n", + "\n", + " assert (len(hist) + 1 == len(bin_edges))\n", + " results = {}\n", + "\n", + " mids = 0.5 * (bin_edges[1:] + bin_edges[:-1])\n", + " cs = np.cumsum(hist)\n", + " median_idx = np.searchsorted(cs, cs[-1] / 2.0)\n", + " results[\"median\"] = mids[median_idx]\n", + " \n", + " mad_mids, mad_h = zip(*sorted(zip(np.abs(mids - results[\"median\"] ), hist)))\n", + " mad_cs = np.cumsum(mad_h)\n", + " mad_idx = np.searchsorted(mad_cs, mad_cs[-1] / 2.0)\n", + "\n", + " \n", + " results[\"mad\"] = mad_mids[mad_idx]\n", + " \n", + " return results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(histogram_robust_stats(h,b))\n", + "\n", + "median = np.median(a_img)\n", + "mad = np.median(np.abs(a_img-median))\n", + "print(f\"median: {median}\\nMAD: {mad}\")" ] }, { diff --git a/pytools/utils/__init__.py b/pytools/utils/__init__.py new file mode 100644 index 0000000..e69de29