diff --git a/manuscript/paper.md b/manuscript/paper.md index 6aa4d37..563ceda 100644 --- a/manuscript/paper.md +++ b/manuscript/paper.md @@ -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 @@ -80,14 +80,14 @@ The `ORMIR_XCT` package contains four modules for processing HR-pQCT data of bon # 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. 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. 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}. +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}. ![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) ## 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. @@ -109,4 +109,8 @@ A common operation performed in IPL is the application of a Gaussian smoothing f ![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) +# 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