diff --git a/ITKIOScanco_UnitConversion.ipynb b/ITKIOScanco_UnitConversion.ipynb new file mode 100644 index 0000000..00e4880 --- /dev/null +++ b/ITKIOScanco_UnitConversion.ipynb @@ -0,0 +1,247 @@ +{ + "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 +}