diff --git a/README.md b/README.md index 428e4ad..7339790 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,12 @@ -# ORMIR_2022 -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. +# ORMIR_XCT + +**By:** [Michael T. Kuczynski](https://www.linkedin.com/in/mkuczyns/), [Nathan J. Neeteson](https://www.linkedin.com/in/nathan-neeteson/), [Kathryn S. Stok](https://www.linkedin.com/in/kstok/), [Andrew J. Burghardt](https://www.linkedin.com/in/aburghardt/), [Michelle A. Espinosa Hernandez](https://www.linkedin.com/in/michelleaespinosah/), [Jared Vicory](https://www.kitware.com/jared-vicory/), [Justin J. Tse](https://www.linkedin.com/in/justin-j-tse/), [Pholpat Durongbhan](https://www.linkedin.com/in/pholpatd/), [Serena Bonaretti](https://sbonaretti.github.io/), [Andy Kin On Wong](https://www.linkedin.com/in/andy-kin-on-wong-76408859/), [Steven K. Boyd](https://mccaig.ucalgary.ca/boyd), [Sarah L. Manske](https://www.linkedin.com/in/sarah-manske-b5402b41/), 2023 +**Version:** 1.0 + +- 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: @@ -15,4 +22,27 @@ If using an Apple M1, M2, or M3 processor, run the following command instead: `pip install -e .` ### Step 4: Run Scripts: -The modules in the `ormir_xct` directory can now be run. Examples for each module are provided in the `examples` directory. \ No newline at end of file +The modules in the `ormir_xct` directory can now be run. Examples for each module are provided in the `examples` directory. + +--- + +## Example Jupyter Notebooks +- Example Jupyter Notebooks demonstrating the major functionality of the ORMIR_XCT package are provided in the *[examples](https://github.com/SpectraCollab/ORMIR_XCT/tree/main/examples)* directory. + +--- + +## Ways to Contribute +### Reporting Bugs +- Bugs can be reported by creating a new GitHub issue in this repository. For each bug, please provide details on how to reproduce the bug and the specific error message (if possible). + +--- + +### Contributing New Features +- To add a new feature, expand existing functionality, add documentation, or other contributions, please submit a new GitHub issue outlining your contribution in detail. +- When submitting a new pull request, ensure you outline what you have changed and why it is necessary to make this change. + +--- + +## Citation +When using the ORMIR_XCT package, please use the following citation: +- *Kuczynski, M.T., et al. "ORMIR_XCT: A Python package for high resolution peripheral quantitative computed tomography image processing." arXiv preprint arXiv:2309.04602 (2023).* \ No newline at end of file diff --git a/environment.yml b/environment.yml index 34acd1f..23af62e 100644 --- a/environment.yml +++ b/environment.yml @@ -21,7 +21,7 @@ dependencies: - itk-core==5.3.0 - itk-filtering==5.3.0 - itk-io==5.3.0 - - itk-ioscanco==0.9.2 + - itk-ioscanco==0.10.0 - itk-numerics==5.3.0 - itk-registration==5.3.0 - itk-segmentation==5.3.0 diff --git a/examples/Automatic_Contour_Example.ipynb b/examples/Automatic_Contour_Example.ipynb index a44ee8b..86db471 100644 --- a/examples/Automatic_Contour_Example.ipynb +++ b/examples/Automatic_Contour_Example.ipynb @@ -7,9 +7,37 @@ "source": [ "# Automatic Contour Example\n", "\n", - "This Jupyter Notebook provides an example of running the automatic periosteal contouring script from the ORMIR_XCT package. An example periosteal contour generated using the standard IPL workflow is used as comparison. The DICE coefficient, Jaccard index, and Hausdorff distance are used to compare between the two contour methods.\n", + "- **By:** [Michael T. Kuczynski](https://www.linkedin.com/in/mkuczyns/), 2024 \n", + "- **License:** CC-BY \n", + "- **How to cite:** Cite the ORMIR_XCT publication: *Kuczynski, M.T., et al. \"ORMIR_XCT: A Python package for high resolution peripheral quantitative computed tomography image processing.\" arXiv preprint arXiv:2309.04602 (2023).*" + ] + }, + { + "cell_type": "markdown", + "id": "3a3b8d0f", + "metadata": {}, + "source": [ + "---\n", + "# Aims\n", "\n", - "The ORMIR_XCT automatic contouring script performs image segmentation on HR-pQCT images of joints to separate bone from other tissues. The script outputs the proximal, distal, and full joint mask for each image." + "- This Jupyter Notebook provides an example of running the automatic periosteal contouring script from the ORMIR_XCT package. \n", + "- An example periosteal contour generated using the standard IPL workflow is used as comparison. \n", + "- The DICE coefficient, Jaccard index, and Hausdorff distance are used to compare between the two contour methods.\n", + "- The ORMIR_XCT automatic contouring script performs image segmentation on HR-pQCT images of joints to separate bone from other tissues. The script outputs the proximal, distal, and full joint mask for each image.\n", + "\n", + " **Table of contents** \n", + " [Step 1: Imports](#imports) \n", + " [Step 2: Automatic Contour](#contour) \n", + " [Step 3: Display Results](#results) \n", + " [Step 4: Compare IPL to ORMIR_XCT](#compare)\n" + ] + }, + { + "cell_type": "markdown", + "id": "7ae0270b", + "metadata": {}, + "source": [ + "---" ] }, { @@ -17,7 +45,8 @@ "id": "bcdc9bff", "metadata": {}, "source": [ - "# Step 1: Imports\n", + "\n", + "## *Step 1: Imports:*\n", "\n", "Import modules/packages and set the input image path. " ] @@ -36,7 +65,10 @@ "from matplotlib import pyplot as plt\n", "\n", "from ormir_xct.autocontour.autocontour import autocontour\n", - "from ormir_xct.util.segmentation_evaluation import calculate_dice_and_jaccard, hausdorff_sitk" + "from ormir_xct.util.segmentation_evaluation import (\n", + " calculate_dice_and_jaccard,\n", + " hausdorff_sitk,\n", + ")" ] }, { @@ -54,12 +86,21 @@ "ipl_mask = sitk.ReadImage(joint_seg_ipl_path, sitk.sitkUInt8)" ] }, + { + "cell_type": "markdown", + "id": "21195166", + "metadata": {}, + "source": [ + "---" + ] + }, { "cell_type": "markdown", "id": "fbdc7d4f", "metadata": {}, "source": [ - "# Step 2: Run the ORMIR_XCT Automatic Contour\n", + "\n", + "## *Step 2: Run the ORMIR_XCT Automatic Contour:*\n", "\n", "Run the ORMIR_XCT automatic periosteal contour script on the input grayscale joint image. This script will return the distal, proximal, and full joint mask. The full joint mask will be used for comparison with the IPL periosteal contour workflow results.\n", "\n", @@ -80,10 +121,20 @@ "outputs": [], "source": [ "mu_water = 0.24090\n", - "rescale_slope = 1603.51904 \n", + "rescale_slope = 1603.51904\n", "rescale_intercept = -391.209015\n", "\n", - "dst_mask, prx_mask, ormir_mask = autocontour(gray_img, mu_water, rescale_slope, rescale_intercept)" + "dst_mask, prx_mask, ormir_mask = autocontour(\n", + " gray_img, mu_water, rescale_slope, rescale_intercept\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "5fba53c0", + "metadata": {}, + "source": [ + "---" ] }, { @@ -91,7 +142,8 @@ "id": "b81b793e", "metadata": {}, "source": [ - "# Step 3: Display Results\n", + "\n", + "## *Step 3: Display Results:*\n", "\n", "Now display the ORMIR_XCT and IPL segmentations together." ] @@ -117,13 +169,13 @@ "outputs": [], "source": [ "# Get the slices we want to view\n", - "view_slice1 = (slice(None), 45 - ormir_np.shape[1]//2, slice(None))\n", - "view_slice2 = (slice(None), ormir_np.shape[1]//2, slice(None))\n", - "view_slice3 = (slice(None), 45 + ormir_np.shape[1]//2, slice(None))\n", + "view_slice1 = (slice(None), 45 - ormir_np.shape[1] // 2, slice(None))\n", + "view_slice2 = (slice(None), ormir_np.shape[1] // 2, slice(None))\n", + "view_slice3 = (slice(None), 45 + ormir_np.shape[1] // 2, slice(None))\n", "\n", - "slice1 = abs(int(45 - ormir_np.shape[1]//2))\n", - "slice2 = int(ormir_np.shape[1]//2)\n", - "slice3 = int(45 + ormir_np.shape[1]//2)" + "slice1 = abs(int(45 - ormir_np.shape[1] // 2))\n", + "slice2 = int(ormir_np.shape[1] // 2)\n", + "slice3 = int(45 + ormir_np.shape[1] // 2)" ] }, { @@ -147,56 +199,69 @@ ], "source": [ "# Plot the segmentations overlaid onto the grayscale image for 3 slices\n", - "fig, axs = plt.subplots(3, 3, figsize=(50,30))\n", + "fig, axs = plt.subplots(3, 3, figsize=(50, 30))\n", "\n", "# Slice 1\n", - "axs[0][0].set_title('ORMIR_XCT Automatic Contour\\nSlice: {}'.format(slice1), fontsize = 40)\n", - "axs[0][0].imshow(gray_np[view_slice1], cmap='gray')\n", - "axs[0][0].imshow(ormir_np[view_slice1], cmap='rainbow', alpha=0.3)\n", + "axs[0][0].set_title(\n", + " \"ORMIR_XCT Automatic Contour\\nSlice: {}\".format(slice1), fontsize=40\n", + ")\n", + "axs[0][0].imshow(gray_np[view_slice1], cmap=\"gray\")\n", + "axs[0][0].imshow(ormir_np[view_slice1], cmap=\"rainbow\", alpha=0.3)\n", "\n", - "axs[0][1].set_title('IPL Automatic Contour\\nSlice: {}'.format(slice1), fontsize = 40)\n", - "axs[0][1].imshow(gray_np[view_slice1], cmap='gray')\n", - "axs[0][1].imshow(ipl_np[view_slice1], cmap='rainbow', alpha=0.3)\n", + "axs[0][1].set_title(\"IPL Automatic Contour\\nSlice: {}\".format(slice1), fontsize=40)\n", + "axs[0][1].imshow(gray_np[view_slice1], cmap=\"gray\")\n", + "axs[0][1].imshow(ipl_np[view_slice1], cmap=\"rainbow\", alpha=0.3)\n", "\n", - "axs[0][2].set_title('ORMIR_XCT and IPL Autocontour Overlaid\\nSlice: {}'.format(slice1), fontsize = 40)\n", - "axs[0][2].imshow(ormir_np[view_slice1], cmap='gray')\n", - "axs[0][2].imshow(ipl_np[view_slice1], cmap='rainbow', vmin=0, vmax=1, alpha=0.3)\n", + "axs[0][2].set_title(\n", + " \"ORMIR_XCT and IPL Autocontour Overlaid\\nSlice: {}\".format(slice1), fontsize=40\n", + ")\n", + "axs[0][2].imshow(ormir_np[view_slice1], cmap=\"gray\")\n", + "axs[0][2].imshow(ipl_np[view_slice1], cmap=\"rainbow\", vmin=0, vmax=1, alpha=0.3)\n", "\n", "# Slice 2\n", - "axs[1][0].set_title('Slice: {}'.format(slice2), fontsize = 40)\n", - "axs[1][0].imshow(gray_np[view_slice2], cmap='gray')\n", - "axs[1][0].imshow(ormir_np[view_slice2], cmap='rainbow', alpha=0.3)\n", + "axs[1][0].set_title(\"Slice: {}\".format(slice2), fontsize=40)\n", + "axs[1][0].imshow(gray_np[view_slice2], cmap=\"gray\")\n", + "axs[1][0].imshow(ormir_np[view_slice2], cmap=\"rainbow\", alpha=0.3)\n", "\n", - "axs[1][1].set_title('Slice: {}'.format(slice2), fontsize = 40)\n", - "axs[1][1].imshow(gray_np[view_slice2], cmap='gray')\n", - "axs[1][1].imshow(ipl_np[view_slice2], cmap='rainbow', alpha=0.3)\n", + "axs[1][1].set_title(\"Slice: {}\".format(slice2), fontsize=40)\n", + "axs[1][1].imshow(gray_np[view_slice2], cmap=\"gray\")\n", + "axs[1][1].imshow(ipl_np[view_slice2], cmap=\"rainbow\", alpha=0.3)\n", "\n", - "axs[1][2].set_title('Slice: {}'.format(slice2), fontsize = 40)\n", - "axs[1][2].imshow(ormir_np[view_slice2], cmap='gray')\n", - "axs[1][2].imshow(ipl_np[view_slice2], cmap='rainbow', vmin=0, vmax=1, alpha=0.3)\n", + "axs[1][2].set_title(\"Slice: {}\".format(slice2), fontsize=40)\n", + "axs[1][2].imshow(ormir_np[view_slice2], cmap=\"gray\")\n", + "axs[1][2].imshow(ipl_np[view_slice2], cmap=\"rainbow\", vmin=0, vmax=1, alpha=0.3)\n", "\n", "# Slice 3\n", - "axs[2][0].set_title('Slice: {}'.format(slice3), fontsize = 40)\n", - "axs[2][0].imshow(gray_np[view_slice3], cmap='gray')\n", - "axs[2][0].imshow(ormir_np[view_slice3], cmap='rainbow', alpha=0.3)\n", + "axs[2][0].set_title(\"Slice: {}\".format(slice3), fontsize=40)\n", + "axs[2][0].imshow(gray_np[view_slice3], cmap=\"gray\")\n", + "axs[2][0].imshow(ormir_np[view_slice3], cmap=\"rainbow\", alpha=0.3)\n", "\n", - "axs[2][1].set_title('Slice: {}'.format(slice3), fontsize = 40)\n", - "axs[2][1].imshow(gray_np[view_slice3], cmap='gray')\n", - "axs[2][1].imshow(ipl_np[view_slice3], cmap='rainbow', alpha=0.3)\n", + "axs[2][1].set_title(\"Slice: {}\".format(slice3), fontsize=40)\n", + "axs[2][1].imshow(gray_np[view_slice3], cmap=\"gray\")\n", + "axs[2][1].imshow(ipl_np[view_slice3], cmap=\"rainbow\", alpha=0.3)\n", "\n", - "axs[2][2].set_title('Slice: {}'.format(slice3), fontsize = 40)\n", - "axs[2][2].imshow(ormir_np[view_slice3], cmap='gray')\n", - "axs[2][2].imshow(ipl_np[view_slice3], cmap='rainbow', vmin=0, vmax=1, alpha=0.3)\n", + "axs[2][2].set_title(\"Slice: {}\".format(slice3), fontsize=40)\n", + "axs[2][2].imshow(ormir_np[view_slice3], cmap=\"gray\")\n", + "axs[2][2].imshow(ipl_np[view_slice3], cmap=\"rainbow\", vmin=0, vmax=1, alpha=0.3)\n", "\n", "plt.show()" ] }, + { + "cell_type": "markdown", + "id": "df74fe73", + "metadata": {}, + "source": [ + "---" + ] + }, { "cell_type": "markdown", "id": "747bf6db", "metadata": {}, "source": [ - "# Step 4: Compare Between IPL and ORMIR_XCT\n", + "\n", + "## *Step 4: Compare Between IPL and ORMIR_XCT:*\n", "\n", "Now compare segmentations using metrics.\n", "\n", @@ -210,7 +275,9 @@ "metadata": {}, "outputs": [], "source": [ - "resampled_ormir = sitk.Resample(ormir_mask, ipl_mask, interpolator=sitk.sitkNearestNeighbor)\n", + "resampled_ormir = sitk.Resample(\n", + " ormir_mask, ipl_mask, interpolator=sitk.sitkNearestNeighbor\n", + ")\n", "\n", "ormir_np = sitk.GetArrayFromImage(resampled_ormir)" ] @@ -236,11 +303,30 @@ "dice, jaccard = calculate_dice_and_jaccard(ipl_np, ormir_np)\n", "hausdorff = hausdorff_sitk(ipl_mask, resampled_ormir)\n", "\n", - "print('DICE: ', dice)\n", - "print('Jaccard: ', jaccard)\n", - "print('Mean Hausdorff Distance: ', hausdorff[0])\n", - "print('Maximum Hausdorff Distance: ', hausdorff[1])" + "print(\"DICE: \", dice)\n", + "print(\"Jaccard: \", jaccard)\n", + "print(\"Mean Hausdorff Distance: \", hausdorff[0])\n", + "print(\"Maximum Hausdorff Distance: \", hausdorff[1])" ] + }, + { + "cell_type": "markdown", + "id": "e2eb10cb", + "metadata": {}, + "source": [ + "---\n", + "\n", + "\n", + "Notebook created using the [template](https://github.com/ORMIRcommunity/templates/blob/main/ORMIR_nb_template.ipynb) of the [ORMIR community](https://ormircommunity.github.io/) (version 1.0, 2023)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fa654d72", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/examples/Bone_Mineral_Density_Example.ipynb b/examples/Bone_Mineral_Density_Example.ipynb index 031ccaf..c34b877 100644 --- a/examples/Bone_Mineral_Density_Example.ipynb +++ b/examples/Bone_Mineral_Density_Example.ipynb @@ -7,7 +7,35 @@ "source": [ "# Bone Mineral Density Example\n", "\n", - "This Jupyter Notebook provides an example of computing volumetric bone mineral density (vBMD) in an HR-pQCT joint image. vBMD is computed in the distal and proximal segments. The ORMIR_XCT autocontour script is used to generate the bone masks." + "- **By:** [Michael T. Kuczynski](https://www.linkedin.com/in/mkuczyns/), 2024 \n", + "- **License:** CC-BY \n", + "- **How to cite:** Cite the ORMIR_XCT publication: *Kuczynski, M.T., et al. \"ORMIR_XCT: A Python package for high resolution peripheral quantitative computed tomography image processing.\" arXiv preprint arXiv:2309.04602 (2023).*" + ] + }, + { + "cell_type": "markdown", + "id": "133c761a", + "metadata": {}, + "source": [ + "---\n", + "# Aims\n", + "\n", + "- This Jupyter Notebook provides an example of computing volumetric bone mineral density (vBMD) in an HR-pQCT joint image. \n", + "- vBMD is computed in the distal and proximal segments. The ORMIR_XCT autocontour script is used to generate the bone masks.\n", + "\n", + " **Table of contents** \n", + " [Step 1: Imports](#imports) \n", + " [Step 2: Automatic Contour](#contour) \n", + " [Step 3: Compute vBMD](#bmd) \n", + " [Step 4: Compare IPL to ORMIR_XCT](#compare)" + ] + }, + { + "cell_type": "markdown", + "id": "2415563c", + "metadata": {}, + "source": [ + "---" ] }, { @@ -15,7 +43,8 @@ "id": "bcdc9bff", "metadata": {}, "source": [ - "# Step 1: Imports\n", + "\n", + "## *Step 1: Imports:*\n", "\n", "Import modules/packages and set the input image path. " ] @@ -50,12 +79,21 @@ "gray_img = sitk.ReadImage(joint_seg_path, sitk.sitkFloat32)" ] }, + { + "cell_type": "markdown", + "id": "0b10f2a5", + "metadata": {}, + "source": [ + "---" + ] + }, { "cell_type": "markdown", "id": "fbdc7d4f", "metadata": {}, "source": [ - "# Step 2: Run the ORMIR_XCT Automatic Contour\n", + "\n", + "## *Step 2: Run the ORMIR_XCT Automatic Contour:*\n", "\n", "Run the ORMIR_XCT automatic periosteal contour script on the input grayscale joint image. This script will return the distal, proximal, and full joint mask.\n", "\n", @@ -78,10 +116,20 @@ "source": [ "mu_water = 0.24090\n", "mu_scaling = 8192\n", - "rescale_slope = 1603.51904 \n", + "rescale_slope = 1603.51904\n", "rescale_intercept = -391.209015\n", "\n", - "dst_mask, prx_mask, ormir_mask = autocontour(gray_img, mu_water, rescale_slope, rescale_intercept)" + "dst_mask, prx_mask, ormir_mask = autocontour(\n", + " gray_img, mu_water, rescale_slope, rescale_intercept\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "64cda5b5", + "metadata": {}, + "source": [ + "---" ] }, { @@ -89,7 +137,8 @@ "id": "b81b793e", "metadata": {}, "source": [ - "# Step 3: Compute vBMD\n", + "\n", + "## *Step 3: Compute vBMD:*\n", "\n", "Now computed vBMD. We will calculate the vBMD for each bone segment." ] @@ -114,16 +163,40 @@ "\n", "# Distal vBMD\n", "masked_gray_img = sitk.Mask(gray_img, dst_mask)\n", - "dst_bmd = bmd_masked(masked_gray_img, dst_mask, 'hu', mu_scaling, mu_water, rescale_slope, rescale_intercept)\n", + "dst_bmd = bmd_masked(\n", + " masked_gray_img,\n", + " dst_mask,\n", + " \"hu\",\n", + " mu_scaling,\n", + " mu_water,\n", + " rescale_slope,\n", + " rescale_intercept,\n", + ")\n", "\n", - "print('Distal vBMD: {:.2f} +/- {:.2f} mg HA/cm^3'.format(dst_bmd[0], dst_bmd[1]))\n", + "print(\"Distal vBMD: {:.2f} +/- {:.2f} mg HA/cm^3\".format(dst_bmd[0], dst_bmd[1]))\n", "\n", "\n", "# Proximal vBMD\n", "prx_gray_img = sitk.Mask(gray_img, prx_mask)\n", - "prx_bmd = bmd_masked(masked_gray_img, prx_mask, 'hu', mu_scaling, mu_water, rescale_slope, rescale_intercept)\n", + "prx_bmd = bmd_masked(\n", + " masked_gray_img,\n", + " prx_mask,\n", + " \"hu\",\n", + " mu_scaling,\n", + " mu_water,\n", + " rescale_slope,\n", + " rescale_intercept,\n", + ")\n", "\n", - "print('Proximal vBMD: {:.2f} +/- {:.2f} mg HA/cm^3'.format(prx_bmd[0], prx_bmd[1]))" + "print(\"Proximal vBMD: {:.2f} +/- {:.2f} mg HA/cm^3\".format(prx_bmd[0], prx_bmd[1]))" + ] + }, + { + "cell_type": "markdown", + "id": "6d8e7543", + "metadata": {}, + "source": [ + "---" ] }, { @@ -131,7 +204,8 @@ "id": "747bf6db", "metadata": {}, "source": [ - "# Step 4: Compare Between IPL and ORMIR_XCT\n", + "\n", + "## *Step 4: Compare Between IPL and ORMIR_XCT:*\n", "\n", "Now compare vBMD results between IPL and ORMIR_XCT.\n", "\n", @@ -139,6 +213,25 @@ "- **Distal vBMD:** 203.84 +/- 212.96 mg HA/cm^3\n", "- **Proximal vBMD:** 224.20 +/- 240.37 mg HA/cm^3" ] + }, + { + "cell_type": "markdown", + "id": "4755148b", + "metadata": {}, + "source": [ + "---\n", + "\n", + "\n", + "Notebook created using the [template](https://github.com/ORMIRcommunity/templates/blob/main/ORMIR_nb_template.ipynb) of the [ORMIR community](https://ormircommunity.github.io/) (version 1.0, 2023)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "41b94b23", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/examples/Trabecular_Segmentation_Example.ipynb b/examples/Trabecular_Segmentation_Example.ipynb index 04ca702..4921a5d 100644 --- a/examples/Trabecular_Segmentation_Example.ipynb +++ b/examples/Trabecular_Segmentation_Example.ipynb @@ -7,7 +7,34 @@ "source": [ "# Trabecular Segmentation Example\n", "\n", - "This Jupyter Notebook provides an example of running the trabecular segmentation script in the ORMIR_XCT package. The results of this segmentation are compared to results from the standard IPL workflow." + "- **By:** [Michael T. Kuczynski](https://www.linkedin.com/in/mkuczyns/), 2024 \n", + "- **License:** CC-BY \n", + "- **How to cite:** Cite the ORMIR_XCT publication: *Kuczynski, M.T., et al. \"ORMIR_XCT: A Python package for high resolution peripheral quantitative computed tomography image processing.\" arXiv preprint arXiv:2309.04602 (2023).*" + ] + }, + { + "cell_type": "markdown", + "id": "d62af34c", + "metadata": {}, + "source": [ + "---\n", + "# Aims\n", + "\n", + "- This Jupyter Notebook provides an example of running the trabecular segmentation script in the ORMIR_XCT package. The results of this segmentation are compared to results from the standard IPL workflow.\n", + "\n", + " **Table of contents** \n", + " [Step 1: Imports](#Imports) \n", + " [Step 2: Trabecular Segmentation](#Segmentation) \n", + " [Step 3: Display Results](#Results) \n", + " [Step 4: Compare IPL to ORMIR_XCT](#Compare)" + ] + }, + { + "cell_type": "markdown", + "id": "27fbbd97", + "metadata": {}, + "source": [ + "---" ] }, { @@ -15,7 +42,8 @@ "id": "bcdc9bff", "metadata": {}, "source": [ - "# Step 1: Imports\n", + "\n", + "## *Step 1: Imports*\n", "\n", "Import modules/packages and set the input image path. " ] @@ -34,7 +62,10 @@ "from matplotlib import pyplot as plt\n", "\n", "from ormir_xct.segmentation.ipl_seg import ipl_seg, threshold_dict\n", - "from ormir_xct.util.segmentation_evaluation import calculate_dice_and_jaccard, hausdorff_sitk" + "from ormir_xct.util.segmentation_evaluation import (\n", + " calculate_dice_and_jaccard,\n", + " hausdorff_sitk,\n", + ")" ] }, { @@ -52,12 +83,21 @@ "ipl_mask = sitk.ReadImage(joint_seg_ipl_path, sitk.sitkUInt8)" ] }, + { + "cell_type": "markdown", + "id": "d3a6851d", + "metadata": {}, + "source": [ + "---" + ] + }, { "cell_type": "markdown", "id": "680547ef", "metadata": {}, "source": [ - "# Step 2: Run the ORMIR_XCT Trabecular Segmentation\n", + "\n", + "## *Step 2: Trabecular Segmentation*\n", "\n", "Run the ORMIR_XCT trabecular segmentation script on the input grayscale joint image. \n", "\n", @@ -73,12 +113,26 @@ "source": [ "mu_water = 0.24090\n", "mu_scaling = 8192\n", - "resale_slope = 1603.51904 \n", + "resale_slope = 1603.51904\n", "rescale_intercept = -391.209015\n", - "lower_thresh = threshold_dict['HU_Lower']\n", - "upper_thresh = threshold_dict['HU_Upper']\n", + "lower_thresh = threshold_dict[\"HU_Lower\"]\n", + "upper_thresh = threshold_dict[\"HU_Upper\"]\n", "\n", - "trab_seg_ormir = ipl_seg(gray_img, lower_thresh, upper_thresh, voxel_size=0.0606964, sigma=0.5,)" + "trab_seg_ormir = ipl_seg(\n", + " gray_img,\n", + " lower_thresh,\n", + " upper_thresh,\n", + " voxel_size=0.0606964,\n", + " sigma=0.5,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "f89228de", + "metadata": {}, + "source": [ + "---" ] }, { @@ -86,7 +140,8 @@ "id": "b81b793e", "metadata": {}, "source": [ - "# Step 3: Display Results\n", + "\n", + "## *Step 3: Display Results*\n", "\n", "Now display the ORMIR_XCT and IPL segmentations together." ] @@ -98,7 +153,9 @@ "metadata": {}, "outputs": [], "source": [ - "ipl_mask = sitk.Resample(ipl_mask, trab_seg_ormir, interpolator=sitk.sitkNearestNeighbor)\n", + "ipl_mask = sitk.Resample(\n", + " ipl_mask, trab_seg_ormir, interpolator=sitk.sitkNearestNeighbor\n", + ")\n", "\n", "# Convert the images to numpy arrays for plotting with Matplotlib\n", "gray_np = sitk.GetArrayFromImage(gray_img)\n", @@ -114,13 +171,13 @@ "outputs": [], "source": [ "# Get the slices we want to view\n", - "view_slice1 = (slice(None), 45 - ormir_np.shape[1]//2, slice(None))\n", - "view_slice2 = (slice(None), ormir_np.shape[1]//2, slice(None))\n", - "view_slice3 = (slice(None), 45 + ormir_np.shape[1]//2, slice(None))\n", + "view_slice1 = (slice(None), 45 - ormir_np.shape[1] // 2, slice(None))\n", + "view_slice2 = (slice(None), ormir_np.shape[1] // 2, slice(None))\n", + "view_slice3 = (slice(None), 45 + ormir_np.shape[1] // 2, slice(None))\n", "\n", - "slice1 = abs(int(45 - ormir_np.shape[1]//2))\n", - "slice2 = int(ormir_np.shape[1]//2)\n", - "slice3 = int(45 + ormir_np.shape[1]//2)" + "slice1 = abs(int(45 - ormir_np.shape[1] // 2))\n", + "slice2 = int(ormir_np.shape[1] // 2)\n", + "slice3 = int(45 + ormir_np.shape[1] // 2)" ] }, { @@ -144,56 +201,69 @@ ], "source": [ "# Plot the segmentations overlaid onto the grayscale image for 3 slices\n", - "fig, axs = plt.subplots(3, 3, figsize=(50,30))\n", + "fig, axs = plt.subplots(3, 3, figsize=(50, 30))\n", "\n", "# Slice 1\n", - "axs[0][0].set_title('ORMIR_XCT Trabecular Segmentation\\nSlice: {}'.format(slice1), fontsize = 40)\n", - "axs[0][0].imshow(gray_np[view_slice1], cmap='gray')\n", - "axs[0][0].imshow(ormir_np[view_slice1], cmap='rainbow', alpha=0.3)\n", + "axs[0][0].set_title(\n", + " \"ORMIR_XCT Trabecular Segmentation\\nSlice: {}\".format(slice1), fontsize=40\n", + ")\n", + "axs[0][0].imshow(gray_np[view_slice1], cmap=\"gray\")\n", + "axs[0][0].imshow(ormir_np[view_slice1], cmap=\"rainbow\", alpha=0.3)\n", "\n", - "axs[0][1].set_title('IPL Trabecular Segmentation\\nSlice: {}'.format(slice1), fontsize = 40)\n", - "axs[0][1].imshow(gray_np[view_slice1], cmap='gray')\n", - "axs[0][1].imshow(ipl_np[view_slice1], cmap='rainbow', alpha=0.3)\n", + "axs[0][1].set_title(\n", + " \"IPL Trabecular Segmentation\\nSlice: {}\".format(slice1), fontsize=40\n", + ")\n", + "axs[0][1].imshow(gray_np[view_slice1], cmap=\"gray\")\n", + "axs[0][1].imshow(ipl_np[view_slice1], cmap=\"rainbow\", alpha=0.3)\n", "\n", - "axs[0][2].set_title('ORMIR_XCT and IPL Overlaid\\nSlice: {}'.format(slice1), fontsize = 40)\n", - "axs[0][2].imshow(ormir_np[view_slice1], cmap='gray')\n", - "axs[0][2].imshow(ipl_np[view_slice1], cmap='rainbow', vmin=0, vmax=1, alpha=0.3)\n", + "axs[0][2].set_title(\"ORMIR_XCT and IPL Overlaid\\nSlice: {}\".format(slice1), fontsize=40)\n", + "axs[0][2].imshow(ormir_np[view_slice1], cmap=\"gray\")\n", + "axs[0][2].imshow(ipl_np[view_slice1], cmap=\"rainbow\", vmin=0, vmax=1, alpha=0.3)\n", "\n", "# Slice 2\n", - "axs[1][0].set_title('Slice: {}'.format(slice2), fontsize = 40)\n", - "axs[1][0].imshow(gray_np[view_slice2], cmap='gray')\n", - "axs[1][0].imshow(ormir_np[view_slice2], cmap='rainbow', alpha=0.3)\n", + "axs[1][0].set_title(\"Slice: {}\".format(slice2), fontsize=40)\n", + "axs[1][0].imshow(gray_np[view_slice2], cmap=\"gray\")\n", + "axs[1][0].imshow(ormir_np[view_slice2], cmap=\"rainbow\", alpha=0.3)\n", "\n", - "axs[1][1].set_title('Slice: {}'.format(slice2), fontsize = 40)\n", - "axs[1][1].imshow(gray_np[view_slice2], cmap='gray')\n", - "axs[1][1].imshow(ipl_np[view_slice2], cmap='rainbow', alpha=0.3)\n", + "axs[1][1].set_title(\"Slice: {}\".format(slice2), fontsize=40)\n", + "axs[1][1].imshow(gray_np[view_slice2], cmap=\"gray\")\n", + "axs[1][1].imshow(ipl_np[view_slice2], cmap=\"rainbow\", alpha=0.3)\n", "\n", - "axs[1][2].set_title('Slice: {}'.format(slice2), fontsize = 40)\n", - "axs[1][2].imshow(ormir_np[view_slice2], cmap='gray')\n", - "axs[1][2].imshow(ipl_np[view_slice2], cmap='rainbow', vmin=0, vmax=1, alpha=0.3)\n", + "axs[1][2].set_title(\"Slice: {}\".format(slice2), fontsize=40)\n", + "axs[1][2].imshow(ormir_np[view_slice2], cmap=\"gray\")\n", + "axs[1][2].imshow(ipl_np[view_slice2], cmap=\"rainbow\", vmin=0, vmax=1, alpha=0.3)\n", "\n", "# Slice 3\n", - "axs[2][0].set_title('Slice: {}'.format(slice3), fontsize = 40)\n", - "axs[2][0].imshow(gray_np[view_slice3], cmap='gray')\n", - "axs[2][0].imshow(ormir_np[view_slice3], cmap='rainbow', alpha=0.3)\n", + "axs[2][0].set_title(\"Slice: {}\".format(slice3), fontsize=40)\n", + "axs[2][0].imshow(gray_np[view_slice3], cmap=\"gray\")\n", + "axs[2][0].imshow(ormir_np[view_slice3], cmap=\"rainbow\", alpha=0.3)\n", "\n", - "axs[2][1].set_title('Slice: {}'.format(slice3), fontsize = 40)\n", - "axs[2][1].imshow(gray_np[view_slice3], cmap='gray')\n", - "axs[2][1].imshow(ipl_np[view_slice3], cmap='rainbow', alpha=0.3)\n", + "axs[2][1].set_title(\"Slice: {}\".format(slice3), fontsize=40)\n", + "axs[2][1].imshow(gray_np[view_slice3], cmap=\"gray\")\n", + "axs[2][1].imshow(ipl_np[view_slice3], cmap=\"rainbow\", alpha=0.3)\n", "\n", - "axs[2][2].set_title('Slice: {}'.format(slice3), fontsize = 40)\n", - "axs[2][2].imshow(ormir_np[view_slice3], cmap='gray')\n", - "axs[2][2].imshow(ipl_np[view_slice3], cmap='rainbow', vmin=0, vmax=1, alpha=0.3)\n", + "axs[2][2].set_title(\"Slice: {}\".format(slice3), fontsize=40)\n", + "axs[2][2].imshow(ormir_np[view_slice3], cmap=\"gray\")\n", + "axs[2][2].imshow(ipl_np[view_slice3], cmap=\"rainbow\", vmin=0, vmax=1, alpha=0.3)\n", "\n", "plt.show()" ] }, + { + "cell_type": "markdown", + "id": "c98d7c03", + "metadata": {}, + "source": [ + "---" + ] + }, { "cell_type": "markdown", "id": "747bf6db", "metadata": {}, "source": [ - "# Step 4: Compare Between IPL and ORMIR_XCT\n", + "\n", + "## *Step 4: Compare IPL to ORMIR_XCT*\n", "\n", "Now compare segmentations using metrics.\n", "\n", @@ -221,16 +291,27 @@ "dice, jaccard = calculate_dice_and_jaccard(ipl_np, ormir_np)\n", "hausdorff = hausdorff_sitk(ipl_mask, trab_seg_ormir)\n", "\n", - "print('DICE: ', dice)\n", - "print('Jaccard: ', jaccard)\n", - "print('Mean Hausdorff Distance: ', hausdorff[0])\n", - "print('Maximum Hausdorff Distance: ', hausdorff[1])" + "print(\"DICE: \", dice)\n", + "print(\"Jaccard: \", jaccard)\n", + "print(\"Mean Hausdorff Distance: \", hausdorff[0])\n", + "print(\"Maximum Hausdorff Distance: \", hausdorff[1])" + ] + }, + { + "cell_type": "markdown", + "id": "31129519", + "metadata": {}, + "source": [ + "---\n", + "\n", + "\n", + "Notebook created using the [template](https://github.com/ORMIRcommunity/templates/blob/main/ORMIR_nb_template.ipynb) of the [ORMIR community](https://ormircommunity.github.io/) (version 1.0, 2023)" ] }, { "cell_type": "code", "execution_count": null, - "id": "61ff81aa", + "id": "7b576ffa", "metadata": {}, "outputs": [], "source": [] diff --git a/examples/finger_joint_comparison.ipynb b/examples/finger_joint_comparison.ipynb index 9047181..6891642 100644 --- a/examples/finger_joint_comparison.ipynb +++ b/examples/finger_joint_comparison.ipynb @@ -7,22 +7,47 @@ "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", + "- **By:** [Michael T. Kuczynski](https://www.linkedin.com/in/mkuczyns/), [Nathan Neeteson](https://www.linkedin.com/in/nathan-neeteson/), 2023 \n", + "- **License:** CC-BY \n", + "- **How to cite:** Cite the ORMIR_XCT publication: *Kuczynski, M.T., et al. \"ORMIR_XCT: A Python package for high resolution peripheral quantitative computed tomography image processing.\" arXiv preprint arXiv:2309.04602 (2023).*" + ] + }, + { + "cell_type": "markdown", + "id": "253e760c", + "metadata": {}, + "source": [ + "---\n", + "# Aims\n", "\n", - "The following distance transform image was generated in IPL:\n", + "- This notebook compares results from the IPL joint space width (JSW) analysis workflow to the results from the ORMIR_XCT package.\n", + "- 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", + "- The following distance transform image was generated in IPL:\n", "\n", "\n", - "The following JSW results were obtained from the above distance transform image:\n", - "\n", + "- The following JSW results were obtained from the above distance transform image:\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 |" + "| JS Asymmetry (JS Max / JS Min) | 1.9545 |\n", + "\n", + " **Table of contents** \n", + " [Step 1: Imports](#imports) \n", + " [Step 2: JSW Workflow](#jsw) \n", + " [Step 3: JSW Parameters](#params) \n", + " [Step 4: Display Results](#results)\n" + ] + }, + { + "cell_type": "markdown", + "id": "00519f9e", + "metadata": {}, + "source": [ + "---" ] }, { @@ -30,7 +55,8 @@ "id": "ebee8652", "metadata": {}, "source": [ - "# Step 1: Imports\n", + "\n", + "## *Step 1: Imports:*\n", "\n", "Import modules/packages and set the input image path. " ] @@ -47,7 +73,12 @@ "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" + "from ormir_xct.joint_space_analysis.jsw_morphometry import (\n", + " jsw_pad,\n", + " jsw_dilate,\n", + " jsw_erode,\n", + " jsw_parameters,\n", + ")" ] }, { @@ -61,12 +92,21 @@ "output_path = \"images\"" ] }, + { + "cell_type": "markdown", + "id": "3c05c22d", + "metadata": {}, + "source": [ + "---" + ] + }, { "cell_type": "markdown", "id": "79816632", "metadata": {}, "source": [ - "# Step 2: Run the JSW Workflow\n", + "\n", + "## *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", @@ -109,7 +149,17 @@ "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\"))" + "sitk.WriteImage(\n", + " dilated_js_mask, os.path.join(output_path, str(basename) + \"_DILATED_JS_MASK.nii\")\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "1b6a9cc2", + "metadata": {}, + "source": [ + "---" ] }, { @@ -117,7 +167,8 @@ "id": "44fc2d09", "metadata": {}, "source": [ - "# Step 3: Compute the JSW Parameters\n", + "\n", + "## *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." ] @@ -132,36 +183,73 @@ "# 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", + "dt_img1, jsw_params1 = jsw_parameters(\n", + " pad_image,\n", + " dilated_js_mask,\n", + " basename,\n", + " output_path,\n", + " js_mask,\n", + " pad_image.GetSpacing()[0],\n", + " oversamp=False,\n", + " skel=False,\n", + ")\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", + "dt_img2, jsw_params2 = jsw_parameters(\n", + " pad_image,\n", + " dilated_js_mask,\n", + " basename,\n", + " output_path,\n", + " js_mask,\n", + " pad_image.GetSpacing()[0],\n", + " oversamp=True,\n", + " skel=False,\n", + ")\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", + "dt_img3, jsw_params3 = jsw_parameters(\n", + " pad_image,\n", + " dilated_js_mask,\n", + " basename,\n", + " output_path,\n", + " js_mask,\n", + " pad_image.GetSpacing()[0],\n", + " oversamp=False,\n", + " skel=True,\n", + ")\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", + "dt_img4, jsw_params4 = jsw_parameters(\n", + " pad_image,\n", + " dilated_js_mask,\n", + " basename,\n", + " output_path,\n", + " js_mask,\n", + " pad_image.GetSpacing()[0],\n", + " oversamp=True,\n", + " skel=True,\n", + ")\n", "sitk.WriteImage(dt_img1, os.path.join(output_path, str(basename) + \"_DT4.nii\"))" ] }, + { + "cell_type": "markdown", + "id": "16f78384", + "metadata": {}, + "source": [ + "---" + ] + }, { "cell_type": "markdown", "id": "68b9841c", "metadata": {}, "source": [ - "# Step 4: Display Results\n", + "\n", + "## *Step 4: Display Results:*\n", "\n", "Use Pandas to print a table of the results we generated." ] @@ -283,26 +371,46 @@ "\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", + "data = {\n", + " \"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(\n", + " data,\n", + " index=[\n", + " \"JS Volume (mm^3)\",\n", + " \"Mean JS Width (mm)\",\n", + " \"JS Max (mm)\",\n", + " \"JS Min (mm)\",\n", + " \"JS Asymmetry (JS Max / JS Min)\",\n", + " ],\n", + ")\n", "\n", "# Set table default options\n", - "pd.set_option('display.max_colwidth', None)\n", - "pd.set_option('colheader_justify', 'center')\n", + "pd.set_option(\"display.max_colwidth\", None)\n", + "pd.set_option(\"colheader_justify\", \"center\")\n", "\n", "df.style.format(precision=2)" ] }, + { + "cell_type": "markdown", + "id": "e9818261", + "metadata": {}, + "source": [ + "---\n", + "\n", + "\n", + "Notebook created using the [template](https://github.com/ORMIRcommunity/templates/blob/main/ORMIR_nb_template.ipynb) of the [ORMIR community](https://ormircommunity.github.io/) (version 1.0, 2023)" + ] + }, { "cell_type": "code", "execution_count": null, - "id": "21bce143", + "id": "f7cc9181", "metadata": {}, "outputs": [], "source": [] @@ -324,7 +432,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.8" + "version": "3.8.5" } }, "nbformat": 4, diff --git a/examples/thickness_of_shapes.ipynb b/examples/thickness_of_shapes.ipynb index 911456a..925fa15 100644 --- a/examples/thickness_of_shapes.ipynb +++ b/examples/thickness_of_shapes.ipynb @@ -7,7 +7,42 @@ "source": [ "# Thickness of Shapes\n", "\n", - "Investigate the local thickness of structures in Python using the ORMIR_XCT package." + "- **By:** [Michael T. Kuczynski](https://www.linkedin.com/in/mkuczyns/), [Nathan Neeteson](https://www.linkedin.com/in/nathan-neeteson/), 2023 \n", + "- **License:** CC-BY \n", + "- **How to cite:** Cite the ORMIR_XCT publication: *Kuczynski, M.T., et al. \"ORMIR_XCT: A Python package for high resolution peripheral quantitative computed tomography image processing.\" arXiv preprint arXiv:2309.04602 (2023).*" + ] + }, + { + "cell_type": "markdown", + "id": "e8797811", + "metadata": {}, + "source": [ + "---\n", + "# Aims\n", + "\n", + "- Investigate the local thickness of structures in Python using the ORMIR_XCT package.\n", + "\n", + " **Table of contents** \n", + " [Step 1: Imports](#imports) \n", + " [Step 2: Generating Shapes](#shapes) \n", + " [Step 3: Check Estimated Thickness](#check) " + ] + }, + { + "cell_type": "markdown", + "id": "174e2413", + "metadata": {}, + "source": [ + "---" + ] + }, + { + "cell_type": "markdown", + "id": "e802e6a9", + "metadata": {}, + "source": [ + "\n", + "## *Step 1: Imports:*" ] }, { @@ -24,12 +59,21 @@ "from ormir_xct.util.hildebrand_thickness import calc_structure_thickness_statistics" ] }, + { + "cell_type": "markdown", + "id": "d2b9bda9", + "metadata": {}, + "source": [ + "---" + ] + }, { "cell_type": "markdown", "id": "86bddd4f", "metadata": {}, "source": [ - "Synthetic shape mask functions:" + "\n", + "## *Step 2: Synthetic shape mask functions:*" ] }, { @@ -41,22 +85,29 @@ "source": [ "def create_shape(shape, voxel_widths, thickness, shape_type=\"sphere\"):\n", " center = (\n", - " voxel_widths[0] * (shape[0] // 2), \n", - " voxel_widths[1] * (shape[1] // 2), \n", - " voxel_widths[2] * (shape[2] // 2)\n", + " voxel_widths[0] * (shape[0] // 2),\n", + " voxel_widths[1] * (shape[1] // 2),\n", + " voxel_widths[2] * (shape[2] // 2),\n", + " )\n", + " x, y, z = np.meshgrid(\n", + " *[voxel_widths[i] * np.arange(0, shape[i]) for i in range(3)], indexing=\"ij\"\n", " )\n", - " x, y, z = np.meshgrid(*[voxel_widths[i] * np.arange(0, shape[i]) for i in range(3)], indexing=\"ij\")\n", " if shape_type == \"sphere\":\n", - " mask = (((x-center[0])**2 + (y-center[1])**2 + (z-center[2])**2) < (thickness/2)**2).astype(int)\n", + " mask = (\n", + " ((x - center[0]) ** 2 + (y - center[1]) ** 2 + (z - center[2]) ** 2)\n", + " < (thickness / 2) ** 2\n", + " ).astype(int)\n", " elif shape_type == \"cylinder\":\n", - " mask = (((x-center[0])**2 + (y-center[1])**2) < (thickness/2)**2).astype(int)\n", + " mask = (\n", + " ((x - center[0]) ** 2 + (y - center[1]) ** 2) < (thickness / 2) ** 2\n", + " ).astype(int)\n", " elif shape_type == \"plate\":\n", - " mask = (np.abs(x-center[0]) < thickness/2).astype(int)\n", + " mask = (np.abs(x - center[0]) < thickness / 2).astype(int)\n", " else:\n", - " raise ValueError(f\"`shape_type` can be `sphere`, `cylinder`, `plate`; got {shape_type}\")\n", - " return mask\n", - "\n", - " " + " raise ValueError(\n", + " f\"`shape_type` can be `sphere`, `cylinder`, `plate`; got {shape_type}\"\n", + " )\n", + " return mask" ] }, { @@ -74,8 +125,8 @@ "metadata": {}, "outputs": [], "source": [ - "shape = tuple([50]*3)\n", - "voxel_widths = tuple([0.0607]*3)\n", + "shape = tuple([50] * 3)\n", + "voxel_widths = tuple([0.0607] * 3)\n", "radius = 1\n", "\n", "sphere = create_shape(shape, voxel_widths, radius, shape_type=\"sphere\")\n", @@ -83,12 +134,21 @@ "plate = create_shape(shape, voxel_widths, radius, shape_type=\"plate\")" ] }, + { + "cell_type": "markdown", + "id": "40f1ba7a", + "metadata": {}, + "source": [ + "---" + ] + }, { "cell_type": "markdown", "id": "9383b9fa", "metadata": {}, "source": [ - "Check their estimated mean thicknesses:" + "\n", + "## *Step 3: Check their estimated mean thicknesses:*" ] }, { @@ -99,7 +159,9 @@ "outputs": [], "source": [ "sphere_thickness_stats = calc_structure_thickness_statistics(sphere, voxel_widths, 0)\n", - "cylinder_thickness_stats = calc_structure_thickness_statistics(cylinder, voxel_widths, 0)\n", + "cylinder_thickness_stats = calc_structure_thickness_statistics(\n", + " cylinder, voxel_widths, 0\n", + ")\n", "plate_thickness_stats = calc_structure_thickness_statistics(plate, voxel_widths, 0)" ] }, @@ -120,9 +182,15 @@ } ], "source": [ - "print(f\"Sphere thickness is {sphere_thickness_stats[0]:0.3f} +/- {sphere_thickness_stats[1]:0.3f}\")\n", - "print(f\"Cylinder thickness is {cylinder_thickness_stats[0]:0.3f} +/- {cylinder_thickness_stats[1]:0.3f}\")\n", - "print(f\"Plate thickness is {plate_thickness_stats[0]:0.3f} +/- {plate_thickness_stats[1]:0.3f}\")" + "print(\n", + " f\"Sphere thickness is {sphere_thickness_stats[0]:0.3f} +/- {sphere_thickness_stats[1]:0.3f}\"\n", + ")\n", + "print(\n", + " f\"Cylinder thickness is {cylinder_thickness_stats[0]:0.3f} +/- {cylinder_thickness_stats[1]:0.3f}\"\n", + ")\n", + "print(\n", + " f\"Plate thickness is {plate_thickness_stats[0]:0.3f} +/- {plate_thickness_stats[1]:0.3f}\"\n", + ")" ] }, { @@ -175,15 +243,21 @@ "\n", "for thickness in true_thicknesses:\n", " print(f\"thickness: {thickness}\")\n", - " \n", + "\n", " sphere = create_shape(shape, voxel_widths, thickness, shape_type=\"sphere\")\n", - " est_sph_thicknesses.append(calc_structure_thickness_statistics(sphere, voxel_widths, 0)[0])\n", - " \n", + " est_sph_thicknesses.append(\n", + " calc_structure_thickness_statistics(sphere, voxel_widths, 0)[0]\n", + " )\n", + "\n", " cylinder = create_shape(shape, voxel_widths, thickness, shape_type=\"cylinder\")\n", - " est_cyl_thicknesses.append(calc_structure_thickness_statistics(cylinder, voxel_widths, 0)[0])\n", - " \n", + " est_cyl_thicknesses.append(\n", + " calc_structure_thickness_statistics(cylinder, voxel_widths, 0)[0]\n", + " )\n", + "\n", " plate = create_shape(shape, voxel_widths, thickness, shape_type=\"plate\")\n", - " est_plt_thicknesses.append(calc_structure_thickness_statistics(plate, voxel_widths, 0)[0])" + " est_plt_thicknesses.append(\n", + " calc_structure_thickness_statistics(plate, voxel_widths, 0)[0]\n", + " )" ] }, { @@ -204,7 +278,7 @@ } ], "source": [ - "max_size = max(shape)//2\n", + "max_size = max(shape) // 2\n", "\n", "fig, axs = plt.subplots(1, 3, figsize=(15, 5))\n", "\n", @@ -231,6 +305,25 @@ "\n", "plt.show()" ] + }, + { + "cell_type": "markdown", + "id": "242db030", + "metadata": {}, + "source": [ + "---\n", + "\n", + "\n", + "Notebook created using the [template](https://github.com/ORMIRcommunity/templates/blob/main/ORMIR_nb_template.ipynb) of the [ORMIR community](https://ormircommunity.github.io/) (version 1.0, 2023)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6d2417b9", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/manuscript/paper.md b/manuscript/paper.md index 8c5008a..563ceda 100644 --- a/manuscript/paper.md +++ b/manuscript/paper.md @@ -19,7 +19,7 @@ authors: - name: Andrew J. Burghardt orcid: 0000-0002-6343-4944 affiliation: 4 - - name: Michelle A.E. Espinosa + - name: Michelle A. Espinosa Hernandez orcid: 0000-0002-8279-9842 affiliation: "3, 5, 6" - name: Jared Vicory @@ -51,13 +51,13 @@ affiliations: index: 2 - name: Department of Biomedical Engineering, The University of Melbourne, Parkville, Australia index: 3 - - name: Department of Radiology and Biomedical Imaging, University of California, San Francisco, USA + - name: Department of Radiology and Biomedical Imaging, University of California, San Francisco, United States of America index: 4 - name: Rehabilitation Sciences Institute, The University of Toronto, Toronto, Canada index: 5 - name: Joint Department of Medical Imaging, University Health Network, Toronto, Canada index: 6 - - name: Kitware, Inc., Carrboro, North Carolina, USA + - name: Kitware, Inc., Carrboro, North Carolina, United States of America index: 7 - name: Department of Radiology, Cumming School of Medicine, University of Calgary, Calgary, Canada index: 8 @@ -69,44 +69,48 @@ date: 30 August 2023 bibliography: paper.bib --- -# 1 – Statement of need: +# Statement of need High resolution peripheral quantitative computed tomography (HR-pQCT) is an imaging technique with a nominal isotropic voxel size of 61µm capable of imaging trabecular bone *in-vivo*. HR-pQCT has a wide range of applications, primarily focused on bone to improve our understanding of musculoskeletal diseases [@Brunet:2020], assess epidemiological associations [@Burt:2016], and evaluate the effects of pharmaceutical interventions [@Brunet:2021]. Processing HR-pQCT images has largely been supported using the scanner manufacturer’s scripting language (Image Processing Language, IPL, Scanco Medical). However, by expanding image processing workflows outside of the scanner manufacturer’s software environment, users have the flexibility to apply more advanced mathematical techniques and leverage modern software packages to improve image processing. The `ORMIR_XCT` Python package was developed to reimplement some existing IPL workflows and provide an open and reproducible package allowing for the development of advanced HR-pQCT data processing workflows. -# 2 – Package summary: +# Package summary The development of this package began during the Jupyter Community Workshop in Maastricht, Netherlands in June 2022, hosted by the Open and Reproducible Musculoskeletal Imaging Research (ORMIR) group. During this workshop, the conceptualization and initial development of the `ORMIR_XCT` package began with support from Kitware. The `ORMIR_XCT` package is currently being used by members of the ORMIR community to assess bone mineral density (BMD) and joint space width (JSW) changes in HR-pQCT scans of patients with osteoarthritis or rheumatoid arthritis. The `ORMIR_XCT` package contains four modules for processing HR-pQCT data of bones and joints: 1) automatic contouring of the periosteal surface, 2) JSW analysis, 3) BMD calculation, and 4) segmentation of trabecular bone. Analyses have been performed to compare outputs from the `ORMIR_XCT` package to results generated using IPL. Jupyter Notebook examples are provided to describe the workflows implemented in the `ORMIR_XCT` package. -# 3 – Comparison to IPL: -## 3.1 - Autocontour: -Automatic periosteal contouring (autocontour) was performed on a sample dataset of HR-pQCT images of the 2nd and 3rd distal interphalangeal (DIP) and trapeziometacarpal (TMC) joints (n = 59). The images came from 10 participants in an osteoarthritic group and 10 participants in a control group. This dataset is representative of typical data that is collected using HR-pQCT. Binarized images were obtained using both the IPL standard method for automatic periosteal contouring [@Burghardt:2010] and the `ORMIR_XCT` package. Segmentation were compared by computing DICE coefficients, Jaccard indices, and Hausdorff distances between images. The IPL method for automatic periosteal contouring was taken as the ground truth. Results of the segmentation comparison are shown in \autoref{fig:Table1}. +# Comparison to IPL +## Autocontour +Automatic periosteal contouring (autocontour) was performed on a sample dataset of HR-pQCT images of the 2nd and 3rd distal interphalangeal (DIP) and trapeziometacarpal (TMC) joints (n = 59). The images came from 10 participants in an osteoarthritic group and 10 participants in a control group. This dataset is representative of typical data that is collected using HR-pQCT. Binarised images were obtained using both the IPL standard method for automatic periosteal contouring [@Burghardt:2010] and the `ORMIR_XCT` package. Segmentation were compared by computing DICE coefficients, Jaccard indices, and Hausdorff distances between images. The IPL method for automatic periosteal contouring was taken as the ground truth. DICE coefficients and Jaccard indices range from 0.0 to 1.0, with 1.0 indicating perfect overlap between segmentations. Results of the segmentation comparison are shown in \autoref{fig:Table1}. -![Table 1: DICE coefficients, Jaccard indices, and Hausdorff distances (voxels) between IPL and `ORMIR_XCT` automatic contouring implementations. All scans were obtained from a second-generation HR-pQCT scanner (XtremeCT2, Scanco Medical) with an isotropic voxel size of 61µm. DIP2: 2nd distal interphalangeal joint, DIP3: 3rd distal interphalangeal joint, TMC: trapeziometacarpal joint.\label{fig:Table1}](figures/Table1.png) +![DICE coefficients, Jaccard indices, and Hausdorff distances (voxels) between IPL and `ORMIR_XCT` automatic contouring implementations. All scans were obtained from a second-generation HR-pQCT scanner (XtremeCT2, Scanco Medical) with an isotropic voxel size of 61µm. DICE coefficients and Jaccard indices range from 0.0 to 1.0, and 1.0 indicates perfect overlap between segmentations. DIP2: 2nd distal interphalangeal joint, DIP3: 3rd distal interphalangeal joint, TMC: trapeziometacarpal joint.\label{fig:Table1}](figures/Table1.png) -## 3.2 - Thickness: +## Thickness The `ORMIR_XCT` package contains an open-source implementation of the algorithm for model-independent thickness estimation by Hildebrand *et al* [@Hildebrand:1996]. This thickness estimation is used to compute several morphological measurements from HR-pQCT images, including trabecular thickness, trabecular separation, and cortical thickness. Given a binary mask describing a three-dimensional object as input, the algorithm first computes the distance map (via a distance transform), describing at each voxel the radius of the largest sphere centered at that voxel that is entirely within the object. Next, the local thickness map is computed by iterating through all voxels in the mask and assigning them a local thickness value that is the diameter of the largest sphere in the distance map that contains that voxel. The thickness of a structure can then be computed as the volume-weighted mean of the local thickness map, with optional adjustments for a minimum thickness value if desired to reduce the effect of surface noise. -@Hildebrand:1996 note an optimization for the algorithm where first, the distance ridge (also known as skeletonization or medial axis) is computed and only the distance map values on the distance ridge are used for local thickness map computation. We implemented this optimized version of the algorithm with NumPy (v1.23.3) [@Harris:2020] and Numba (v0.56.3) [@Lam:2015]. Finally, we have noted and compensated for a discrepancy between the definition of the distance map by @Hildebrand:1996 and the outputs of many common open-source distance map algorithms (including those available in SimpleITK (v2.0.2) [@Lowekamp:2013] and scikit-image (v0.19.2) [@Van-der-Walt:2014]). @Hildebrand:1996 define their distance map as "the radius of the largest sphere centered at the point and still completely inside the structure". In contrast, most other distance transforms compute the distance from the center of a voxel inside the mask to the nearest voxel center outside the mask. +@Hildebrand:1996 note an optimisation for the algorithm where first, the distance ridge (also known as skeletonization or medial axis) is computed and only the distance map values on the distance ridge are used for local thickness map computation. We implemented this optimised version of the algorithm with NumPy (v1.23.3) [@Harris:2020] and Numba (v0.56.3) [@Lam:2015]. Finally, we have noted and compensated for a discrepancy between the definition of the distance map by @Hildebrand:1996 and the outputs of many common open-source distance map algorithms (including those available in SimpleITK (v2.0.2) [@Lowekamp:2013] and scikit-image (v0.19.2) [@Van-der-Walt:2014]). @Hildebrand:1996 define their distance map as "the radius of the largest sphere centered at the point and still completely inside the structure". In contrast, most other distance transforms compute the distance from the center of a voxel inside the mask to the nearest voxel center outside the mask. To compensate for this, we developed an oversampling distance transform function (\autoref{fig:DT_Figure}). In this function, an input mask is up-sampled to double the voxel size, and a sequence of morphological filters are applied to ensure that the newly introduced voxels on the boundary between the structure and the background are assigned as background voxels. A distance transform is applied to this mask, resulting in a distance map where the spheres fit entirely within the structure defined by the mask, and the resulting distance map is down-sampled to the original voxel size. -![Figure 1: Examples of the thickness estimation algorithm implemented in the `ORMIR_XCT` package. Top row: thickness calculation on a 1 voxel thick plate demonstrating how thickness estimation is improved using the oversampling distance transform function (thickness = 1 voxel using oversampling, thickness = 2 voxels without oversampling). Bottom row: thickness calculation applied to a sample finger joint.\label{fig:DT_Figure}](figures/DT_Figure.png) +![Examples of the thickness estimation algorithm implemented in the `ORMIR_XCT` package. Top row: thickness calculation on a 1 voxel thick plate demonstrating how thickness estimation is improved using the oversampling distance transform function (thickness = 1 voxel using oversampling, thickness = 2 voxels without oversampling). Bottom row: thickness calculation applied to a sample finger joint.\label{fig:DT_Figure}](figures/DT_Figure.png) The user is given a choice, via a Boolean flag parameter in the compute_local_thickness_from_mask and calc_structure_thickness_statistics functions, to use the oversampling distance transform for thickness estimation. The choice is given because different users may have variable preferences in applying alternate definitions regarding the bounds of their structures (whether they be contained within the mask voxels or extend to the center of the neighbouring background voxels) and because using the normal (*i.e.*, not oversampled) distance transform will provide better congruence with thickness calculations produced using IPL. To evaluate the thickness computation implemented in the `ORMIR_XCT` package, a set of synthetic shapes were generated to simulate bones of varying sizes. These shapes include solid and hollow spheres, solid and hollow cylinders, and plates of varying thicknesses [@ORMIR_XCT_Shapes:2023]. The standard IPL joint space width (JSW) analysis workflow [@Stok:2020] was performed on the synthetic shape dataset using the following parameters: `ridge_epsilon = 0.9`, `assign_epsilon = 1.8`, `suppress_boundary = 0`, and `version = 3`. This set of parameters is currently used for JSW analysis of *in vivo* finger joints. ORMIR_XCT thickness was calculated using the oversampling distance transform function described above. Bland-Altman plots were generated to compare thickness results between algorithms (\autoref{fig:IPL_vs_ORMIR_XCT-DT}). -![Figure 2: Bland-Altman and regression plots of mean joint space thickness computed using IPL and the `ORMIR_XCT` package . Top row: black lines indicate mean difference and red dotted line shows the limits of agreement. Bottom row: solid black line indicates the linear regression between IPL and ORMIR_XCT computed thickness and red dotted lines indicate a line of unity.\label{fig:IPL_vs_ORMIR_XCT-DT}](figures/IPL_vs_ORMIR_XCT-DT.png) +![Bland-Altman and regression plots of mean joint space thickness computed using IPL and the `ORMIR_XCT` package . Top row: black lines indicate mean difference and red dotted line shows the limits of agreement. Bottom row: solid black line indicates the linear regression between IPL and ORMIR_XCT computed thickness and red dotted lines indicate a line of unity.\label{fig:IPL_vs_ORMIR_XCT-DT}](figures/IPL_vs_ORMIR_XCT-DT.png) -## 3.3 - Bone mineral density: -Images obtained from HR-pQCT scanners are saved using the manufacturer’s proprietary file format (AIM). To process HR-pQCT images using Python, these files need to be converted to a file format that can be interpreted by common image processing Python libraries. The `ITKIOScanco` module is used in the `ORMIR_XCT` package to convert HR-pQCT images to other common medical image file types. However, the `ITKIOScanco` module converts HR-pQCT images to Hounsfield Units (HU) by default. Measures of BMD require the image to be calibrated in density units (mg of hydroxyapatite per cm3). Several image unit conversions have been included in the `ORMIR_XCT` package to allow for conversion between HR-pQCT native units (termed Scanco units) and each of HU, density, and linear attenuation units. A sample dataset of HR-pQCT images of second and third metacarpal phalangeal (MCP) joints (n = 292) was used to compare BMD results between IPL and the `ORMIR_XCT` package. BMD was reported separately for the distal (DST) and proximal (PRX) segments of each joint (\autoref{fig:IPL_vs_ORMIR_XCT-BMD}). +## Bone mineral density +Images obtained from HR-pQCT scanners are saved using the manufacturer’s proprietary file format (AIM). To process HR-pQCT images using Python, these files need to be converted to a file format that can be interpreted by common image processing Python libraries. The `ITKIOScanco` module is used in the `ORMIR_XCT` package to convert HR-pQCT images to other common medical image file types. However, the `ITKIOScanco` module converts HR-pQCT images to Hounsfield Units (HU) by default. Measures of BMD require the image to be calibrated in density units (mg of hydroxyapatite per cm3). Several image unit conversions have been included in the `ORMIR_XCT` package to allow for conversion between HR-pQCT native units (termed Scanco units) and each of HU, density, and linear attenuation units. A sample dataset of HR-pQCT images of second and third metacarpal phalangeal (MCP) joints (n = 292) was used to compare BMD results between IPL and the `ORMIR_XCT` package. BMD was reported separately for the distal (DST) and proximal (PRX) segments of each joint (\autoref{fig:IPL_vs_ORMIR_XCT-BMD}). -![Figure 3: Bone mineral density (mg HA/cm3) linear regression and Bland-Altman plots comparing IPL and the ORMIR_XCT computed bone mineral density (BMD ). All scans were obtained from a second-generation HR-pQCT scanner (XtremeCT2, Scanco Medical) with an isotropic voxel size of 61µm. MCP2: 2nd metacarpal phalangeal joint, MCP3: 3rd metacarpal phalangeal joint, DST: distal, PRX: proximal.\label{fig:IPL_vs_ORMIR_XCT-BMD}](figures/IPL_vs_ORMIR_XCT-BMD.png) +![Bone mineral density (mg HA/cm3) linear regression and Bland-Altman plots comparing IPL and the ORMIR_XCT computed bone mineral density (BMD ). All scans were obtained from a second-generation HR-pQCT scanner (XtremeCT2, Scanco Medical) with an isotropic voxel size of 61µm. MCP2: 2nd metacarpal phalangeal joint, MCP3: 3rd metacarpal phalangeal joint, DST: distal, PRX: proximal.\label{fig:IPL_vs_ORMIR_XCT-BMD}](figures/IPL_vs_ORMIR_XCT-BMD.png) -## 3.4 - Trabecular segmentation: +## Trabecular segmentation A common operation performed in IPL is the application of a Gaussian smoothing filter followed by a global threshold-based image segmentation. The `gauss_seg` function in IPL has been reimplemented in the `ORMIR_XCT` package and the similarity of IPL and `ORMIR_XCT` trabecular segmentations were computed using DICE coefficients, Jaccard indices, and Hausdorff distances (\autoref{fig:Table2}). The same dataset used to test the automatic contouring module was used to compare IPL and `ORMIR_XCT` trabecular segmentation. -![Table 2: DICE scores, Jaccard indices, and Hausdorff distances (voxels) for comparison of trabecular segmentation performed using IPL versus the `ORMIR_XCT` package. All scans were obtained from a second-generation HR-pQCT scanner (XtremeCT2, Scanco Medical) with an isotropic voxel size of 61µm. For each algorithm, Gaussian smoothing was performed with sigma = 0.5. Subsequently, segmentation was performed with lower and upper thresholds of 1,170 and 10,000 HU, respectively. DIP2: second distal interphalangeal joint, DIP3: third distal interphalangeal joint, TMC: trapeziometacarpal joint.\label{fig:Table2}](figures/Table2.png) +![DICE scores, Jaccard indices, and Hausdorff distances (voxels) for comparison of trabecular segmentation performed using IPL versus the `ORMIR_XCT` package. All scans were obtained from a second-generation HR-pQCT scanner (XtremeCT2, Scanco Medical) with an isotropic voxel size of 61µm. For each algorithm, Gaussian smoothing was performed with sigma = 0.5. Subsequently, segmentation was performed with lower and upper thresholds of 1,170 and 10,000 HU, respectively. DIP2: second distal interphalangeal joint, DIP3: third distal interphalangeal joint, TMC: trapeziometacarpal joint.\label{fig:Table2}](figures/Table2.png) -# References: +# Acknowledgements + +M.T.K. was supported by a Natural Sciences and Engineering Research Council (NSERC, Canada) fellowship. N.J.N. was supported by a Canadian Institutes of Health Research grant, Alberta Graduate Excellence Scholarship, and The Arthritis Society Canada Training Graduate Ph.D. Salary Award. S.L.M. was supported by The Arthritis Society Canada Stars Career Development Award. The ORMIR workshop where this package was developed was sponsored by NumFOCUS and financially supported by Bloomberg and Amazon Web Services. + +# References diff --git a/ormir_xct/joint_space_analysis/jsw_morphometry.py b/ormir_xct/joint_space_analysis/jsw_morphometry.py index c5b4701..7254473 100644 --- a/ormir_xct/joint_space_analysis/jsw_morphometry.py +++ b/ormir_xct/joint_space_analysis/jsw_morphometry.py @@ -36,7 +36,7 @@ def jsw_pad(joint_seg_image): Parameters ---------- - joint_seg_image : string + joint_seg_image : SimpleITK.Image Returns ------- @@ -54,6 +54,17 @@ def jsw_pad(joint_seg_image): joint_seg_image, [MISC2, MISC2, 0], [MISC2, MISC2, 0], 0 ) pad_image = sitk.BinaryThreshold(pad_image, 1, 127, 60, 0) + size = [joint_seg_image.GetSize()[0]+ MISC2*2, joint_seg_image.GetSize()[1]+ MISC2*2, joint_seg_image.GetSize()[2]] + origin = joint_seg_image.GetOrigin() + spacing = joint_seg_image.GetSpacing() + direction = joint_seg_image.GetDirection() + pad_image = sitk.Resample(pad_image, + size, + interpolator=sitk.sitkNearestNeighbor, + outputOrigin=origin, + outputSpacing=spacing, + outputDirection=direction + ) return pad_image @@ -115,7 +126,11 @@ def jsw_erode(dilated_image, pad_image): Returns ------- + eroded_image : SimpleITK.Image + js_mask : SimpleITK.Image + + dilated_js_mask : SimpleITK.Image """ # Erode the image, set the eroded mask's value to 30 eroded_image = sitk.BinaryErode( diff --git a/tests/autocontour/test_autocontour_gobj.py b/tests/autocontour/test_autocontour_gobj.py deleted file mode 100644 index e65cb8a..0000000 --- a/tests/autocontour/test_autocontour_gobj.py +++ /dev/null @@ -1,20 +0,0 @@ -""" -test_autocontour.py - -Created by: Michael Kuczynski -Created on: June 29, 2022 - -Description: Test automatic periosteal contour. -""" - -import unittest -import numpy as np -import SimpleITK as sitk - -from ormir_xct.autocontour.autocontour_gobj import autocontour_gobj -from ormir_xct.autocontour.AutocontourKnee import AutocontourKnee - - -class TestAutocontour(unittest.TestCase): - def test_autocontour(self): - pass diff --git a/tests/data/test_dicom/IMG0001.dcm b/tests/data/test_dicom/IMG0001.dcm new file mode 100644 index 0000000..994b448 Binary files /dev/null and b/tests/data/test_dicom/IMG0001.dcm differ diff --git a/tests/data/test_dicom/IMG0002.dcm b/tests/data/test_dicom/IMG0002.dcm new file mode 100644 index 0000000..e3a0871 Binary files /dev/null and b/tests/data/test_dicom/IMG0002.dcm differ diff --git a/tests/data/test_dicom/IMG0003.dcm b/tests/data/test_dicom/IMG0003.dcm new file mode 100644 index 0000000..db6b55c Binary files /dev/null and b/tests/data/test_dicom/IMG0003.dcm differ diff --git a/tests/data/test_dicom/IMG0004.dcm b/tests/data/test_dicom/IMG0004.dcm new file mode 100644 index 0000000..3d7b14a Binary files /dev/null and b/tests/data/test_dicom/IMG0004.dcm differ diff --git a/tests/data/test_dicom/IMG0005.dcm b/tests/data/test_dicom/IMG0005.dcm new file mode 100644 index 0000000..e2e9a60 Binary files /dev/null and b/tests/data/test_dicom/IMG0005.dcm differ diff --git a/tests/data/test_dicom/IMG0006.dcm b/tests/data/test_dicom/IMG0006.dcm new file mode 100644 index 0000000..e1b5089 Binary files /dev/null and b/tests/data/test_dicom/IMG0006.dcm differ diff --git a/tests/data/test_dicom/IMG0007.dcm b/tests/data/test_dicom/IMG0007.dcm new file mode 100644 index 0000000..28617a1 Binary files /dev/null and b/tests/data/test_dicom/IMG0007.dcm differ diff --git a/tests/data/test_dicom/IMG0008.dcm b/tests/data/test_dicom/IMG0008.dcm new file mode 100644 index 0000000..5ac48b9 Binary files /dev/null and b/tests/data/test_dicom/IMG0008.dcm differ diff --git a/tests/data/test_dicom/IMG0009.dcm b/tests/data/test_dicom/IMG0009.dcm new file mode 100644 index 0000000..6cb6249 Binary files /dev/null and b/tests/data/test_dicom/IMG0009.dcm differ diff --git a/tests/data/test_dicom/IMG0010.dcm b/tests/data/test_dicom/IMG0010.dcm new file mode 100644 index 0000000..42b8fd1 Binary files /dev/null and b/tests/data/test_dicom/IMG0010.dcm differ diff --git a/tests/data/test_dicom/IMG0011.dcm b/tests/data/test_dicom/IMG0011.dcm new file mode 100644 index 0000000..f8877ec Binary files /dev/null and b/tests/data/test_dicom/IMG0011.dcm differ diff --git a/tests/data/test_dicom/IMG0012.dcm b/tests/data/test_dicom/IMG0012.dcm new file mode 100644 index 0000000..a2ce045 Binary files /dev/null and b/tests/data/test_dicom/IMG0012.dcm differ diff --git a/tests/data/test_dicom/IMG0013.dcm b/tests/data/test_dicom/IMG0013.dcm new file mode 100644 index 0000000..fa9c23e Binary files /dev/null and b/tests/data/test_dicom/IMG0013.dcm differ diff --git a/tests/data/test_dicom/IMG0014.dcm b/tests/data/test_dicom/IMG0014.dcm new file mode 100644 index 0000000..1a4616a Binary files /dev/null and b/tests/data/test_dicom/IMG0014.dcm differ diff --git a/tests/data/test_dicom/IMG0015.dcm b/tests/data/test_dicom/IMG0015.dcm new file mode 100644 index 0000000..2ab4fd0 Binary files /dev/null and b/tests/data/test_dicom/IMG0015.dcm differ diff --git a/tests/data/test_dicom/IMG0016.dcm b/tests/data/test_dicom/IMG0016.dcm new file mode 100644 index 0000000..d757dd5 Binary files /dev/null and b/tests/data/test_dicom/IMG0016.dcm differ diff --git a/tests/data/test_dicom/IMG0017.dcm b/tests/data/test_dicom/IMG0017.dcm new file mode 100644 index 0000000..bf524f3 Binary files /dev/null and b/tests/data/test_dicom/IMG0017.dcm differ diff --git a/tests/data/test_dicom/IMG0018.dcm b/tests/data/test_dicom/IMG0018.dcm new file mode 100644 index 0000000..ff3e425 Binary files /dev/null and b/tests/data/test_dicom/IMG0018.dcm differ diff --git a/tests/data/test_dicom/IMG0019.dcm b/tests/data/test_dicom/IMG0019.dcm new file mode 100644 index 0000000..b5399cb Binary files /dev/null and b/tests/data/test_dicom/IMG0019.dcm differ diff --git a/tests/data/test_dicom/IMG0020.dcm b/tests/data/test_dicom/IMG0020.dcm new file mode 100644 index 0000000..b696595 Binary files /dev/null and b/tests/data/test_dicom/IMG0020.dcm differ diff --git a/tests/data/test_dicom/IMG0021.dcm b/tests/data/test_dicom/IMG0021.dcm new file mode 100644 index 0000000..56a402c Binary files /dev/null and b/tests/data/test_dicom/IMG0021.dcm differ diff --git a/tests/data/test_dicom/IMG0022.dcm b/tests/data/test_dicom/IMG0022.dcm new file mode 100644 index 0000000..90ac4e7 Binary files /dev/null and b/tests/data/test_dicom/IMG0022.dcm differ diff --git a/tests/data/test_dicom/IMG0023.dcm b/tests/data/test_dicom/IMG0023.dcm new file mode 100644 index 0000000..6d374d2 Binary files /dev/null and b/tests/data/test_dicom/IMG0023.dcm differ diff --git a/tests/data/test_dicom/IMG0024.dcm b/tests/data/test_dicom/IMG0024.dcm new file mode 100644 index 0000000..5aa4f5f Binary files /dev/null and b/tests/data/test_dicom/IMG0024.dcm differ diff --git a/tests/data/test_dicom/IMG0025.dcm b/tests/data/test_dicom/IMG0025.dcm new file mode 100644 index 0000000..14ca4ff Binary files /dev/null and b/tests/data/test_dicom/IMG0025.dcm differ diff --git a/tests/data/test_joint.nii b/tests/data/test_joint.nii new file mode 100644 index 0000000..1b5fec7 Binary files /dev/null and b/tests/data/test_joint.nii differ diff --git a/tests/data/test_joint_mask.nii b/tests/data/test_joint_mask.nii new file mode 100644 index 0000000..4325933 Binary files /dev/null and b/tests/data/test_joint_mask.nii differ diff --git a/tests/joint_space_analysis/test_bone_connectivity.py b/tests/joint_space_analysis/test_bone_connectivity.py index 6798ae2..3a5212d 100644 --- a/tests/joint_space_analysis/test_bone_connectivity.py +++ b/tests/joint_space_analysis/test_bone_connectivity.py @@ -1,9 +1,47 @@ import unittest +import numpy as np +import SimpleITK as sitk from ormir_xct.joint_space_analysis.connected_check import connected_check +def create_sphere_mask(shape, voxel_width, radius): + """ + Creates a binary NumPy sphere with given radius and voxel size. + """ + center = ( + voxel_width[0] * (shape[0] // 2), + voxel_width[1] * (shape[1] // 2), + voxel_width[2] * (shape[2] // 2), + ) + x, y, z = np.meshgrid( + *[voxel_width[i] * np.arange(0, shape[i]) for i in range(3)], indexing="ij" + ) + mask = np.zeros(shape, dtype=bool) + mask[ + (x - center[0]) ** 2 + (y - center[1]) ** 2 + (z - center[2]) ** 2 < radius**2 + ] = 1 + return mask + + class TestBoneConnectivityCheck(unittest.TestCase): - @unittest.skip("unimplemented") - def test_placeholder(self): - pass + def test_connected_check(self): + voxel_width = (1, 1, 1) + shape = (30, 30, 30) + radius = 4 + sphere_array = create_sphere_mask(shape, voxel_width, radius).astype(float) + translated_sphere_x = np.roll(sphere_array, 5, axis=[1, 0, 0]) + + empty = sitk.Image(shape, sitk.sitkUInt8) + one_sphere = sitk.GetImageFromArray(sphere_array) + two_sphere = sitk.GetImageFromArray(sphere_array + translated_sphere_x) + one_sphere = sitk.Cast(one_sphere, sitk.sitkUInt8) + two_sphere = sitk.Cast(two_sphere, sitk.sitkUInt8) + + empty_labels = connected_check(empty) + one_sphere_labels = connected_check(one_sphere) + two_sphere_labels = connected_check(two_sphere) + + self.assertAlmostEqual(empty_labels, 0, 8) + self.assertAlmostEqual(one_sphere_labels, 1, 8) + self.assertAlmostEqual(two_sphere_labels, 2, 8) diff --git a/tests/joint_space_analysis/test_jsw_main.py b/tests/joint_space_analysis/test_jsw_main.py deleted file mode 100644 index d0ab33c..0000000 --- a/tests/joint_space_analysis/test_jsw_main.py +++ /dev/null @@ -1,9 +0,0 @@ -import unittest - -from ormir_xct.joint_space_analysis.jsw_main import main - - -class TestJointSpaceWidthMain(unittest.TestCase): - @unittest.skip("unimplemented") - def test_placeholder(self): - pass diff --git a/tests/joint_space_analysis/test_jsw_morphometry.py b/tests/joint_space_analysis/test_jsw_morphometry.py index 7f56f75..ccb407b 100644 --- a/tests/joint_space_analysis/test_jsw_morphometry.py +++ b/tests/joint_space_analysis/test_jsw_morphometry.py @@ -1,14 +1,188 @@ +import os +import shutil +import tempfile import unittest +import numpy as np +import SimpleITK as sitk from ormir_xct.joint_space_analysis.jsw_morphometry import ( jsw_pad, jsw_dilate, jsw_erode, jsw_parameters, + MISC2, + CALC, ) +def create_plate_mask(shape, voxel_width, radius): + """ + Creates a binary NumPy plate with given radius and voxel size. + """ + center = ( + voxel_width[0] * (shape[0] // 2), + voxel_width[1] * (shape[1] // 2), + voxel_width[2] * (shape[2] // 2), + ) + x, y, z = np.meshgrid( + *[voxel_width[i] * np.arange(0, shape[i]) for i in range(3)], indexing="ij" + ) + mask = np.zeros(shape, dtype=bool) + mask = (np.abs(x - center[0]) < radius / 2).astype(int) + return mask + + class TestJointSpaceWidthMorphometry(unittest.TestCase): - @unittest.skip("unimplemented") - def test_placeholder(self): - pass + def spacing_check(self, input_image, ref_image): + input_spacing = input_image.GetSpacing() + output_spacing = ref_image.GetSpacing() + + self.assertAlmostEqual(input_spacing[0], output_spacing[0], 4) + self.assertAlmostEqual(input_spacing[1], output_spacing[1], 4) + self.assertAlmostEqual(input_spacing[2], output_spacing[2], 4) + + def dimension_check(self, input_image, ref_image): + input_size = input_image.GetSize() + output_size = ref_image.GetSize() + + self.assertAlmostEqual(input_size[0], output_size[0], 4) + self.assertAlmostEqual(input_size[1], output_size[1], 4) + self.assertAlmostEqual(input_size[2], output_size[2], 4) + + def position_check(self, input_image, ref_image): + input_origin = input_image.GetOrigin() + output_origin = ref_image.GetOrigin() + + self.assertAlmostEqual(input_origin[0], output_origin[0], 4) + self.assertAlmostEqual(input_origin[1], output_origin[1], 4) + self.assertAlmostEqual(input_origin[2], output_origin[2], 4) + + def direction_check(self, input_image, ref_image): + input_direction = input_image.GetDirection() + output_direction = ref_image.GetDirection() + + self.assertAlmostEqual(input_direction[0], output_direction[0], 4) + self.assertAlmostEqual(input_direction[1], output_direction[1], 4) + self.assertAlmostEqual(input_direction[2], output_direction[2], 4) + + def test_jsw_pad(self): + voxel_width = (1, 1, 1) + shape = (30, 30, 30) + radius = 4 + plate_array = create_plate_mask(shape, voxel_width, radius).astype(int) + + plate_image = sitk.GetImageFromArray(plate_array) + + pad_image = jsw_pad(plate_image) + + image_size = plate_image.GetSize() + pad_image_size = pad_image.GetSize() + + # Image is padded by MISC2 (27) on each side in the X and Y direction. + # I.e., 54 in X and 54 in Y + plate_x = image_size[0] + plate_y = image_size[1] + plate_z = image_size[2] + pad_x = pad_image_size[0] + pad_y = pad_image_size[1] + pad_z = pad_image_size[2] + + self.assertAlmostEqual(plate_x + MISC2 * 2, pad_x, 8) + self.assertAlmostEqual(plate_y + MISC2 * 2, pad_y, 8) + self.assertAlmostEqual(plate_z, pad_z, 8) + + self.spacing_check(pad_image, plate_image) + self.position_check(pad_image, plate_image) + self.direction_check(pad_image, plate_image) + + def test_jsw_dilate(self): + voxel_width = (1, 1, 1) + shape = (30, 30, 30) + radius = 4 + + # Image is dilated using a kernel of [MISC2, MISC2, MISC2], where MISC2 = 27 + plate_array_1 = create_plate_mask(shape, voxel_width, radius).astype(int) + plate_array_2 = create_plate_mask(shape, voxel_width, radius + MISC2).astype( + int + ) + + plate_1 = sitk.GetImageFromArray(plate_array_1) + plate_2 = sitk.GetImageFromArray(plate_array_2) + + dilated = jsw_dilate(plate_1) + dilated = sitk.BinaryThreshold(dilated, 127, 127, 1, 0) + dilated_array = sitk.GetArrayFromImage(dilated) + + np.testing.assert_array_equal(plate_array_2, dilated_array) + self.spacing_check(dilated, plate_2) + self.position_check(dilated, plate_2) + self.direction_check(dilated, plate_2) + + def test_jsw_erode(self): + voxel_width = (1, 1, 1) + shape = (200, 200, 200) + radius = CALC + + # Image is eroded using a kernel of [CALC, CALC, CALC], where CALC = 35 + plate_array_1 = create_plate_mask(shape, voxel_width, 3 * radius).astype(int) + plate_array_2 = create_plate_mask(shape, voxel_width, radius).astype(int) + + plate_1 = sitk.GetImageFromArray(plate_array_1) + plate_2 = sitk.GetImageFromArray(plate_array_2) + + plate_1 = sitk.BinaryThreshold(plate_1, 1, 1, 127, 0) + eroded_image, js_mask, dilated_js_mask = jsw_erode(plate_1, jsw_pad(plate_1)) + eroded_image = sitk.BinaryThreshold(eroded_image, 30, 30, 1, 0) + eroded_array = sitk.GetArrayFromImage(eroded_image) + + np.testing.assert_array_equal(plate_array_2, eroded_array) + self.spacing_check(eroded_image, plate_2) + self.position_check(eroded_image, plate_2) + self.direction_check(eroded_image, plate_2) + + def test_jsw_parameters(self): + voxel_width = (1, 1, 1) + shape = (100, 100, 100) + radius = 10 + plate_array = create_plate_mask(shape, voxel_width, radius).astype(int) + + plate_array = 1 - plate_array + plate_image = sitk.GetImageFromArray(plate_array) + + # Step 1: Image padding + pad = jsw_pad(plate_image) + + # Step 2: Image dilation + dilate = jsw_dilate(pad) + + # Step 3: Image erosion + erode, js_mask, dilated_js_mask = jsw_erode(dilate, pad) + + # Step 4: Compute parameters + test_dir = tempfile.mkdtemp() + test_csv = os.path.join(test_dir, "jsw_params.csv") + + dt_img, jsw_params = jsw_parameters( + pad, + dilated_js_mask, + test_csv, + test_dir, + js_mask=js_mask, + voxel_size=1, + oversamp=False, + skel=False, + minimum=0.0, + ) + shutil.rmtree(test_dir) + + space = js_mask.GetSpacing() + voxel = np.prod(space) + img = sitk.GetArrayFromImage(js_mask) + vol = voxel * np.sum(img) + + self.assertAlmostEqual(jsw_params[0][2], vol, 8) + self.assertAlmostEqual(jsw_params[0][3], radius, 8) + self.assertAlmostEqual(jsw_params[0][4], 0, 8) + self.assertAlmostEqual(jsw_params[0][5], radius, 8) + self.assertAlmostEqual(jsw_params[0][7], radius, 8) + self.assertAlmostEqual(jsw_params[0][8], radius / radius, 8) diff --git a/tests/util/test_hildebrand_thickness.py b/tests/util/test_hildebrand_thickness.py index e6f5663..7f3b6b1 100644 --- a/tests/util/test_hildebrand_thickness.py +++ b/tests/util/test_hildebrand_thickness.py @@ -25,12 +25,6 @@ def create_sphere_mask(shape, voxel_width, radius): return mask -class TestComputeLocalThicknessFromMask(unittest.TestCase): - @unittest.skip("unimplemented") - def test_placeholder(self): - pass - - class TestCalcStructureThicknessStatistics(unittest.TestCase): def test_get_sphere_thickness(self): voxel_width = (1, 1, 1) @@ -39,7 +33,7 @@ def test_get_sphere_thickness(self): min_thickness = 0 sphere = create_sphere_mask(shape, voxel_width, radius) s = calc_structure_thickness_statistics( - sphere, voxel_width, min_thickness, pad_amount=10, oversample=False + sphere, voxel_width, min_thickness, oversample=False ) self.assertAlmostEqual(2 * radius, s[0], places=1) diff --git a/tests/util/test_vtk_sitk_conversion.py b/tests/util/test_vtk_sitk_conversion.py index c6c0169..851f7a4 100644 --- a/tests/util/test_vtk_sitk_conversion.py +++ b/tests/util/test_vtk_sitk_conversion.py @@ -20,7 +20,7 @@ def create_test_image(dimensions, value): """ Creates a numpy array containing a single value of given dimension. """ - array = np.full(dimensions, value, np.int32) + array = np.full(dimensions, value, np.int16) return array