SPHARM-PDM based image preprocessing pipeline for quantitative morphometric analysis (QMA) for in situ joint assessment in rabbit and rat models

The accessibility of quantitative measurements of joint morphometry depends on appropriate tibial alignment and volume of interest (VOI) selection of joint compartments; often a challenging and time-consuming manual task. In this work, we developed a novel automatic, efficient, and model-invariant image preprocessing pipeline that allows for highly reproducible 3D quantitative morphometric analysis (QMA) of the joint. The pipeline addresses the problem by deploying two modules: an alignment module and a subdivision module. Alignment is achieved by representing the tibia in its basic form using lower degree spherical harmonic basis functions and aligning using principal component analysis. The second module subdivides the joint into lateral and medial VOIs via a watershedding approach based on persistence homology. Multiple repeated micro-computed tomography scans of small (rat) and medium (rabbit) animal knees were processed using the pipeline to demonstrate model invariance. Existing QMA was performed to evaluate the pipeline’s ability to generate reproducible measurements. Intraclass correlation coefficient and mean-normalised root-mean-squared error of more than 0.75 and lower than 9.5%, respectively, were achieved for joint centre of mass, joint contact area under virtual loading, joint space width, and joint space volume. Processing time and technical requirements were reduced compared to manual processing in previous studies.

its medial and lateral components for separate joint QMA measurements of each side. The common orientation is defined as the position where the long axis of the tibia is aligned with the vertical z-axis while the rest of the joint is rigidly transformed accordingly to preserve the original relative pose of all joint components. Previously, this was achieved by manually aligning one joint sample as a reference. All other samples were then registered and transformed onto that reference using B-spline interpolation 13 to reduce rotational errors 14 . Selection of the volume of interest (VOI) of the joint's medial and lateral sides was subsequently done manually by an expert through visual inspection of the 3D images according to appropriate anatomical features.
However, the manual nature of the initial alignment process, as well as the selection of the medial and lateral joint VOI has meant the quality of the joint QMA measurements is highly dependent on the skills and training of the operators to achieve high reproducibility and disease discriminating quality. Moreover, the diversity of animal models involved in preclinical OA research introduces a variety of imaging resolutions [15][16][17][18] . Consequently, settings for image processing operations such as dilation, erosion, closing, opening, and others, in the workflow must also be manually adjusted to achieve a reliable result. This is a non-trivial task that prohibits the use of a simple adaptation between models 19 , and is a time-consuming process that can occupy skilled operators for several weeks per study and present significant challenges to the robust application of joint QMA measurements in studies involving larger volumes of data and new joint models.
To tackle these challenges, this study presents an image processing pipeline that can automatically perform appropriate joint alignment and subdivision for efficient and high-quality 3D joint QMA while remaining robust across multiple small animal models. This study aims to estimate the shape and pose of the tibia and, by extension, the joint, based on a shape analysis technique called the spherical harmonic description method (SPHARM) 20 , which employs spherical harmonic descriptors in distinction to the traditional registration-based approach. Medial and lateral subdivision of the VOI is implemented by locating the dividing point using a watershedding approach based on persistent homology 21 , which is a topological data analysis method that identifies points with strong features. This study hypothesises that the pipeline will allow measurements of 3D joint QMA parameters with precision and reproducibility comparable to that of earlier studies that used manual alignment 11,12 , while achieving faster implementation and maintaining precision across species. Specifically, this study aims to compare measurement precision and reproducibility, as well as the overall processing time between the presented novel method and the manual processing approach on the same datasets which were published previously 11,12 .

Materials and methods
Animals and microCT scan protocols. Datasets were obtained from previous studies of intact ex vivo joints from rat 12 and rabbit 11 and are described in detail in the respective publications. The rat dataset consists of microCT scans (SCANCO Medical AG, Brüttisellen, Switzerland) of 21 tibio-femoral joints from 11 age-matched, Wistar rats. With a typical medial-lateral width of the tibia of 8.09 ± 0.62 mm, they were scanned with an isotropic nominal resolution of 10 µm. The rat dataset has an average image size in x, y, and z dimension of 1652.36 ± 221.39, 1416.59 ± 242, 943.5 ± 200.33 voxels and an average total number of voxels of 2.13 × 10 9 ± 3.09 × 10 8 , respectively. The rabbit dataset consists of microCT scans of 6 tibio-femoral joints from 6 age-matched, New Zealand white rabbits. With a typical medial-lateral width of the tibia around 19.10 ± 0.59 mm, they were obtained with an isotropic voxel size of 18 μm, average image size in the x, y, and z dimension of 2017.00 ± 43.14, 2017.00 ± 43.14, 1891.83 ± 462.15 voxels and an average total number of voxels of 7.65 × 10 9 ± 1.59 × 10 9 , respectively. During image acquisition of both datasets, initial joint positioning was controlled by placing a wedge (approx. 160° angle) behind the knee to control flexion-extension. The femoral condyles and the tibial plateau were included in the volume of interest, with the upper limit defined as the epiphyseal bone of the femur and the lower limit defined as the epiphyseal bone of the tibia (approx. 35-40 mm in rabbits and 12-15 mm in rats). Both datasets were filtered using a constrained 3D Gaussian filter (window, σ = 1.2, support, s = 1) after image reconstruction.
For each joint, four scans were performed. PRE scans were obtained by scanning joints without any contrast agent, while three repeat scans (labelled HEX1, HEX2, and HEX3) were obtained after a single intraarticular injection of SiO2-microbeads (0-20 μm diameter) (SWARCO Vestglass GmbH, Recklinghausen, Germany) to allow the visualisation of the cartilage. The three repeat scans were carried out to test for reproducibility with re-positioning between each scan.
Preprocessing pipeline for joint QMA. The developed image preprocessing pipeline consists of 2 main components: an alignment module and a subdivision module. Filtered and segmented 3D binarised images of the femur, tibia, and their respective cartilage volumes are required inputs for the workflow. These images served as inputs for the alignment module. The resulting aligned, binarised images are used as inputs for the subdivision module which is the final module of the workflow. The images output from the pipeline are aligned and split into medial and lateral VOI and are ready for submission to the existing joint QMA analysis modules. An overview of the pipeline and the joint QMA measurements used to evaluate each process, are summarised in Fig. 1.
Alignment module. Joint alignment method in this work estimated the tibia's elementary shape and pose using SPHARM 20 . In brief, SPHARM provides a hierarchical and multi-scale boundary description of objects with spherical topology. It computes an area-preserving mapping of the object's 3D voxel mesh onto a sphere in a separate parameter domain. The basis functions of the sphere are spherical harmonics, so that SPHARM describes the object as a set of basis function weights. By truncating the number of basis functions used in the description, different levels of object detail can be represented. where Y m * l (θ, ϕ) denotes the complex conjugate of Y m l and P m l describes the associated Legendre polynomials. Some of the low order real spherical harmonics as derived from the above equation are visualised in Supplementary Figure S1 online. This highlights the hierarchical representation of spherical objects where a higher spherical harmonics degree of l would lead to more complex forms of θ and ϕ and, thus, allows more details of the objects to be represented. Vice versa, by limiting the degree l to lower degree, the object detail would be reduced, up to the point where the object is represented by a sphere at l = 0. www.nature.com/scientificreports/ To express the object's surface using the described spherical harmonics, the three coordinate functions are decomposed and the surface V (θ, ϕ) = (x(θ , ϕ), y(θ, ϕ), z(θ, ϕ)) T takes the form of: where the coefficients c m l are 3D vectors of the three coordinates functions, obtained using a minimum squared error approach. Therefore, the values of the basis functions at each point in the discretised parameterised domain are gathered in the matrix, z = z i,j(l,m) = Y m l (θ i , ϕ i ) where j(l, m) is a function assigning an index to every pair (l, m) and i denotes the indices of n vert points to be approximated. The coordinates of these points are arranged in the vector v = (V 1 , V 2 , . . . , V n vert ) T and all coefficients are gathered in the vector c = (c 0 0 , c −1 1 , c 0 1 , . . . ) T . The coefficients that best approximate the points from a least-square perspective are obtained by: Using spherical harmonic basis functions, a hierarchical surface description is obtained that includes further details as more basis functions are added according to degree l and order m (see Fig. 2a,b).
After obtaining the tibia's basic form from its truncated SPHARM description, principal component analysis (PCA) was performed. The smallest component of the tibia, which always aligns with its vertical axis, was used www.nature.com/scientificreports/ as the surrogate for the tibial vertical orientation, while the horizontal orientation was controlled by the secondsmallest component. Subsequently, a series of transformation matrices aligning the smallest component with the z-axis and the second-smallest component with the x-axis was calculated, as seen in Fig. 2c, and was used to transform the full joint image using a B-spline interpolation method as the subsequent joint alignment step.
Subdivision module. Images aligned in the previous module were further processed through the subdivision module to divide the image of the joint into its medial and lateral VOI. The module operates by performing a watershedding operation based on persistence homology to identify the local minimum between the intercondylar tubercles of the tibia's intercondylar eminence, as indicated in Fig. 3a,b, and using that location as the dividing point for defining the tibia, and subsequently the joint, medial and lateral side. The module achieves this by analysing a series of projections created from the aligned microCT dataset and taking advantage of the controlled alignment achieved in the previous module, where the axial and coronal planes are aligned with the z-and x-axis, respectively. The subdivision module created a 2-dimensional (2D) topographic map of the tibial plateau's features by projecting along the coronal plane as seen in Fig. 3c. Subsequently, the map was used to create a 1-dimensional (1D) projection on the coronal plane as illustrated by the blue profile in Fig. 3d. To further limit the impact of osteophytes on the dividing point detection algorithm, the 1D projection was cropped to exclude values not in the centre of the object which is the area where the intercondylar tubercle is located.
From the processed 1D projection map, the dividing point was located using persistence homology, which is an algebraic method for capturing topological features like local extrema 23,24 , computed using the persistence1d 25 algorithm. In brief, the algorithm detects persistence homology from a 1D signal, which corresponds to all pairs www.nature.com/scientificreports/ of its local minima and maxima. Using this method, the global maxima is always paired with the background and has the highest persistence. The module's watershedding approach removes all persistence pairs except the second-largest one, whose maxima corresponds to the minimum point between intercondylar eminences; and thus, defined as the watershed point used to divide and select the sub-images.
3D joint QMA, reproducibility, and accuracy. The joint centre of mass, defined as a vector with orientation (α, degree; β, degree; and γ, degree) measured along the three principle Cartesian axes, and connecting the centres of mass of the two bones, was used to measure the relative position of the joint to evaluate the performance of the alignment module 11 . To assess the subdivision module, as well as the pipeline as a whole, several 3D joint metrics were measured from the resulting medial (M) and lateral (L) sub-images of the joint. Joint contact area under virtual loading was also measured. This was done by incrementally shifting the tibia onto the femur along the vertical z-axis. The contact area is, then, defined as the distance travelled to first contact ( χ L/M , μm) and the rate at which contact area increases ( m L/M , mm 2 step ) 11 . Alongside these metrics, 3D measurements of the joint space width (JSW): joint space volume ( JSV L/M , mm 3 ), JSW ( JSW L/M , mm), minimum JSW ( JSW L/M .min , mm), and maximum JSW ( JSW L/M .max , mm) were also measured 12,26 .
Reproducibility of the measured 3D QMA was tested using RStudio (RStudio: Integrated development environment for R, Version 0.95.258, Boston, MA, USA). Precision error (PE(SD)) and precision error expressed as coefficients of variation (PE(%CV)), as well as an intraclass correlation coefficient (two-way random effects, consistency, multiple measurements ICC 27 ) with 95% confidence interval (CI), were calculated to test the absolute agreement between each measurement 28 . ICC values range from 0 to 1, with those close to 1 indicating high similarity between measurements from the same group. Reproducibility of measurements is classified as excellent when ICC values are higher than 0.75 29 . Good, fair, and poor reproducibility is classified when ICC values range from 0.6 to 0.74, 0.40 to 0.59, and 0 to 0.4, respectively 29 .
Reproducibility results were compared with reproducibility values obtained in previous studies 11,12 using the same dataset but has performed the tasks in this workflow manually. Additionally, to measure accuracy, the joint QMA results obtained in this study were compared with results from earlier studies 11,12 . Measurement differences were evaluated as root-mean-squared error expressed in both absolute (RMSE) and mean-normalised (NRMSE) form by using the original results as reference.
Pipeline implementation and performance. The pipeline was implemented in two separate parts.
SPHARM processing (statistical shape analysis module 30 on 3DSlicer 31 ), alignment determination, and dividing point location were done on Windows 10 running on Intel i7-10875H processor (8 cores, 16 MB cache, 2.30 GHz to 5.10 GHz) with Nvidia Quadro P620 GPU and 16 GB RAM, while the subsequent joint alignment and subdivision were done on networked OpenVMS V8.4 I64 running on HP Integrity rx2800 i2 platform (8 Intel Itanium CPUs, 1.6 GHz, 5 MB cache) with 32 GB of allocated RAM. The source code is available from the corresponding author upon request.
Pipeline performance was recorded as the average CPU time (in seconds) required to perform each process. Processing time was measured to enable benchmarking against manual processing time in earlier studies (in hours) on the same dataset 11,12 whose estimates were obtained through examination of processing logs. Additionally, to benchmark the alignment module against another solution, the automatic joint alignment module was applied to the original images, without using SPHARM as part of the process.

Results
To select the optimal SPHARM degree for the workflow, a parametric study where the ICCs of the resulting joint centre of mass were evaluated for an increasing number of harmonics. The results are available in Table S1 of the supplementary document which shows that a SPHARM degree of 5 yielded the highest ICC result and was used for all the subsequent results in this manuscript.   Table 2). For rabbit data, all except m M with ICC of 0.841, has ICC of more than 0.9 ( Table 2). For all contact area measurements of both models, small precision errors are reported with PE(%CV) lower than 10% for all values.
For rat data, excellent reproducibility in was achieved for all joint space measurements, except JSW.max where the reproducibility is very low ( JSW.max L : 0.384, JSW.max M : 0.536). Among the rest, ICCs for JSW.min are slightly lower ( JSW.min L : 0.859, JSW.min M : 0.874) than others which have ICCs of more than 0.9, as shown in Table 3. Similarly, as seen in Table 4, excellent ICCs were also achieved for rabbit joint space measurements with values ranging from 0.780 (for JSW.min M ) to 0.955 (for JSV L ) except for JSW.max L , which have an ICC of 0.608. Small precision error is also reported for both models, with PE(%CV) lower than 6 from the manual measures of less than 9.5% (Table 4). For rabbit data, 3D joint space measurements were not performed, so that accuracy results are available only for the rat dataset which have differences from the manual measures of less than 9.5% (Table 4). JSV and JSW.max have particularly low differences; less than 3.5% and less than 2.5%, respectively (Table 4).
Pipeline performance. As seen in Table 5, the average time the pipeline used to finish processing a sample is shorter (rat: 693 s, rabbit: 2,112 s) than manual processing done on the same datasets in previous rat and rabbit studies (7,200 s for initial alignment and an additional 2,700 s for each rat sample and an additional 5,400 s for each rabbit sample). It is noted that a relatively high standard deviation was observed (rat: 212 s, rabbit: 276 s) due to networked connections. The most time-consuming process of the pipeline was the subsequent joint alignment using B-spline interpolation (rat: 248 ± 69 s, rabbit: 883 ± 261 s) and SPHARM processing (rat: 184 ± 77 s, rabbit: 1,019 ± 61 s) while the subdivision location was the fastest (rat: 19 ± 10 s, rabbit: 10 ± 3 s). Table 1. Reproducibility of the rat and rabbit joint centre of mass in terms of intraclass correlation coefficient (ICC) and precision errors (PE) expressed in absolute and a percentage of the coefficient of variation of the repeated measure (α: angle with respect to x-axis, β: angle with respect to y-axis, γ: angle with respect to z-axis).  Table 2. Reproducibility of the rat and rabbit joint contact area under virtual loading in terms of intraclass correlation coefficient (ICC) and precision errors (PE) expressed in absolute and a percentage of the coefficient of variation of the repeated measure ( χ : distance travelled to first contact, m : rate at which contact area increases).

Discussions
In this study, an image processing pipeline, that automatically prepares 3D joints images for sensitive and reproducible joint QMA measurements, does not require adaptation between animal models and is computationally efficient, was developed. Using SPHARM spherical harmonics modelling, the pipeline describes the tibia's basic form and determines the transformation needed to align the joint. The second module subdivides the joint into medial and lateral compartments by locating the tibia's intercondylar eminence as the watershed point using persistence homology to divide the VOI. The novelty of this work lies in the development of a new image processing pipeline that has been proven to prepare images efficiently and automatically for high-quality QMA, while remaining robust across two different preclinical animal models. This allows researchers to effectively perform 3D joint QMA in studies involving a larger number of samples, such as in longitudinal studies, with less expertise and a reduced time requirements compared to manual processing. Coupled with its robustness across rat and rabbit species, two common preclinical models for OA research 32 , the pipeline has the potential to make 3D joint QMA a more accessible technique. The alignment and subdivision pipeline, used in conjunction with 3D joint QMA, could be applied to studies involving other musculoskeletal diseases beyond OA to reveal previously unknown quantitative changes to the joint. Moreover, although the proposed pipeline was developed and validated on microCT data, it should be noted that the workflow starts after segmentation. Consequently, if the segmented components of a joint are available with appropriate resolution (i.e., two segmented bone images for the joint centre of mass and joint space width measurements, and additional cartilage images for contact area under virtual loading), the framework could be applied to calculate 3D joint QMA for other imaging modalities (such as magnetic resonance imaging or CT), in addition to other joint sites.
The pipeline also utilises SPHARM in a novel application. From its original proposal 20 , SPHARM has been widely used in statistical shape modelling for analysis of many organs, while others have found applications in using SPHARM descriptors for efficient image rotation estimation to take advantage of its hierarchical representation property 37,38 . However, few works have taken advantage of SPHARM's ability to represent objects in Table 3. Reproducibility of the rat and rabbit joint space measurements in terms of intraclass correlation coefficient (ICC) and precision errors (PE) expressed in absolute and a percentage of the coefficient of variation of the repeated measure (JSW: mean joint space width, JSV: joint space volume, JSW.min: minimum joint space width, JSW.max: maximum joint space width). www.nature.com/scientificreports/ Table 4. 3D joint QMA results for rat and rabbit datasets obtained in previous studies through manual processing, and in this study using the proposed automatic workflow. Accuracy of the automatically processed rat and rabbit measurements are shown in terms of root-mean-squared error expressed in absolute (RMSE) and mean-normalised (NRMSE) form (JSW: mean joint space width, JSV: joint space volume, JSW.min: minimum joint space width, JSW.max: maximum joint space width, α: angle with respect to x-axis, β: angle with respect to y-axis, γ: angle with respect to z-axis, χ : distance travelled to first contact, m : rate at which contact area increases). www.nature.com/scientificreports/ hierarchical levels of details in the object space for image processing purposes. The approach used in this work eliminated complex computation needed to perform image registration with O(n 2 ) complexity 33 and focused instead on finding a single solution that aligned the principal components to the desired pose with O(n) complexity. This performance improvement can be highlighted in Table 5, where the average time needed to calculate the transformation matrix for both rat and rabbit dataset were approximately 27 s, while the B-spline image transformation performed in the subsequent joint alignment step took much longer to complete. Moreover, the pipeline's use of an efficient algorithm for subdivision-locating persistence homology pairs by Kozlov and Weinkauf 25 , which has O(nlogn) complexity-allows for highly efficient processing overall. On average, the pipeline finished processing each sample in 693 s for rat and 2,112 s for rabbit, which is a great reduction from the 2,700 s for each rat and 5,400 s for each rabbit needed in earlier studies 11,12 with most of the reduction coming from automation of the alignment. As shown in Table 5, the major difference in processing time between both animals can be traced back to the SPHARM modelling and subsequent joint alignment processes, where the average time needed for a rabbit sample is significantly higher than that of a rat. This is likely due to the differences in animal size and imaging scale (imaged volume and resolution) as described in the imaging protocol and resulted in rabbit images containing significantly more voxels than rat images (5.52 × 10 9 more voxels on average). Therefore, more time was needed to perform SPHARM modelling and transformation. Other processes did not rely on the resolution of the images and resulted in similar processing times across animal models.
With regards to the quality of the QMA measurements, the pipeline is able to produce QMA measurement results with excellent reproducibility (ICC > 0.75) and precision errors (< 2% PE(%CV)) for the centre of mass and < 10% PE(%CV) for contact area but have comparable or slightly lower ICC compared to the gold-standard of manual processing by experts. As seen in Table 1, reproducibility values for the joint centre of mass from the pipeline are, generally, very slightly lower than that obtained from manually processed images from previous studies 11,12 . Reproducibility of contact area measurements, however, showed a mix of superior and inferior values when comparing with previous works. In general, manual processing in the earlier rat study 12 produced slightly better ICC than those from the pipeline, while manual processing in the earlier rabbit study 11 showed more mixed results as seen in Table 2. This could point to the increased experience of the operators in the follow-up study in the rat model 12 as compared to the pioneering 3D joint QMA study of the rabbit model 11 .
It should be noted that the joint space measurements performed in the rat study 12 and the rabbit study 11 were not the same JSW modules. In the original QMA work 11 , JSW was measured as the distance of the joint space between the centre of the femoral condyle and the tibia; while, in the later work 12 , the JSW was directly measured in 3D using the SPECTRA consensus approach 26 with JSV also being evaluated. In this study, the SPECTRA consensus approach was used to calculate JSW and JSV on both rat and rabbit dataset. As with the centre of mass and contact area measurements, reproducibility of the rat JSW and JSV showed excellent (ICC > 0.9), though slightly lower, values compared to results obtained with manual preprocessing. However, the pipeline results for JSW.min M/L have significantly lower precision error PE(%CV) as seen in Table 3. For rabbit data, this study also presents novel 3D JSW and JSV measurements with corresponding reproducibility values in Table 3 which shows excellent reproducibility for all joint space measures, except JSW.max L with ICC of 0.608. This generally highquality result for previously unperformed joint space QMA for the rabbit dataset highlights the pipeline's ability to correctly subdivide the images into the appropriate VOI and can support new QMA parameters in the future. With regards to low reproducibility of JSW.max , this study's result affirms earlier study's recommendation 12 not to use in further analysis for rat knee. Moreover, the lower ICC for JSW.max L in rabbit sample also lead this study to recommend JSW.max not be used for further analysis in rabbit knee as well.
The 3D joint QMA values measured using the proposed automatic workflow were shown to have similar values to those from earlier studies in Table 4. For both datasets, all measurements have errors of up to 9.5%, highlighting the sensitivity of 3D joint QMA measures to acquisition alignment. However, it should be further noted that there are no standard values of acceptable precision (PE(SD) and PE(%CV) and accuracy (RMSE and NRMSE) as they are highly dependent on the measurement context. Values of error for these measurements should be considered acceptable when the measurement approach provides results which are sensitive to the effects of disease or treatment, while maintaining consistent reproducibility values.
Since SPHARM modelling was used as the basis for its alignment algorithm, one of the essential requirements is that the input tibia must be spherical in nature. This presents one of the main limitations of the pipeline as any samples such as those bone with convex surfaces at the edge of the image due to the VOI cropping through Table 5. CPU time for each process of the pipeline expressed in mean (± SD) seconds. Total time for manual processing performed in earlier studies is shown on the last column for comparison.

Alignment module Subdivision module Framework
Manual processing from earlier studies 11 www.nature.com/scientificreports/ the medullary canal will not meet this requirement and will fail during the SPHARM processing step. Some manual closing, and smoothing operations will have to be performed on these samples prior to being processed through the pipeline. Additionally, even though the pipeline has solved the issue of joint alignment in the image processing stage, it is clear that systematic positioning of the joint during the acquisition of the images is fundamental to ensure that all joints are captured in a similar pose. In joint imaging where multiple rigid structures are of interest, it is not possible to simultaneously align both the femur and tibia without altering their original relative position which provides important pathological information 34 . Similar challenges have been noted in clinical measurement of radiological JSW, where reproducible patient position is required for reliable tracking of joint changes 35,36 . Future improvements in this direction would be in the form of a standardised positioning for acquisition that would allow the joints to be scanned and rescanned with minimal difference in pose between each. Further investigation, in strict quantitative terms, the definition of a gold standard alignment from which joint QMA can be precisely measured should also be done. Together, it is expected that even higher measurement reproducibility could be achieved.
Moreover, further in vivo experiments with regards to sensitivity and reproducibility of the proposed workflow and the novel 3D joint QMA parameters is required to determine whether these ex vivo results can be replicated in vivo. Movement artifact in live animal may reduce the precision and accuracy of these measurements, requiring further development of the workflow.

Conclusion
In previous work, quantitative measurements of joint morphometry using QMA has shown potential as a platform to quantify disease-based morphometric features for joint research using multiple preclinical animal models. However, its accessibility and usage are limited by high sensitivity to alignment and joint subdivision which are technically challenging and time-consuming to implement manually. In this work, we developed an automatic, efficient, and model-invariant image processing pipeline. The software was found to allow 3D joint QMA measurements with excellent reproducibility comparable to those obtained from manual processing in earlier studies.

Data availability
Datasets and all source code generated, used, and/or analysed during the current study are available from the corresponding author on reasonable request.