diff --git a/notebooks/AutoMinMaxStats.ipynb b/notebooks/AutoMinMaxStats.ipynb new file mode 100644 index 0000000..180b4c3 --- /dev/null +++ b/notebooks/AutoMinMaxStats.ipynb @@ -0,0 +1,294 @@ +{ + "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": [ + "import SimpleITK as sitk\n", + "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": null, + "metadata": {}, + "outputs": [], + "source": [ + "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}\")" + ] + }, + { + "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, 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", + " 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-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-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", + " 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'{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} {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'{sigma_mult} Sigma')\n", + " plt.axvline(x=sigma_range[1], color='y', alpha=0.5)\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", + " 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", + " 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(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'{mad_mult} Median Absolute Deviation', color='r', fontsize=5)\n", + " \n", + " ax = fig.add_subplot(2, 4, 6)\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(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'{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(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", + " \n", + " plt.show()\n", + " if fig_filename:\n", + " plt.savefig(fig_filename)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "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))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_range_options?" + ] + }, + { + "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", + "\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)" + ] + } + ], + "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 +} 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 diff --git a/pytools/utils/__init__.py b/pytools/utils/__init__.py new file mode 100644 index 0000000..e69de29