Jmcox9 unit conversion #74
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,247 @@ | ||
| { | ||
|
There was a problem hiding this comment. Line #9. del all_AIM_filenames, all_ISQ_filenames We can remove the Reply via ReviewNB There was a problem hiding this comment. Line #1. filename = 'C0005757_SCAP.AIM' Can you download the file with the Reply via ReviewNB There was a problem hiding this comment. For now I've implemented this using the example data referred to here. I'm not entirely sure how this is supposed to function in a general sense, since the URL provided in the example may not be the location where these images will be found. The hash value seems to be an important argument for pooch.retrieve as well, but it's not clear to me how this will be obtained. There was a problem hiding this comment. Line #4. print(metadata) Can we save this cell's output in the notebook so it can be browsed on GitHub? Reply via ReviewNB There was a problem hiding this comment. I'm approaching this as saving the metadata to a .txt file. If you had something else in mind, please let me know. There was a problem hiding this comment. Line #6. mu_threshHigh = 1.0 per Python naming conversions, snake_case is used for variable names. So Reply via ReviewNB There was a problem hiding this comment. Line #7. reader.SetFileName(filename) We want to use the more python Reply via ReviewNB There was a problem hiding this comment. Line #17. BWoutname = filename[0:-4] + '_BWmask.nrrd' from pathlib import Path Path(filename)
-> remove the extension with this. Reply via ReviewNB There was a problem hiding this comment. There was a problem hiding this comment. Line #5. print(metadata) can we save the output here into the notebook so it can be viewed on GitHub? Reply via ReviewNB There was a problem hiding this comment. There was a problem hiding this comment. Line #29. HU_values = [HU_pixel_low, HU_pixel_high] Please print out all the relevant scalar values here at the end of the cell and save the output in the cell so people browsing the notebook online can see examples of typical values. Reply via ReviewNB |
||
| "cells": [ | ||
| { | ||
| "attachments": {}, | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "This workbook loads and converts a SCANCO .ISQ or .AIM file in the directory to several \n", | ||
| ".nrrd formatted images with the following units: density (g/cc), Hounsfield Units, and a\n", | ||
| "binary thresholded image mask \n", | ||
| "\n", | ||
| "Created by: Jason Cox, 04/07/2023\n", | ||
| "\n", | ||
| "Citations: McCormick, Matthew, et al. \"ITK: Enabling Reproducible \n", | ||
| "Research and Open Science.” Frontiers in Neuroinformatics, vol. 8, \n", | ||
| "Frontiers Media SA, 2014, doi:10.3389/fninf.2014.00013." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# Install non-default packages\n", | ||
| "import sys\n", | ||
| "!{sys.executable} -m pip install itk pooch" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 4, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# Import required packages\n", | ||
| "from pathlib import Path\n", | ||
| "import re\n", | ||
| "import itk\n", | ||
| "import pooch\n", | ||
| "import glob" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 2, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# ISQ and AIM metadata returned by dict (cell 3) contain the same fields. Therefore, this program could be modified to loop over all files\n", | ||
| "# in the directory, if that is the desired implementation. The variable all_filenames is unused in the current version, but can be used in a\n", | ||
| "# modified version that implements a loop.\n", | ||
| "\n", | ||
| "all_AIM_filenames = glob.glob(\"*.AIM\")\n", | ||
| "all_ISQ_filenames = glob.glob(\"*.ISQ\")\n", | ||
| "all_filenames = all_AIM_filenames + all_ISQ_filenames\n" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 15, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# This cell originally operated on the image C0005757_SCAP.AIM, though for the purposes of employing pooch, a sample file stored online at the provided URL is now used.\n", | ||
| "# This section will need to be modified to handle general implementation.\n", | ||
| "\n", | ||
| "# filename = 'C0005757_SCAP.AIM'\n", | ||
| "filename = 'C0004255.ISQ'\n", | ||
| "file_url = 'https://data.kitware.com/api/v1/file/591e56178d777f16d01e0d20/download'\n", | ||
| "file_sha256 = 'c2a3750c75826cb077d92093d43976cc0350198b55edecd681265eebabfb438b'\n", | ||
| "file_path = pooch.retrieve(file_url, file_sha256, fname=filename, progressbar=False)\n", | ||
| "image = itk.imread(file_path)\n", | ||
| "metadata = dict(image)\n", | ||
| "HU_pixels = itk.array_view_from_image(image)\n", | ||
| "# del image\n", | ||
| "# HU_pixels = np.array(HU_pixels) # Array arithmetic below seems to work without this line, but certain operations were not tested as they require too much memory" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# As discussed below, some required constants are missing from metadata. Compare the output of this cell with the file header at:\n", | ||
| "# https://docs.google.com/document/d/1L2GCnHa4jmK87p0pGI9MNN10MWRMwrk2WqPbPeklP_I/edit\n", | ||
| "# This section does, however, collect some metadata and write it to a .txt file.\n", | ||
| "\n", | ||
| "print(metadata)\n", | ||
| "metadata_str = str(metadata)\n", | ||
| "\n", | ||
| "metadata_outname = Path(filename).stem + \"_metadata.txt\"\n", | ||
| "with open(metadata_outname, \"w\") as file:\n", | ||
| " file.write(metadata_str)" | ||
| ] | ||
| }, | ||
| { | ||
| "attachments": {}, | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "The next cell provides a template file header from C0005757_SCAP.AIM, but this needs to be read from the .AIM or .ISQ.\n", | ||
| "\n", | ||
| "Then, the following cell identifies numeric constants in the header necessary for converting the image in HU (units given by itk.imread)\n", | ||
| "to density, or for thresholding/masking the image.\n", | ||
| "\n", | ||
| "The only constants not readable from the header which must be provided by the user are the per-milli values for low and high thresholds. In basic terms, the low/high threshold pixel values are computed as the respective low/high per-milli value multiplied by the highest pixel value in the image (a quanitity that is given by the header)." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 17, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "header = '{Global Density: intercept=-1.96600006e+02, Global Calib. default unit type=2 (Density), Global Position Slice 1 [um]=47050, Global HU: mu water=0.51400, Global Reconstruction-Alg.=3, Global Parameter units=[1/cm], Global Average data value=322.51215, Global Minimum data value=-2891.00000, Global Created by=ISQ_TO_AIM (IPL), Global Patient Name=JC-NBPI-Rat-May17, Global Energy [V]=70000, Global Scanner ID=2135, Global Standard deviation=0.32443, Global Time=29-JUN-2022 12:39:11.75, Global Orig-ISQ-Dim-p=1792 1792 1781, Global Original Creation-Date=20-AUG-2018 23:18:59.99, Global Parameter name=Linear Attenuation, Global Orig-ISQ-Dim-um=17918 17918 17808, Global Scanner type=10, Global Scan Distance [um]=20480, Global Index Patient=136, Global Index Measurement=7866, Global Minimum value=-0.70581, Global Maximum data value=32767.00000, Global Intensity [uA]=114, Global Maximum value=7.99976, Global Density: slope=3.55519989e+02, Global Mu_Scaling=4096, Global Reference line [um]=47050, Global No. samples=2048, Global Default-Eval=15, Global Site=4, Global Average value=0.07874, Global Calibration Data=Unknown. Backconverted from AIM., Global Standard data deviation=1328.87610, Global No. projections per 180=1000, Global Original file=DK0:[MICROCT.DATA.00000136.00007866]B0005757.ISQ, Global Scaled by factor=4096, Global Angle-Offset [mdeg]=0, Global Integration time [us]=800000, Global Density: unit=mg HA/ccm}'" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 18, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# Gather numeric constants needed for unit conversion. Currently, the file metadata returned with dict does not contain all of the\n", | ||
| "# necessary constants. For an example of a full file header with highlighted variables of interest, see the google doc \n", | ||
| "# \"SCANCO micro-CT Units Conversion\" at: https://docs.google.com/document/d/1L2GCnHa4jmK87p0pGI9MNN10MWRMwrk2WqPbPeklP_I/edit\n", | ||
| "\n", | ||
| "# Thresholds must be entered by the user, they are not contained in the file header\n", | ||
| "mu_threshold_low = 0.250\n", | ||
| "mu_threshold_high = 1.0\n", | ||
| "\n", | ||
| "scanco_mu = float(re.search('Global Scaled by factor=(.*?),', header).group(1))\n", | ||
| "# scanco_mu = metadata['MuScaling']\n", | ||
| "max_val = float(re.search('Global Maximum data value=(.*?),', header).group(1))\n", | ||
| "# max_val = metadata['**********']\n", | ||
| "density_slope = float(re.search('slope=(.*?),', header).group(1))\n", | ||
| "# density_slope = metadata['**********']\n", | ||
| "density_intercept = float(re.search('intercept=(.*?),', header).group(1))\n", | ||
| "# density_intercept = metadata['**********']\n", | ||
| "water_mu = float(re.search('water=(.*?),', header).group(1))\n", | ||
| "# water_mu = metadata['MuWater']\n", | ||
| "\n", | ||
| "# del header # when the code base is modified such that metadata gives all relevant constants, this line and others that reference \"header\" can be deleted\n", | ||
| "# del metadata\n", | ||
| "\n", | ||
| "mu_low = mu_threshold_low * max_val\n", | ||
| "mu_pixel_low = mu_low/scanco_mu\n", | ||
| "mu_high = mu_threshold_high * max_val\n", | ||
| "mu_pixel_high = mu_high/scanco_mu\n", | ||
| "HU_pixel_low = (-1000) + mu_pixel_low * (1000/water_mu)\n", | ||
| "HU_pixel_high = (-1000) + mu_pixel_high * (1000/water_mu)\n", | ||
| "HU_values = [HU_pixel_low, HU_pixel_high]" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 25, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# Create a binarized version of the input image by applying low and high thresholds. Threshold values are accepted as int type variables, so some rounding is necessary.\n", | ||
| "lower_threshold = int(round(HU_values[0]))\n", | ||
| "upper_threshold = int(round(HU_values[1]))\n", | ||
| "inside_value = 1\n", | ||
| "outside_value = 0\n", | ||
| "\n", | ||
| "BW_image = itk.binary_threshold_image_filter(image,lower_threshold=lower_threshold,upper_threshold=upper_threshold,inside_value=inside_value,outside_value=outside_value)\n", | ||
| "\n", | ||
| "BW_outname = Path(filename).stem + \"_BWmask.nrrd\"\n", | ||
| "\n", | ||
| "itk.imwrite(BW_image, BW_outname)\n", | ||
| "\n", | ||
| "# The following can be enabled to free up memory\n", | ||
| "# del image\n", | ||
| "# del BW_image" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 26, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# Create unit-converted outputs. MUpixels is not generally desired as an output due to representing non-standard units. However, creation of MUpixels\n", | ||
| "# can be skipped, as RHOpixels can be created by combining the functions needed to convert HU to MU and MU to RHO.\n", | ||
| "\n", | ||
| "# MUpixels = (HUpixels+1000)/(1000/water_mu)\n", | ||
| "# RHOpixels = (MUpixels*densitySlope) + densityIntercept\n", | ||
| "RHO_pixels = ((((HU_pixels+1000)/(1000/water_mu))*density_slope) + density_intercept)/1000 # Dividing by 1000 at the end changes units from mg/cc to g/cc, which is more commonly published\n" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 29, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# Create .nrrd file output with units of density (mg/cc or g/cc)\n", | ||
| "\n", | ||
| "RHO_out = itk.image_view_from_array(RHO_pixels)\n", | ||
| "for k, v in metadata.items():\n", | ||
| " RHO_out[k] = v\n", | ||
| "RHO_outname = Path(filename).stem + \"_RHO.nrrd\"\n", | ||
| "itk.imwrite(RHO_out, RHO_outname)" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 30, | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# Create .nrrd file output with Hounsfield Units (equivalent to loading the file in 3D Slicer)\n", | ||
| "\n", | ||
| "HU_out = itk.image_view_from_array(HU_pixels)\n", | ||
| "for k, v in metadata.items():\n", | ||
| " HU_out[k] = v\n", | ||
| "HU_outname = Path(filename).stem + \"_HU.nrrd\"\n", | ||
| "itk.imwrite(HU_out, HU_outname)" | ||
| ] | ||
| } | ||
| ], | ||
| "metadata": { | ||
| "kernelspec": { | ||
| "display_name": "kitfisto", | ||
| "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.11.2" | ||
| }, | ||
| "orig_nbformat": 4 | ||
| }, | ||
| "nbformat": 4, | ||
| "nbformat_minor": 2 | ||
| } | ||
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Before importing, can we first have a cell that installs all dependencies so the notebook is more reproducible? See, e.g. https://github.com/KitwareMedical/ITKIOScanco/blob/master/examples/ReadISQ.ipynb
Reply via ReviewNB