Morphometry-based measurements of the structural response to whole-brain radiation

Abstract

Purpose

Morphometry techniques were applied to quantify the normal tissue therapy response in patients receiving whole-brain radiation for intracranial malignancies.

Methods

Pre- and Post-irradiation magnetic resonance imaging (MRI) data sets were retrospectively analyzed in N = 15 patients. Volume changes with respect to pre-irradiation were quantitatively measured in the cerebrum and ventricles. Measurements were correlated with the time interval from irradiation. Criteria for inclusion included craniospinal irradiation, pre-irradiation MRI, at least one follow-up MRI, and no disease progression. The brain on each image was segmented to remove the skull and registered to the initial pre-treatment scan. Average volume changes were measured using morphometry analysis of the deformation Jacobian and direct template registration-based segmentation of brain structures.

Results

An average cerebral volume atrophy of $-$0.2 and $-$3 % was measured for the deformation morphometry and direct segmentation methods, respectively. An average ventricle volume dilation of 21 and 20 % was measured for the deformation morphometry and direct segmentation methods, respectively.

Conclusion

The presented study has developed an image processing pipeline for morphometric monitoring of brain tissue volume changes as a response to radiation therapy. Results indicate that quantitative morphometric monitoring is feasible and may provide additional information in assessing response.

Appendix: Gaussian approximation of the Green’s Kernel

Following a diffeomorphic image registration framework [4, 7, 14], consider two images, $I_0$ and $I_1$, defined on an Eulerian reference domain, $\varOmega \subset \mathbb {R}^3$. The goal of the image registration is to determine the motion, $\mathbf {\varphi }(\varvec{x},t): \varOmega \times [0,1] \rightarrow \varOmega$, that maps the reference image, $I_0$, to the current image, $I_1$.

$$\begin{aligned} I_0: \varOmega \rightarrow \mathbb {R} \qquad I_1: \varOmega \rightarrow \mathbb {R} \qquad I_1 = I_0 \circ \mathbf {\varphi }(\cdot ,1) \end{aligned}$$

A symmetric diffeomorphic registration is optimal with respect to a given image similarity metric, $d: \varOmega \times \varOmega \rightarrow \mathbb {R}$, and penalized by the velocity of the transformation.

$$\begin{aligned}&\inf _{\mathbf {\varphi }_1} \inf _{\mathbf {\varphi }_2} f(\mathbf {\varphi }_1,\mathbf {\varphi }_2) \quad f(\mathbf {\varphi }_1,\mathbf {\varphi }_2)\nonumber \\&\quad = \left( \begin{array}{ll} \int _{0}^{0.5} \left( \left\| \frac{d \mathbf {\varphi }_1}{dt}(t)\right\| ^2_L + \left\| \frac{d \mathbf {\varphi }_2}{dt}(t)\right\| ^2_L\right) \;dt \\ + \int _\varOmega d(I_0 \circ \mathbf {\varphi }_1^{-1}(\cdot ,.5),I_1 \circ \mathbf {\varphi }_2^{-1}(\cdot ,.5)) \end{array} \right) \nonumber \\&\qquad \mathbf {\varphi }= \mathbf {\varphi }_1 \circ \mathbf {\varphi }_2^{-1} \end{aligned}$$

(3)

Here, the deformations, $\mathbf {\varphi }_i, i=1,2$, are defined with respect to the midpoint, $t=.5$, of the transformation. Time is parameterized in opposite directions between $\mathbf {\varphi }_1$ and $\mathbf {\varphi }_2$. The operator norm $\Vert \cdot \Vert _L$ is induced by a differential operator of the type, $L=\alpha \varDelta + Id$, $\alpha \in \mathbb {R}$ [7]. The symmetric diffeomorphic formulation mappings are constructed sufficiently smooth such that the inverse of the motion is well defined, $\mathbf {\varphi }^{-1}(\varvec{x},t): \varOmega \times [0,1] \rightarrow \varOmega$, and gives a consistent solution for the forward and inverse mapping, $\mathbf {\varphi }\circ \mathbf {\varphi }^{-1} = Id$.

The Euler–Lagrange equations provide necessary conditions for which a solution of the optimization formulation (3) must satisfy.

$$\begin{aligned} \nabla _{v_i} f = L v_i + \nabla d (v_i) =0 \qquad i \;=\; 1,2 \end{aligned}$$

(4)

Here, the gradient of the objective function (3) is with respect to the velocity of the transformation. The operator $L$ represents a physics constraint on the deformation solution. Assuming the deformation behaves as a viscous fluid provides an intuitive solution field that may be understood to adhere to first principle conservation laws.

$$\begin{aligned} Lv = \mu \varDelta v + (\mu + \lambda ) \nabla (\nabla \cdot v) \end{aligned}$$

(5)

Further, compared with linear elastic displacement models that constrain the accuracy of large deformations because of internal elastic strain, accurate large deformations may be achieved within this viscous fluid model because internal forces disappear over time and the desired deformation can be fully achieved [8]. However, this approach leads to computationally expensive numerical solution schemes that couple the individual components of the deformation. Alternatively, assuming each deformation component is decoupled and diffuses along the respective gradient of the deformation field yields the algorithmically and numerically tractable Gaussian convolution kernel used in SyN [4]. The greedy update at the midpoint in time simplifies to a fixed point iteration on the velocity field [36], and the transformation field is iteratively updated through a finite difference approximation of the velocity

$$\begin{aligned} L&= - \varDelta + \frac{1}{\varDelta t} Id \nonumber \\&\Rightarrow - \varDelta v_i + \frac{1}{\varDelta t} v_i + \nabla d (v_i) =0\nonumber \\&\Leftrightarrow \left\{ \begin{array}{l} \varDelta v_i = \frac{\partial }{\partial t} v_i \approx \frac{v_i (\varDelta t) - v_i (0)}{\varDelta t}\\ v_i(0) = \varDelta t \; \nabla d (v_i) \end{array} \right. \nonumber \\&\Rightarrow K \star \varDelta t \; \nabla d (v_i) \approx v_i (\varDelta t) \approx \frac{\mathbf {\varphi }_i (\varDelta t) - \mathbf {\varphi }_i (0) }{\varDelta t} \end{aligned}$$

(6)

Here, $K \star$ represents the Gaussian convolution operation and is obtained from a Fourier transform solution of the heat equation. Solutions of the Euler–Lagrange equations may be interpreted as a decoupled component-wise solution to a isotropic heat transfer equation with initial conditions given by the gradient of the similarity metric.