An Analysis by Synthesis Approach for Automatic Vertebral Shape Identification in Clinical QCT

TL;DR

This paper introduces a model-based method for automatic cortical surface identification in clinical QCT scans, improving accuracy in osteoporosis assessment despite limited spatial resolution.

Contribution

It develops a statistical bone model and an optimization approach to accurately locate the vertebral cortex in QCT images, demonstrating effectiveness both ex-vivo and in-vivo.

Findings

Accurately identified cortical surface ex-vivo using phantom and cadaveric data.

Demonstrated in-vivo applicability with manual surface comparison.

Improved cortical detection despite QCT's limited resolution.

Abstract

Quantitative computed tomography (QCT) is a widely used tool for osteoporosis diagnosis and monitoring. The assessment of cortical markers like cortical bone mineral density (BMD) and thickness is a demanding task, mainly because of the limited spatial resolution of QCT. We propose a direct model based method to automatically identify the surface through the center of the cortex of human vertebra. We develop a statistical bone model and analyze its probability distribution after the imaging process. Using an as-rigid-as-possible deformation we find the cortical surface that maximizes the likelihood of our model given the input volume. Using the European Spine Phantom (ESP) and a high resolution \mu CT scan of a cadaveric vertebra, we show that the proposed method is able to accurately identify the real center of cortex ex-vivo. To demonstrate the in-vivo applicability of our method we use manually obtained surfaces for comparison.

Tables

Table 1. Table 1 : Results of the phantom experiment. Given are for each of the three ESP vertebrae the radii and heights of the estimated surface. For both variables the mean, the standard deviation (SD) and the difference to the ground truth values (Diff) are shown. All measures are in mm. N denotes the number of samples on each surface (vertical cortex for radius, endplates for height).

	Radius					Height
Vertebra	N	Mean	Diff	SD		N	Mean	Diff	SD
Low	10072	17.73	-0.02	0.14		6101	23.95	-0.05	0.28
Medium	9897	17.52	0.02	0.16		6239	23.64	-0.36	0.42
High	9602	17.23	-0.02	0.16		23306	23.00	0.00	0.30

Equations

$$Y(t)=Y_0+(Y_1-Y_0)\cdot H(t+W)+(Y_2-Y_1)\cdot H(t-W),$$

$$g_{O,h}(z):={\tilde{g}}_O\left(\frac{z}{h}\right),$$

$$g_I(x,y):=\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{x^2+y^2}{2{\sigma }^2}\right).$$

$$\tilde{Z}(r,z)$$

$$Z(t)$$

$${\mu }_{\zeta }(t|w)$$

$${\sigma }_{\zeta }^2(t|w)$$

$$f_{\zeta }(z,t|w)=N\left(z|{\mu }_{\zeta }(t|w),{\sigma }_{\zeta }^2(t|w)\right).$$

$$f_{Z,W}(z,w,t)$$

$$f_W(w)$$

$$f_Z(z,t)=\int f_{\zeta }(z,t|w)f_W(w)\mathrm{d}w.$$

$$E(\hat{\mathcal{S}})=\underset{R_i\in SO(3)}{\min }\left\{{\sigma }_E^{-2}\sum _{i=1}^N{\gamma }_i\sum _{j\in \mathcal{N}(i)}{\omega }_{ij}\Vert ({\hat{\vec{x}}}_i-{\hat{\vec{x}}}_j)-{\vec{R}}_i({\vec{x}}_i-{\vec{x}}_j){\Vert }_2^2\right\},$$

$${\mathcal{L}}_Z({\hat{\vec{x}}}_i|{\hat{\vec{\rho }}}_i)=\sum _{j=1}^{2K+1}\log f_Z({\hat{\rho }}_{ij},t_j).$$

$${\mathcal{L}}_{E|Z}({\hat{\mathcal{S}}}_i|{\hat{\vec{\rho }}}_i)\propto {\mathcal{L}}_Z({\hat{\vec{x}}}_i|{\hat{\vec{\rho }}}_i)-E(\hat{\mathcal{S}})-\log \int f_{Z|E}({\hat{\vec{\rho }}}_i|\hat{\mathcal{S}})\mathrm{d}\hat{\mathcal{S}}.$$

$${\mathcal{L}}_{Z|s_i}({\vec{x}}_i|{\vec{\rho }}_i,s_i)={\mathcal{L}}_Z({\vec{x}}_i+s_i{\vec{n}}_i|{\vec{\rho }}_i).$$

$${\mathcal{L}}_{s_i|Z}(s_i|{\vec{x}}_i,{\vec{\rho }}_i)=$$

$${\mathcal{L}}_{Z|s_i}({\vec{x}}_i|{\vec{\rho }}_i,s_i)+{\mathcal{L}}_{s_i}(s_i)-\log \int f_{Z|s_i}({\vec{\rho }}_i|{\vec{t}}_i+\tau )f_{s_i}(\tau )\mathrm{d}\tau .$$

$$E_{\text{shape}}(\hat{\mathcal{S}})=\sum _{i=1}^N{\gamma }_i{\left(({\hat{\vec{x}}}_i-{\vec{y}}_i{)}^T\overrightarrow{n_i}\right)}^2+E(\hat{\mathcal{S}}),$$

$$\left(2{\sigma }_E^{-2}L\otimes I_3+B\right)\vec{x}=2{\sigma }_E^{-2}\vec{c}+\vec{d},$$

Full text

1 Introduction

Osteoporosis is a systematic skeletal disease that is characterized by low bone mass and deterioration of bone microstructure resulting in high fracture risk [6]. Its high prevalence of 24% in women beyond age of 65 makes osteoporosis a wide spread disease and a highly relevant research topic [7]. Quantitative computed tomography (QCT) has become a reliable tool for osteoporotic fracture risk prediction [1, 11] and monitoring [10, 8]. Volumetric trabecular bone mineral density (BMD) was identified as a good marker for bone strength. Still, the high under-diagnosis rate of 84% [28] indices that trabecular BMD alone is not sufficient as a bone strength marker. In osteoporotic patients the cortex takes the main load [25]. Therefore, the vertebral cortical bone is a worthwhile subject of study [12, 13]. However, the assessment of cortical markers like cortical thickness or cortical BMD is a challanging task [22]: the spatial resolution of clinical QCT ranges from 0.3-0.5 mm in-plane and from 1-3 mm out-of-plane. The thickness of the cortex of a vertebral body is reported [27, 24] to range from 0.25 mm to 0.4 mm and is therefore clearly below the Nyquist-Frequency, resulting in tremendous overestimation of cortical thickness in clinical QCT. Using high resolution QCT (HR-QCT) an in-plane resolution of up to 0.15 mm and an out-of-plane resolution of about 0.3 mm can be achieved at the expense of higher radiation dose, but cortical thickness is still clearly overestimated [10].

In this paper we address the problem of identifying the center of the cortex of a vertebral body from clinical QCT scans without any user interaction. As can be seen in figure 1 the apparent ridge of the cortical bone, i.e. the surface of maximum intensity, moves when the ratio of cortical to trabecular BMD changes. The same is true for different scanner resolutions [22] and cortical thicknesses. The strength of this effect does also vary with scanner, resolution and reconstruction kernel, making results from different scanners hard to compare [9]. A direct deconvolution of the resulting image is not applicable due to low signal to noise ratio.

1.0.1 Related Work

Shape identification and segmentation of vertebral bodies has been a research topic for a long time. Kang et al. [17] use a region growing approach followed by a refinement step based on relative thresholding to distinguish soft tissue from bone. While the method yields accurate results for thick cortices (5 times the voxel size and above), it shows the typical overestimation for thin cortices as present in human vertebrae. Mastmayer et al. [19, 18] proposed a multi step semi-automatic segmentation method based on the Euler-Lagrange equation and local adaptive volume growing that yields promising results, but again suffers from partial volume effects for thin cortices. The graph cut method proposed by Aslan et al. [3, 2] shows superior performance over previous work, but since the segmentation is voxel based, voxel size, especially out-of-plane, remains the limiting factor.

All of the methods mentioned yield voxel masks as the segmentation result. Therefore, the voxel size constraints the accuracy of cortex identification. Treece et al. [30] developed a mathematical model of the bone anatomy and applied a simplified model of the imaging system. The resulting measurement model is then fitted to the data given an initial segmentation of the target bone. This way they were able to provide an unbiased cortical thickness estimate down to $0.3$ mm on proximal femur.

1.0.2 Our Approach

We propose an analysis by synthesis (AbS) based approach (also called direct method) to accurately fit a template surface to the center of the cortical bone underlying a clinical QCT scan. The idea behind AbS based image analysis is not to analyze image features directly but instead synthesize an artificial image of a parametrizable model and find the parameters for which the synthetic and the input image match best. This way no image derivatives are required, making the process robust to noise. This model based approach can, to a certain degree, compensate for loss of information in the imaging process by incorporating prior knowledge into the model. However, for CT a full synthesis would require a radon transform of the full model followed by a CT reconstruction. Since this full synthesis is very expensive computationally, we simplify the scanning process to a blur, implemented by a convolution with the point spread functions (PSF) of the system. To simplify the synthesis further, we make use of sparse synthesis [14, 23] where not the whole image is synthesized, but only a sparse subset of it. Here, the result of the sparse synthesis is a set of one dimensional profiles orthogonal to the cortical surface, equivalent to a sparse sampling of a full synthesis111The equivalence is actually not given for the full volume. In areas where two cortices are close together, the sparse synthesis differs from the full synthesis. However, these regions are not critical for clinical routine.. Our model consists of a closed genus 0 surface representing the center of the cortex. As in medial representations [26], we assign a thickness and a BMD value to every point on the surface. The trabecular region inside the bone and the soft tissue region outside are represented by BMD distributions. Using a maximum a posteriori (MAP) estimation, we find the surface that maximizes the likelihood of the model parameters given measurements of the input volume.

In chapter 2 we give a detailed insight into the proposed method. We develop a statistical model of bone by extending the deterministic model of Treece et al. [30] in chapter 2.1. The sparse synthesis yields a statistical measurement model which we derive in chapter 2.2. In chapter 2.3 we show how the MAP estimation can be carried out in a data parallel process by employing the as-rigid-as-possible (ARAP [29]) deformation scheme. In chapter 3, we evaluate the accuracy of our method by using the European Spine Phantom (ESP, [16]). In addition, we further evaluate the accuracy using a µCT scan of a cadaveric human vertebra. To show the applicability of our method to clinical data we compare our results on 100 in-vivo QCT scans with manually obtained annotations. Finally, in chapter 4 we conclude this article and discuss future work.

2 Method

The input to our method is a calibrated222Houndsfield units are converted to bone mineral equivalents using known densities of a calibration phantom which is simultaneously scanned with the patient. QCT scan of a vertebra, a pre-estimated statistical measurement model and a labeled three-dimensional sketch of the target bone as depicted in figure 2(a). Each label of the template corresponds to a differently parameterized statistical bone model. The output is a triangle mesh representing the unbiased surface of the cortex center.

2.1 Statistical Bone Model

We model the surface of the cortex center as a closed genus 0 triangle mesh $\mathcal{S}$ with $N$ vertices. Its piecewise linear embedding is given by the vertex positions $V\subseteq\mathbb{R}^{3}$. When we look at a one dimensional profile orthogonal to the surface at any point on $\mathcal{S}$ (figure 2(b)), then the BMD graph is a piecewise constant function of the signed distance $t$ to the cortex center. By modeling the three density levels (soft tissue density, cortical BMD, trabecular BMD) as gaussian random variables $Y_{i}\sim\mathcal{N}(\mu_{Y_{i}},\sigma_{Y_{i}}^{2})\,,i=0,1,2$ and by modeling the cortical thickness333The cortical thickness is $2\cdot W$. as a random variable $W$ with444$W$ is a non-negative size quantity which is commonly modeled as a log-normal distribution. $\log W\sim\mathcal{N}(\mu_{W},\sigma_{W}^{2})$, we obtain a stochastic process $Y(t)$ of the profile:

[TABLE]

where $H$ denotes the Heaviside step function. Note that (1) takes a similar form as [30, eq. 1] with all unknowns replaced by random variables.

Since the density distribution of the soft tissue is not equal everywhere, for example the intervertebral discs have higher density than muscles, we use differently parameterized random processes for different sections of the bone. Figure 2(a) depicts the three different regions used in our model: vertical cortex (green), endplates (red), foramen (yellow). The ”cut pedicles” region (blue in fig. 2(a)) is not used in the model.

2.2 Measurement Model

The next step in the AbS framework is the synthesis step where the imaging system is simulated by applying the in-plane and out-of-plane PSF to the bone model. By sampling the resulting image one can then acquire observable profiles. However, with sparse synthesis we are able to simplify the process by not performing the convolutions in the global coordinate system, but in the local system of the profile.

2.2.1 Combined PSF

The slice sensitivity profile (SSP) of a spiral CT can be approximated by the convolution of the rectangular profile of helical CT by the triangular table movement function [15]. The width of the rectangular profile is determined by the collimation, while the width of the movement function is determined by the table feed. We therefore approximate the out-of-plane PSF (perpendicular to a CT slice) for pitch factor of 1 with

[TABLE]

where $\tilde{g}_{O}(z)\coloneqq\mathrm{\Pi}(z)*\mathrm{\Lambda}(2z)$ and $h$ is the slice width. Like [30] we approximate the in-plane PSF with a rotational invariant gaussian of width $\sigma$:

[TABLE]

If we slice through the cortex at $\vec{x}\in\mathcal{S}$ with a plane defined by the z-axis and the surface normal $\vec{n}$ at $\vec{x}$, we can define a local coordinate system centered on $\vec{x}$ with the z-axis as the ordinate and the projection of $\vec{n}$ onto the x-y-plane as the abscissa (r-axis) as depicted in figure 3. Since the profile now lies inside r-z-plane and because $g_{I}$ is rotational invariant, we can obtain the measurement process $\tilde{Z}(r,z)$ by convolving $\tilde{Y}(r,z)\coloneqq Y(r\sin\theta+z\cos\theta)$ with $g_{O,h}$ along the z-axis and with $g_{I}$ along the r-axis:

[TABLE]

where $t\coloneqq r\sin\theta+z\cos\theta$ and $\theta$ is the angle between $\vec{n}$ and the z-axis. By substituting $\phi\coloneqq\tau\sin\theta$ and $\psi\coloneqq\lambda\cos\theta$ we get

[TABLE]

Hence, we can simplify the two dimensional convolution by a single convolution of the one-dimensional process from (1) with an angle dependent PSF $g_{\theta}$.

Although, we used a gaussian for the in-plane PSF in the derivation of the combined PSF, we note that any symmetric square integrable PSF can be used here. Ohkubo et al. [21] determined the PSFs for several reconstruction kernels. Based on their measurements one can observe that a gaussian is a good approximation for smooth kernels like the Siemens B40, but for sharper kernels like the Siemens B80 it is not: the B80 amplifies some higher frequencies in a narrow band to enhance edges while the gaussian damps all high frequencies.

2.2.2 Stochastic Measurement Process

Let $G_{\theta}\coloneqq\int g_{\theta}(t)\mathop{}\!\mathrm{d}t$ be the primitive function of the angle dependent PSF $g_{\theta}$ from (2.2.1). Since $g_{\theta}$ has finite energy, we know that $\displaystyle\lim_{t\rightarrow-\infty}G_{\theta}(t)=0$ and $\displaystyle\lim_{t\rightarrow+\infty}G_{\theta}(t)=1$. With this we can define the optimal (i.e. noise free) stochastic measurement process $Z(t)$ as

[TABLE]

Since $Y_{0},Y_{1},Y_{2}$ are stochastically independent normally distributed random variables the conditional stochastic process $\zeta(t)\coloneqq Z(t)|W$ is gaussian with mean and variance function given as:

[TABLE]

Thus the probability density function (PDF) $f_{\zeta}$ of $\zeta(t)$ is given by

[TABLE]

With the joint probability density function of $Z(t)$ and $W$

[TABLE]

the probability density function $f_{Z}$ of $Z(t)$ is obtained by marginalizing out $W$:

[TABLE]

2.3 Deformation Model and Optimization

Given a profile through a point $\vec{x}\in\mathcal{S}$ of an input volume $I$, we can express the degree of conformity of the profile with a synthetic profile by means of the likelihood of the measurement model. The optimal surface $\mathcal{S}^{\prime}$ is then defined by the linear embedding $V^{\prime}\subset\mathbb{R}^{3}$ for which the likelihood of all profiles through $\mathcal{S}^{\prime}$ are maximal. However, this is an ill posed problem, because it is not guaranteed that the maximizer yields a realistic shape. By constraining the deformation between an initial surface $\mathcal{S}$ and $\mathcal{S}^{\prime}$ to be as rigid as possible, shape degeneration can be avoided and the problem becomes well formed.

We modify the as-rigid-as-possible (ARAP) energy term from [29] to act as a shape prior in the MAP framework as follows:

[TABLE]

where $\omega_{ij}$ are the cotangent weights of edge $(i,j)$ [20], $\gamma_{i}$ are per-vertex weights and $\mathcal{N}(i)$ contains all vertices of the 1-ring of vertex $i$. $\vec{R}_{i}\in SO(3)$ are per-vertex rotation matrices adding additional $3N$ degrees of freedom to the deformation. $E(\hat{\mathcal{S}})$ can be seen as a measure of how far the deformation is from an isometry [4]. With the scale $\sigma_{E}^{-2}$ equation (11) is equivalent (up to constants) to the negative log likelihood function of a gaussian random variable $E$ modeling the ARAP energy of $\mathcal{S}^{\prime}$.

Now, let $\vec{t}=(t_{i})_{1,\dots,2K+1}$ be a vector of $2K+1$ samples from $[-t_{0},t_{0}]\subset\mathbb{R}$ with a sampling interval $\delta_{t}$. For all $i=1,\dots,N$ let $\hat{\vec{l}}_{i}\coloneqq\{\hat{\vec{x_{i}}}+t_{j}\cdot\hat{\vec{n}}_{i}\,|\,j=1,\dots,2K+1\}$ be a vector of discrete samples of a profile through position $\hat{\vec{x}}_{i}\in\hat{V}$ and let $\hat{\vec{\rho}}_{i}\coloneqq I(\hat{\vec{l}}_{i})$ be the densities obtained by sampling the interpolated555In our implementation we use trilinear interpolation for performance reasons. input volume $I$ at positions $\hat{\vec{l}}_{i}$. The log likelihood $\mathcal{L}(\hat{\vec{x}}_{i}|\hat{\vec{\rho}}_{i})$ of $\hat{\vec{x}}_{i}$ given measurements $\hat{\vec{\rho}}_{i}$ is given with (10) by666To keep the derivation simple, we assume independence here. In general the samples are correlated by the interpolation method.

[TABLE]

The posterior log likelihood is then proportional to

[TABLE]

The direct optimization of (13) is analytically intractable and numerical optimization would either require the computation of image gradients, which we are tying to avoid, or gradient less optimization methods which do not perform very well on high dimensional optimization problems. We therefore split the optimization into two sub problems that can be efficiently solved: first the optimal displacement along the surface normal for each profile is estimated. Afterwards, the surface is optimized to match the displaced positions under ARAP constraints. This two step optimization scheme is iterated until convergence.

2.3.1 Optimal Displacements

We first modify equation (12) by adding a latent variable $s_{i}\sim\mathcal{N}(0,\sigma_{s_{i}}^{2})$ as follows:

[TABLE]

Using MAP estimation we can find the optimal displacement $\hat{s}_{i}$ by maximizing the posterior likelihood

[TABLE]

In our implementation we pre-evaluate $f_{Z}$ for discrete samples of $t,\theta$ and $z$ to a three dimensional histogram and approximate $f_{Z}$ by histogram lookup with trilinear interpolation. This way, equation (15) can be maximized very efficiently using data parallel exhaustive search along $s_{i}$ at discrete steps.

2.3.2 Surface Fitting

Now that the optimal displacements are known, the surface needs to be fitted to the displaced positions $\vec{y}_{i}\coloneqq\vec{x}_{i}+\hat{s}_{i}\vec{n}_{i}$ under ARAP constraints. We therefore minimize the following energy term:

[TABLE]

where the first term addresses the point to plane distances and the second term the ARAP constraints. There might be profiles $\vec{\rho}_{i}$ where no cortex can be observed. Since the estimated displacements from those profiles are meaningless, we use the per-vertex weights $\gamma_{i}$ to down-weight those estimates by using the posterior PDF: $\gamma_{i}\coloneqq f_{s_{i}|Z}(\hat{s}_{i}|\vec{x}_{i},\vec{\rho}_{i})$.

Like in [29] we use an alternating iterative optimization scheme to optimize (16): we first keep the rotations $\vec{R}_{i}$ fixed and optimize for the positions $\hat{\vec{x}}_{i}$ and then optimize for the rotations while keeping the positions fixed. For fixed positions, the optimal rotations can be found by SVD (refer to [29, eqs. 5,6] for details). The optimal positions $\hat{\vec{x}}_{i}$ can be found by setting the partial derivatives of $E_{\text{shape}}$ to zero, which results in the following sparse linear system of equations:

[TABLE]

where $\vec{x}=\left(\vec{x}_{1}^{T},\dots,\vec{x}_{N}^{T}\right)^{T}$, $L$ is the Laplacian matrix of $\mathcal{S}$, $B$ is a block diagonal matrix with entries $(\gamma_{i}\vec{n}_{i}\vec{n}^{T})_{ii}$ and $\otimes$ denotes the Kronecker product. The vector $\vec{c}$ contains the ARAP constraints as a concatenations of vectors $\vec{c}_{i}=\sum_{j}L_{ij}\frac{(\vec{R}_{i}+\vec{R}_{j})}{2}(\vec{x}_{i}-\vec{x}_{j})$ and $\vec{d}$ consist of concatenations of vectors $\vec{d}_{i}=\gamma_{i}\vec{n}_{i}\vec{n}_{i}^{T}(\vec{x}_{i}+\hat{s}_{i}\vec{n}_{i})$. Equation (17) can be efficiently solved using a preconditioned conjugate gradient solver (PCG). The alternating optimization scheme is iterated until convergence.

3 Experiments and Results

We implemented the proposed method using MATLAB and C++. The measurement model requires the slice spacing $h$ and the width $\sigma$ of the in-plane PSF as input parameters. The value of $\sigma$ can be easily estimated from phantom scans777Since the calibration phantom is present in all QCT scans, no separate scanning process is required.. The parameters of the BMD priors $Y_{i}$, can be estimated from single measurements or prior knowledge. The parameters of the width prior $W$ are set to reflect the range of reported cortical thicknesses. For each scanner configuration $(h,\sigma)$ we pre-computed $f_{Z}$ for discrete samples from $(t,\theta,z)\in[-2,2]\times[0^{\circ},90^{\circ}]\times[-1000,2000]$ (41, 91 and 3001 samples, respectively) data parallel on a GPGPU and saved the result as a 3D histogram as described in section 2.3.1. We set $\sigma_{s_{i}}=2$ and $\sigma_{E}=2$ for all experiments. The optimization usually provides good results after a few iterations and takes about $1$ to $5$ minutes to complete on our system (Intel® Core™ i7-4790 CPU @ 3.60GHz, 4 Cores), depending on the image size and number of iterations888In our prototype implementation, the optimization does not utilize the GPU, yet, but we note that all operations can be easily ported to the GPU..

3.1 Ex-Vivo

To evaluate the accuracy of the proposed method, we scanned the ESP with a clinical CT scanner using a low dose protocol (Siemens Somatom 64, 120 kV, 80 mAs, kernel B40s) at an in-plane resolution of 0.4 mm and a slice spacing of 1 mm. The ESP consists of three geometrical phantom vertebra (low, medium and high), each with different wall thicknesses and densities. The proposed method was used to acquire the cortical surfaces $\mathcal{S}_{i}$ of the three phantom vertebrae. To gain a dense surface, we used random mesh sampling of $\mathcal{S}_{i}$ using the method described in [5]. For the accuracy evaluation of the vertical cortex, we first fit a cylinder to the point cloud belonging to the vertical cortex and then computed the radius for each point separately. Since the vertebral bodies of the ESP have a diameter of 36 mm, the optimal radii of the cylindrical surfaces through the cortex centers are $17.75$ mm ($0.5$ mm wall), $17.5$ mm ($1$ mm wall) or $17.25$ mm ($1.5$ mm wall), for the low, medium and high vertebra, respectively. To evaluate the accuracy of the endplate surface we fit two planes to our point cloud: one for the upper endplate and one for the lower endplate. Hereafter, the point to plane distance of each sample to its opposite plane was computed. All bodies of the phantom have a height of 25 mm, so the optimal distances are 24 mm for the low and medium vertebra ($1$ mm wall) and 23 mm ($2$ mm wall) for the high vertebra.

Table 1 summarizes the results of the ESP experiment. For all vertebra levels, the difference between the mean radius of the estimated surface and the ground truth radius is near zero. The standard deviation is below one half of the in-plane resolution. The same applies for the estimated heights, noting that the out-of-plane resolution is $2.5$-times lower than the in-plane resolution.

Since the shape of the ESP is very simple, we used a cadaveric vertebra embedded in resin for a more realistic reference. The embedded vertebra was scanned with a µCT system (SCANCO Medical, 70 kV, 360 mAs) with an isotropic resolution of 31 µm and with a clinical QCT system (Siemens Somatom 64, 120 kV, 100 mAs, kernel B40s) with an in-plane resolution of 0.2 mm and a slice spacing of 1 mm. Both scans where calibrated and the QCT scan was resampled and rigidly registered to the µCT scan. The resulting rigid transformation matrix was saved for later use. We applied our method to the unregistered QCT scan of the embedded vertebra and sampled the resulting surface as above. Using the inverse transformation matrix the point cloud was transformed into the coordinate system of the µCT scan. To evaluate the accuracy of the cortical surface we first need to identify the ground-truth cortex in the µCT scan. We sampled the µCT scan along lines orthogonal to the acquired surface at every sample point, from 5 mm outside to 5 mm inside the volume defined by the surface using $2001$ samples per line ($5$ µm spacing). We binarized the profile using a threshold of $500$ mg/cc. To fill small cavities in the cortex (see figure 5), we applied a morphological closing operation to the binarized signals. We defined the periosteal (outer) surface as the first rising edge and the endosteal (inner) surface as the first falling edge. The center of the cortex is then the midpoint between the periosteal and the endosteal surface. Since the surface estimated by the proposed method is located at the center of the sampled line ($t=0$), the signed distance to the real cortex center is simply the location of the midpoint.

We were able to identify the cortex center for $8506$ samples for the vertical cortex and for $8273$ samples for the endplates using this method. The average distance to the center of the vertical cortex was $0.0662\pm 0.2327$ mm (mean $\pm$ standard deviation) and $0.0607\pm 0.2347$ mm for the endplates. The absolute error was below $0.075$ mm for $25$%, below $0.16$ mm for $50$% and below $0.28$ mm for $75$% of the samples. Figure 4 depicts the spatial distribution of the absolute errors on the surface. The highest errors are scattered around the surface, but do not form larger clusters. We manually inspected the profiles with the highest errors and found that most of them are located at positions where the real cortex is hard to identify by the thresholding method. Figure 5 shows an example of such a location. While the cortex centers identified by the thresholding method vary quite widely, the estimated surface stays smooth between the periosteal and the endosteal surfaces.

3.2 In-Vivo

For the in-vivo evaluation of the proposed method we used 100 clinical QCT scans (Phillips Brilliance 16, 120 kV, 90-420 mAs) of the first lumbar vertebra from an osteoporosis study with an in-plane resolution of 0.7 mm and a slice thickness of 0.8 mm (increment 0.4 mm). The scans were calibrated and segmented by an expert operator. In the segmentation process, a triangle mesh representing the apparent cortical ridge surface was exported. We used those meshes as a reference to evaluate the applicability of the proposed method to in-vivo data. However, since the manually obtained meshes do not identify the center of the cortex but the apparent cortical ridge, a bias is unavoidable. After applying the proposed method to the in-vivio scans, we did the same random sampling as in the ex-vivo case and afterwards computed for each sample the distance to the reference mesh using [5].

For the vertical cortex we acquired $98000$ samples ($980$ samples per patient) with an average distance of $0.25\pm 0.52$ mm (mean $\pm$ standard deviation). For the endplates we acquired $101600$ samples ($1016$ samples per patient) with an average distance of $0.19\pm 0.56$ mm. By visual inspection of the spatial distribution of absolute errors, we found that the majority of high errors are located in the transition area between the vertical cortex and the endplates.

To check if the bias could be explained by the displacement between ridge and cortex center (cf. figure 1), we evaluated the 99% confidence intervals of $\zeta$ (8) for a cortical thickness of $0.3$ mm and for priors $Z_{0}$ and $Z_{2}$ estimated from the 100 in-vivo scans. For an angle of $\theta=90^{\circ}$, the displacement between the cortical ridge and the cortex center lies between $0.1$ mm and $0.3$ mm and for an angle of $\theta=0^{\circ}$ between $0.09$ mm and $2.3$ mm. This could explain the observed bias, but for the in-vivo case it cannot be directly verified. To substantiate our hypothesis, we did a statistical analysis based on measurements obtained using the original manual segmentation of the 100 in-vivo scans. There were no significant correlations between the bias and the integral BMD, cortical thickness, the total volume or the ratio of cortical to trabecular BMD. However, the cortical thickness measure that was used is based on voxel distances and therefore it is not very expressive. Therefore, we also investigated the density weighted cortical thickness (cortical thickness multiplied by cortical BMD divided by $1200$ mg/cc). We found a significant negative linear relationship between the bias and density weighted cortical thickness (wCt.Th). Figure 6 shows a scatter plot with the regression line. There is a significant negative slope of $-0.94$ ($p=0.019$). Therefore, if the cortex gets thicker, the bias gets lower. This is consistent with the cortex center shift as depicted in figure 1: a thicker cortex is equivalent to a higher resolution resulting in a better agreement between the cortical ridge and the cortex center.

4 Conclusion

We presented an analysis by synthesis approach for automatic vertebral shape identification in clinical QCT. The foundation of our method is a statistical model of bone which is convolved with the in-plane and out-of-plane PSF of spiral CT resulting in a statistical measurement model. We use an as-rigid-as-possible deformation scheme to find the surface of the center of cortical bone that maximizes the posterior likelihood of the model given the image. Since the posterior likelihood is analytically intractable, we propose an alternating optimization scheme to find an approximate solution using an efficient data parallel process.

The evaluation of the proposed method using a clinical QCT scan of the ESP and a µCT scan of an embedded vertebra demonstrates its high sub-voxel accuracy. The applicability to in-vivo data was shown by comparing the estimated surfaces to manual annotations of 100 in-vivo QCT scans from an osteoporosis study. The remaining bias to the manually obtained meshes might be explained by the displacement of the apparent cortical ridge to the cortex center due to low spatial resolution. We substantiated this hypothesis by a statistical correlation analysis.

We think the proposed method is a good starting point for further assessment of cortical bone markers. The estimation of the (local) cortical thickness and the cortical BMD should be possible by maximizing the posterior probability of the respective parameters, like we did for the displacements. We may also note that, although we presented our method here for vertebral shape identification, it is not limited to vertebrae. Using a different template mesh the method can also be applied to other bones or bone parts, e.g. proximal femur.

There are of course limitations of our method: it fails to accurately estimate the cortex center at the transition area between the vertical cortex and the endplates. In these areas the assumptions we make for our model are not fulfilled. If two cortices are close together, both cortices appear in a single profile, leading to a low likelihood for both cortex positions. However, for clinical assessment, those transition areas can be easily excluded. We also note that we did not evaluate the precision of our method, yet. Since the method is deterministic, it produces the same result for multiple runs on the same input image, but the determination of the precision using multiple re-located scans of the same object is still to be done.

To promote the comparability of QCT analysis algorithms we are making our software publicly available999https://github.com/ithron/CortidQCT under open source license.

4.0.1 Acknowledgments.

This work was part of the Diagnostik Bilanz Study which is part of the BioAsset project. BioAsset is funded by a grant of the Bundesministerium für Bildung und Forschung (BMBF), Germany, Förderkennzeichen 01EC1005. This work was also supported by the German Research Foundation, DFG, No. KO2044/9-1.