Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add Python Jupyter notebook with prototype of visualization range #10

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Add Python Jupyter notebook with prototype of visualization range #10

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
14ccb7e
Add Python Jupyter notebook with prototype of visualization range
blowekamp Jun 9, 2021
a362b09
Update notebook with streaming image to histogram prototype
blowekamp Jun 14, 2021
74371c6
Updates to batch process files and improvements in figure
blowekamp Jun 18, 2021
13ab315
More experimentation with robust statistics
blowekamp Jun 22, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
294 changes: 294 additions & 0 deletions notebooks/AutoMinMaxStats.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
}
4 changes: 4 additions & 0 deletions pytools/ng/build_histogram.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
import numpy as np

def build_histogram(file_name: str):
pass
Empty file added pytools/utils/__init__.py
View file
Open in desktop
Empty file.