Merge pull request #12 from SpectraCollab/mkuczyns/cleanup

Mkuczyns/cleanup

.github/workflows/draft-pdf.yml

@@ -0,0 +1,19 @@
+on: [push]
+
+jobs:
+  paper:
+    runs-on: ubuntu-latest
+    name: ORMIR_XCT Paper Draft
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v3
+      - name: Build draft PDF
+        uses: openjournals/openjournals-draft-action@master
+        with:
+          journal: joss
+          paper-path: manuscript/paper.md
+      - name: Upload
+        uses: actions/upload-artifact@v1
+        with:
+          name: ORMIR_XCT Paper Draft
+          path: manuscript/paper.pdf

README.md

@@ -1,2 +1,9 @@
 # ORMIR_2022
-This project was partially developed during the Jupyter Community Workshop “Building the Jupyter Community in Musculoskeletal Imaging Research”.
+ORMIR_XCT is a Python package for processing high resolution peripheral computed tomography (HR-pQCT) scans. Development of this project began during the “Building the Jupyter Community in Musculoskeletal Imaging Research” workshop hosted by the Open and Reproducible Musculoskeletal Imaging Research (ORMIR) group.
+
+## Installation
+### Step 1: Install the ORMIR_XCT Anaconda Environment:
+`conda env create -f environment.yml`
+
+### Step 2: Install the Package:
+`pip install -e .`

examples/finger_joint_comparison.ipynb

@@ -0,0 +1,332 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "78149392",
+   "metadata": {},
+   "source": [
+    "# Finger Joint Comparison\n",
+    "\n",
+    "This notebook compares results from the IPL joint space width (JSW) analysis workflow to the results from the ORMIR_XCT package. To provide a more direct comparison, the same joint segmentation was used as input to the JSW analysis (generated using IPL). ORMIR_XCT results are shown with all combonations of oversampling and skeletonization.\n",
+    "\n",
+    "## IPL JSW Results\n",
+    "\n",
+    "The following distance transform image was generated in IPL:\n",
+    "<img src=\"images/IPL_JSW_DT_IMAGE.png\" width=\"300\" height=\"200\">\n",
+    "\n",
+    "The following JSW results were obtained from the above distance transform image:\n",
+    "\n",
+    "| JSW Parameter | Value |\n",
+    "|:--------:|:--------:|\n",
+    "|  JS Volume (mm^3)  |  78.2068   |\n",
+    "|  Mean JS Width (mm)  |  1.8472   |\n",
+    "|  JS Max (mm)   |  2.6101   |\n",
+    "|  JS Min (mm)   |  1.3354   |\n",
+    "|  JS Asymmetry (JS Max / JS Min)   |  1.9545   |"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ebee8652",
+   "metadata": {},
+   "source": [
+    "# Step 1: Imports\n",
+    "\n",
+    "Import modules/packages and set the input image path. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "93e98f00",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import SimpleITK as sitk\n",
+    "\n",
+    "from ormir_xct.joint_space_analysis.jsw_morphometry import jsw_pad, jsw_dilate, jsw_erode, jsw_parameters"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "7b4af5a3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "joint_seg_path = os.path.join(\"images\", \"IPL_vs_ORMIR_XCT_JOINT.nii\")\n",
+    "output_path = \"images\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "79816632",
+   "metadata": {},
+   "source": [
+    "# Step 2: Run the JSW Workflow\n",
+    "\n",
+    "1. **Padding:** Pad the binary joint image with black space (zeros) to ensure the outside space is greater than the inside (joint) space.\n",
+    "2. **Dilation:** Next we dilate the binary image by a fixed constant (see `jsw_morphometry.py`) that is taken from the original IPL script. Here we use the SimpleITK Ball structural element for dilation which should be similar to the Euclidean metric used in IPL. Once dilation is complete, remove any islands and fill any holes in the dilated mask.\n",
+    "3. **Erosion:** Next erode the image and set the image's value to 30. Then, add the erosion and dilation images together. Areas with overlap between images will have a value of 90, and the joint space will have a value of 30. We can then use a binary threshold to extract the joint space. The joint space mask is then dilated to ensure correct distances are measured at the edges of the joint space.\n",
+    "4. **JSW Parameters:** Now we will compute JS volume, JSW mean, standard deviation, minimum, maximum, and JS asymmetry (JSW.Max / JSW.Min). Results are output to a CSV file."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "f460183e",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Padding image...\n",
+      "Dilating image...\n",
+      "Eroding image...\n"
+     ]
+    }
+   ],
+   "source": [
+    "filename = os.path.basename(joint_seg_path)\n",
+    "basename = os.path.splitext(filename)[0]\n",
+    "img = sitk.ReadImage(joint_seg_path, sitk.sitkUInt8)\n",
+    "\n",
+    "# 1. Padding\n",
+    "print(\"Padding image...\")\n",
+    "pad_image = jsw_pad(img)\n",
+    "\n",
+    "# 2. Dilation\n",
+    "print(\"Dilating image...\")\n",
+    "dilated_image = jsw_dilate(pad_image)\n",
+    "sitk.WriteImage(dilated_image, os.path.join(output_path, str(basename) + \"_DILATE.nii\"))\n",
+    "\n",
+    "# 3. Erosion\n",
+    "print(\"Eroding image...\")\n",
+    "eroded_image, js_mask, dilated_js_mask = jsw_erode(dilated_image, pad_image)\n",
+    "sitk.WriteImage(eroded_image, os.path.join(output_path, str(basename) + \"_ERODE.nii\"))\n",
+    "sitk.WriteImage(js_mask, os.path.join(output_path, str(basename) + \"_JS_MASK.nii\"))\n",
+    "sitk.WriteImage(dilated_js_mask, os.path.join(output_path, str(basename) + \"_DILATED_JS_MASK.nii\"))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "44fc2d09",
+   "metadata": {},
+   "source": [
+    "# Step 3: Compute the JSW Parameters\n",
+    "\n",
+    "All combinations of oversampling and skeletonization will be run through the distance transform algorithm to see which set of parameters best matches results from IPL."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "bfbb96f0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# 4. JSW Parameters\n",
+    "\n",
+    "# Oversampling=False, Skeletonization=False\n",
+    "dt_img1, jsw_params1 = jsw_parameters(pad_image, dilated_js_mask, basename,\n",
+    "                                      output_path, js_mask, pad_image.GetSpacing()[0], \n",
+    "                                     oversamp=False, skel=False)\n",
+    "sitk.WriteImage(dt_img1, os.path.join(output_path, str(basename) + \"_DT1.nii\"))\n",
+    "\n",
+    "# Oversampling=True, Skeletonization=False\n",
+    "dt_img2, jsw_params2 = jsw_parameters(pad_image, dilated_js_mask, basename,\n",
+    "                                      output_path, js_mask, pad_image.GetSpacing()[0], \n",
+    "                                     oversamp=True, skel=False)\n",
+    "sitk.WriteImage(dt_img1, os.path.join(output_path, str(basename) + \"_DT2.nii\"))\n",
+    "\n",
+    "# Oversampling=False, Skeletonization=True\n",
+    "dt_img3, jsw_params3 = jsw_parameters(pad_image, dilated_js_mask, basename,\n",
+    "                                      output_path, js_mask, pad_image.GetSpacing()[0], \n",
+    "                                     oversamp=False, skel=True)\n",
+    "sitk.WriteImage(dt_img1, os.path.join(output_path, str(basename) + \"_DT3.nii\"))\n",
+    "\n",
+    "# Oversampling=True, Skeletonization=True\n",
+    "dt_img4, jsw_params4 = jsw_parameters(pad_image, dilated_js_mask, basename,\n",
+    "                                      output_path, js_mask, pad_image.GetSpacing()[0], \n",
+    "                                     oversamp=True, skel=True)\n",
+    "sitk.WriteImage(dt_img1, os.path.join(output_path, str(basename) + \"_DT4.nii\"))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "68b9841c",
+   "metadata": {},
+   "source": [
+    "# Step 4: Display Results\n",
+    "\n",
+    "Use Pandas to print a table of the results we generated."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "id": "f3789891",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<style type=\"text/css\">\n",
+       "</style>\n",
+       "<table id=\"T_e6376\">\n",
+       "  <thead>\n",
+       "    <tr>\n",
+       "      <th class=\"blank level0\" >&nbsp;</th>\n",
+       "      <th id=\"T_e6376_level0_col0\" class=\"col_heading level0 col0\" >ipl_results</th>\n",
+       "      <th id=\"T_e6376_level0_col1\" class=\"col_heading level0 col1\" >oversamp=False, skel=False</th>\n",
+       "      <th id=\"T_e6376_level0_col2\" class=\"col_heading level0 col2\" >oversamp=True, skel=False</th>\n",
+       "      <th id=\"T_e6376_level0_col3\" class=\"col_heading level0 col3\" >oversamp=False, skel=True</th>\n",
+       "      <th id=\"T_e6376_level0_col4\" class=\"col_heading level0 col4\" >oversamp=True, skel=True</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th id=\"T_e6376_level0_row0\" class=\"row_heading level0 row0\" >JS Volume (mm^3)</th>\n",
+       "      <td id=\"T_e6376_row0_col0\" class=\"data row0 col0\" >78.21</td>\n",
+       "      <td id=\"T_e6376_row0_col1\" class=\"data row0 col1\" >80.47</td>\n",
+       "      <td id=\"T_e6376_row0_col2\" class=\"data row0 col2\" >80.47</td>\n",
+       "      <td id=\"T_e6376_row0_col3\" class=\"data row0 col3\" >80.47</td>\n",
+       "      <td id=\"T_e6376_row0_col4\" class=\"data row0 col4\" >80.47</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th id=\"T_e6376_level0_row1\" class=\"row_heading level0 row1\" >Mean JS Width (mm)</th>\n",
+       "      <td id=\"T_e6376_row1_col0\" class=\"data row1 col0\" >1.85</td>\n",
+       "      <td id=\"T_e6376_row1_col1\" class=\"data row1 col1\" >1.92</td>\n",
+       "      <td id=\"T_e6376_row1_col2\" class=\"data row1 col2\" >1.92</td>\n",
+       "      <td id=\"T_e6376_row1_col3\" class=\"data row1 col3\" >1.95</td>\n",
+       "      <td id=\"T_e6376_row1_col4\" class=\"data row1 col4\" >1.73</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th id=\"T_e6376_level0_row2\" class=\"row_heading level0 row2\" >JS Max (mm)</th>\n",
+       "      <td id=\"T_e6376_row2_col0\" class=\"data row2 col0\" >2.61</td>\n",
+       "      <td id=\"T_e6376_row2_col1\" class=\"data row2 col1\" >2.44</td>\n",
+       "      <td id=\"T_e6376_row2_col2\" class=\"data row2 col2\" >2.35</td>\n",
+       "      <td id=\"T_e6376_row2_col3\" class=\"data row2 col3\" >2.41</td>\n",
+       "      <td id=\"T_e6376_row2_col4\" class=\"data row2 col4\" >2.31</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th id=\"T_e6376_level0_row3\" class=\"row_heading level0 row3\" >JS Min (mm)</th>\n",
+       "      <td id=\"T_e6376_row3_col0\" class=\"data row3 col0\" >1.34</td>\n",
+       "      <td id=\"T_e6376_row3_col1\" class=\"data row3 col1\" >1.32</td>\n",
+       "      <td id=\"T_e6376_row3_col2\" class=\"data row3 col2\" >1.35</td>\n",
+       "      <td id=\"T_e6376_row3_col3\" class=\"data row3 col3\" >0.24</td>\n",
+       "      <td id=\"T_e6376_row3_col4\" class=\"data row3 col4\" >0.19</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th id=\"T_e6376_level0_row4\" class=\"row_heading level0 row4\" >JS Asymmetry (JS Max / JS Min)</th>\n",
+       "      <td id=\"T_e6376_row4_col0\" class=\"data row4 col0\" >1.95</td>\n",
+       "      <td id=\"T_e6376_row4_col1\" class=\"data row4 col1\" >1.85</td>\n",
+       "      <td id=\"T_e6376_row4_col2\" class=\"data row4 col2\" >1.75</td>\n",
+       "      <td id=\"T_e6376_row4_col3\" class=\"data row4 col3\" >9.94</td>\n",
+       "      <td id=\"T_e6376_row4_col4\" class=\"data row4 col4\" >12.05</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n"
+      ],
+      "text/plain": [
+       "<pandas.io.formats.style.Styler at 0x208e7b71340>"
+      ]
+     },
+     "execution_count": 11,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# Get the JS Volume\n",
+    "shape_stats = sitk.LabelShapeStatisticsImageFilter()\n",
+    "shape_stats.ComputeOrientedBoundingBoxOn()\n",
+    "shape_stats.Execute(js_mask)\n",
+    "\n",
+    "stats_list = [\n",
+    "    (\n",
+    "        shape_stats.GetPhysicalSize(i),\n",
+    "        shape_stats.GetNumberOfPixels(i),\n",
+    "    )\n",
+    "    for i in shape_stats.GetLabels()\n",
+    "]\n",
+    "\n",
+    "jsv = stats_list[0][0]\n",
+    "\n",
+    "# Oversampling=False, Skeletonization=False\n",
+    "js_mean1 = jsw_params1[0][3]\n",
+    "js_max1 = jsw_params1[0][7]\n",
+    "js_min1 = jsw_params1[0][5]\n",
+    "js_as1 = js_max1 / js_min1\n",
+    "\n",
+    "# Oversampling=True, Skeletonization=False\n",
+    "js_mean2 = jsw_params2[0][3]\n",
+    "js_max2 = jsw_params2[0][7]\n",
+    "js_min2 = jsw_params2[0][5]\n",
+    "js_as2 = js_max2 / js_min2\n",
+    "\n",
+    "# Oversampling=False, Skeletonization=True\n",
+    "js_mean3 = jsw_params3[0][3]\n",
+    "js_max3 = jsw_params3[0][7]\n",
+    "js_min3 = jsw_params3[0][5]\n",
+    "js_as3 = js_max3 / js_min3\n",
+    "\n",
+    "# Oversampling=True, Skeletonization=True\n",
+    "js_mean4 = jsw_params4[0][3]\n",
+    "js_max4 = jsw_params4[0][7]\n",
+    "js_min4 = jsw_params4[0][5]\n",
+    "js_as4 = js_max4 / js_min4\n",
+    "\n",
+    "\n",
+    "# Display a table with the IPL and ORMIR_XCT calculated JS parameters\n",
+    "data = {\"ipl_results\":  [78.2068, 1.8472, 2.6101, 1.3354, 1.9545],\n",
+    "        \"oversamp=False, skel=False\": [jsv, js_mean1, js_max1, js_min1, js_as1],\n",
+    "        \"oversamp=True, skel=False\": [jsv, js_mean1, js_max2, js_min2, js_as2],\n",
+    "        \"oversamp=False, skel=True\": [jsv, js_mean2, js_max3, js_min3, js_as3],\n",
+    "        \"oversamp=True, skel=True\": [jsv, js_mean3, js_max4, js_min4, js_as4],\n",
+    "        }\n",
+    "df = pd.DataFrame(data, index=[\"JS Volume (mm^3)\", \"Mean JS Width (mm)\", \"JS Max (mm)\", \n",
+    "                               \"JS Min (mm)\", \"JS Asymmetry (JS Max / JS Min)\"])\n",
+    "\n",
+    "# Set table default options\n",
+    "pd.set_option('display.max_colwidth', None)\n",
+    "pd.set_option('colheader_justify', 'center')\n",
+    "\n",
+    "df.style.format(precision=2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "21bce143",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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": 5
+}