1 Introduction

Craniosynostosis is a congenital skull malformation in which one or more cranial sutures fuse early, producing head malformations due to the compensatory growth of the brain along the non-fused sutures. These malformations can result in severe conditions such as increased intra-cranial pressure and impaired brain growth [1]. In case of single suture-fusion, depending on which suture is affected, craniosynostosis can be classified as metopic, coronal, sagittal, or lambdoid. The first fully automatic, objective and quantitative method for the diagnosis of craniosynostosis was presented in [2], where the authors used a multi-atlas of healthy cranial shapes to quantify deformities based on two different shape descriptors: (1) malformations, defined as the local Euclidean distances between the mesh representation of the patient’s cranium and its closest normal cranial shape from a normative multi-atlas, and (2) curvature discrepancies, defined as local curvature differences. However, how to correct cranial deformities was beyond the scope of that work.

The optimal repositioning and reconstruction of bone structures is a common problem in reconstructive surgery. Computer aided planning software for cranio-maxillofacial interventions [3] has improved significantly the capabilities of surgical planning. However, most of these tools need substantial manual interaction to create a plan for surgery, so the results are very dependent on the expertise of the specialists. In the specific case of cranial vault remodeling to treat metopic craniosynostosis, both the decision of performing surgery and the procedure itself are very subjective [2].

A fully automatic method to find the optimal post-surgical shape to target during cranial vault reconstruction surgery was presented in [4], where an image registration framework to calculate an interventional plan that minimizes the malformations in the patient’s cranium was presented. That framework had three shortcomings: (1) the frontal bones were considered as two different rigid objects, but cranial vault remodeling to correct metopic craniosynostosis very often involves bone bending/reshaping [5, 6]; (2) the algorithm, based on volumetric image registration, only took into account the Euclidean distance between the patient’s cranium and its closest normal shape, without considering smoothness and curvature (which is critical to obtain esthetically acceptable surgical results); (3) the method did not consider the division of the frontal bones into smaller bone pieces, which is typically needed to surgically correct cases with severe cranial deformities.

We propose an automated framework based on a novel locally affine surface-based diffeomorphic registration method to calculate the optimal surgical plan for the treatment of metopic craniosynostosis. The goal of the proposed work is to objectively define the bone cuts, bends and repositioning that will transform the patient’s cranial shape to best fit its closest normal shape from a multi-atlas. The proposed registration framework estimates one single transformation between a patient’s cranial shape and its closest normal shape, while it models individually the repositioning and bending of each local bone piece. This work overcomes the previously described shortcomings of [4] by: (1) allowing bone bending in addition to translation and rotation; (2) creating a cranial shape that is not only closer to normality in terms of point-to-point distances but also in terms of global shape and curvature; (3) allowing the surgeons to consider different interventional approaches (i.e. different bone cuts) by introducing bone cut templates to guide the subdivision of the cranial bones into smaller bone pieces.

2 Methods

2.1 Template-Based Bone Cut Labeling

We segmented the patient’s cranial bones from computed tomography (CT) images as described in [2]. In summary, the cranial bones were extracted from CT by thresholding the images based on the Hounsfield units of bone tissue. Then, a graph-cut based method was used to label each cranial bone. Based on [2], we created a reference template by manually segmenting the volumetric image of a healthy case and including a subdivision of the frontal bones, as shown in Fig. 1. Then, given the segmented cranial bones of a patient, the following cost function was minimized to create a labeling scheme that includes the subdivision of the patient’s frontal bones:

$$ E\left( f \right) = \sum\nolimits_{p\epsilon P} {\frac{{d_{{f_{p} }} }}{{d_{{f_{p} }} + \mathop {\hbox{max} }\nolimits_{f} (d_{f} |f \ne f_{p} )}}} , $$
(1)

where \( f \) is a labeling scheme assigning label \( f_{p} \) to voxel \( p \), and \( d_{f} \) is the distance to the cranial bone in the template volume with label \( f \). Figure 1(b) shows an example of frontal bone segmentation and bone piece labeling for one patient.

Fig. 1.
figure 1

(a) Anterior view of the cranium of one patient with the left (LF) and right (RF) frontal bones segmented. (b) Bone piece cuts labeled using the proposed bone cut template. Label 1 indicates the supra-orbital bar, labels 2–3 are subdivisions of the LF bone, and labels 4–5 are subdivisions of the RF bone. The rest of the cranium is shown in black.

Full size image

2.2 Poly-affine Transformation Model

The transformation model to create the optimal surgical plan must allow different bone pieces to be manipulated in an affine fashion during surgery. Let \( \varvec{T} = \left( {\varvec{L},\varvec{t}} \right) \) be an affine transformation with linear part \( \varvec{L} \) and translation \( \varvec{t} \). The transformation of a point at coordinates \( \varvec{x} \) expressed in homogeneous coordinates is calculated as \( \varvec{Tx} \). In [7], it was shown that \( \varvec{T} \) can be modeled as a temporal process (transforming the point from time \( s = 0 \) to time \( s = 1 \), \( s \in \left[ {0,1} \right] \)) in which we can associate a family of velocity vector fields to \( \varvec{T} \) by writing \( \varvec{v}\left( {\varvec{x},s} \right) = \varvec{v}\left( \varvec{x} \right) = \log \left( \varvec{T} \right)\varvec{x} \). The transformed coordinates of \( \varvec{x} \) at time \( s \) can then be calculated as

$$ \varphi \left( {\varvec{x},s} \right) = exp\left( {s\,log\left( \varvec{T} \right)} \right)\varvec{x}. $$
(2)

Given \( M \) affine transformations, they can be composed in a smooth and invertible way by combining the velocity associated to them using the following equation:

$$ \varvec{v}\left( \varvec{x} \right) = \frac{{\sum\nolimits_{i = 1}^{M} {w_{i} \left( \varvec{x} \right)\varvec{v}_{i} \left( \varvec{x} \right)} }}{{\sum\nolimits_{i = 1}^{M} {w_{i} \left( \varvec{x} \right)} }}, $$
(3)

where \( w_{i} \left( \varvec{x} \right) \) is the smooth and invertible function defining the contribution of the transformation \( i \) to the total velocity at coordinates \( \varvec{x} \).

As it was shown in [7], given the normalized weight functions \( w_{i} \left( \varvec{x} \right) \), the transformation \( \varphi \left( {\varvec{x},\varvec{ }\frac{1}{{2^{N} }}} \right) \) of a point at coordinates \( \varvec{x} \) to time \( s = \frac{1}{{2^{N} }} \) can be written as

$$ \varphi \left( {\varvec{x},\varvec{ }\frac{1}{{2^{N} }}} \right) = \left( {\sum\nolimits_{i = 1}^{M} {w_{i} \left( \varvec{x} \right)\varvec{T}_{i}^{{\frac{1}{{2^{N} }}}} } } \right)\varvec{x}, $$
(4)

where \( \varvec{T}_{i}^{{\frac{1}{{2^{N} }}}} \) is the \( 2^{N} \) root of the local affine transformation \( \varvec{T}_{i} \). Then, if the temporal interval is divided into \( 2^{N} \) intermediate points, the final transformation can be obtained by composing \( \varphi \left( {.,\varvec{ }\frac{1}{{2^{N} }}} \right) \) exactly \( 2^{N} \) times:

$$ \varphi_{n} \left( \varvec{x} \right) = \left( {\sum\nolimits_{i = 1}^{M} {w_{i} \left( {\varphi_{n - 1} \left( \varvec{x} \right)} \right)\varvec{T}_{i}^{{\frac{1}{{2^{N} }}}} } } \right)\varphi_{n - 1} \left( \varvec{x} \right), $$
(5)

where \( \varphi_{n} \left( \varvec{x} \right) \) is the transformation at intermediate point \( n \).

2.3 Locally Affine Regions

In [8], the weights associated to each local transformation were defined using Gaussian kernels centered at predefined anchor points. In that framework, no single region of the image was constrained to present the same affine behavior, but a weighted combination of them. In [4], the transformation at local bone regions was constrained to be strictly rigid. To that end, signed distance functions (SDF) were calculated for each region on the source image, and the weights \( w_{i} \) associated to each local transformation were defined using a logistic regression function applied to each SDF. However, the weights were estimated statically based on the source image and were not updated at each temporal integration step in Eq. (5). This means that the areas where rigid objects could be transformed without losing their rigidity were predefined. To overcome this limitation, we propose a new weighting scheme that adapts to changes at each temporal integration step. We define the weight function as

$$ w_{i}^{n} \left( \varvec{x} \right) = \frac{1}{{1 + \exp \left( {c D_{i}^{n} \left( \varvec{x} \right)} \right)}}, $$
(6)

where \( c \) controls the slope of the weight function at the region boundaries, \( D_{i}^{n} \left( \varvec{x} \right) \) is the SDF of the region associated to transformation \( i \) at coordinates \( \varvec{x} \) at the temporal integration step \( n \), and \( w_{i}^{n} \left( \varvec{x} \right) \) is the weight of the local transformation \( i \) at the temporal integration step \( n \). The transformation can then be written as

$$ \varphi_{n} \left( \varvec{x} \right) = \left( {\sum\nolimits_{i = 1}^{M} {w_{i}^{n - 1} \left( {\varphi_{n - 1} \left( \varvec{x} \right)} \right)\varvec{T}_{i}^{{\frac{1}{{2^{N} }}}} } } \right)\varphi_{n - 1} \left( \varvec{x} \right). $$
(7)

2.4 Surface-Based Optimization

In [4], the registration framework was designed to reduce malformations (defined as the distance between the patient’s cranium and its closest normal shape from a multi-atlas), while curvature information was not considered. In our work, we introduce a global surface-based dissimilarity measure based on currents [9] to compare the patient’s cranium and its closest normal shape. We chose this kind representation because it quantifies shape differences by considering both distances and curvature.

To obtain a surface representing the cranial shape of the patient, we used the approach described in [2] based on the constrained relaxation of an embedding sphere. Having the transformed patient \( (S) \) and the closest normal meshes \( (C) \), their comparison in the current’s space can be written as

$$ \begin{aligned} M\left( {S,C} \right) & = \sum\nolimits_{f, g} {\varvec{n}_{\varvec{f}}^{T} k_{W} \left( {\varvec{c}_{\varvec{g}} , \varvec{c}_{\varvec{f}} } \right)\varvec{n}_{\varvec{g}} } + \sum\nolimits_{q, r} {\varvec{n}_{\varvec{q}}^{T} k_{W} \left( {\varvec{c}_{\varvec{q}} , \varvec{c}_{\varvec{r}} } \right)\varvec{n}_{\varvec{r}} } \\ & - 2\sum\nolimits_{f, q} {\varvec{n}_{\varvec{f}}^{\varvec{T}} k_{W} \left( {\varvec{c}_{\varvec{q}} , \varvec{c}_{\varvec{f}} } \right)\varvec{n}_{\varvec{q}} } \\ \end{aligned} $$
(8)

where \( f \) and \( g \) represent triangles of \( S \), \( q \) and \( r \) represent triangles of \( C \), \( \varvec{n}_{\varvec{f}} \) is the normal vector of triangle \( f \), \( \varvec{c}_{\varvec{f}} \) is the center of triangle \( f \), and \( k_{W} \) is a kernel function. In our implementation, we used an isotropic Gaussian kernel.

2.5 Surface Preservation and Dissimilarity Measure

Since the affine transformations combined in Eq. (7) can change the surface area of a bone piece, we constrained our framework to preserve bone piece surface areas by including the following surface preservation term in the cost function:

$$ {\text{A}}\left( {S,C} \right) = \sum\limits_{l} {\left( {\frac{{\mathop \sum \nolimits_{f\epsilon l} \left\| {\varvec{n}_{f} } \right\|^{2} }}{{\mathop \sum \nolimits_{q\epsilon l} \left\| {\varvec{n}_{q} } \right\|^{2} }} - 1} \right)^{2} } , $$
(9)

where \( l \) represents a labeled bone piece, \( f\epsilon l \) are the triangles of the transformed bone piece \( l \) in the patient’s cranial shape, and \( q\epsilon l \) are the triangles before transformation.

A regularizing smoothing term was also added to obtain smoother transitions between bone pieces. The final cost function can then be written as

$$ D\left( {S,C} \right) = M\left( {S,C} \right) + \beta A\left( {S,C} \right) + \gamma \sum\nolimits_{p\epsilon P} {\sum\nolimits_{q\epsilon Q} {\left( {1 - \frac{{\varvec{n}_{\varvec{p}} \cdot \varvec{n}_{\varvec{q}} }}{{\left\| {\varvec{n}_{\varvec{p}} } \right\|^{2} \left\| {\varvec{n}_{\varvec{q}} } \right\|^{2} }}} \right)^{2} } } , $$
(10)

where \( P \) represents the triangles in the boundaries between bone pieces, Q represents the neighboring triangles to \( P \), and \( \beta \) and \( \gamma \) are balancing parameters. Finally, a regular gradient descent optimizer was used to minimize the dissimilarity measure.

3 Experiments

We used CT images obtained from 15 patients (age 2.87 ± 2.58 months) with metopic craniosynostosis to evaluate the optimal surgical plan calculated with the proposed framework. The closest normal shape for each patient was obtained from a multi-atlas created from CT images of 100 healthy infants (age 5.80 ± 3.31 months, range 0–12 months old). Since the standard procedure for fronto-orbital advancement to treat metopic craniosynostosis only involves repositioning of the frontal bones and the supra-orbital bar [5], only the local transformations associated to bone pieces in the frontal bones were estimated, while the local transformation in the rest of the cranium was set to identity.

The value for the weighting term \( \gamma \) was set empirically to 10 times the initial value of \( M\left( {S,C} \right) \). The weighting term \( \beta \) was estimated based on the values reported in [2], where the authors showed a mean difference of malformations and curvature discrepancies of 50% on healthy subjects compared to metopic craniosynostosis patients in the frontal bones. Since \( A\left( {S,C} \right) \) is normalized to the interval \( \left[ {0,1} \right] \), and targeting a reduction of 50% on \( M\left( {S,C} \right) \), we estimated \( \beta \) using:

$$ \beta = 0.5 \times M\left( {S,C} \right) \times \frac{{{\text{frontal}}\,{\text{bones}}\,{\text{surface}}}}{{{\text{total}}\,{\text{cranial}}\,{\text{surface}}}} \times \frac{1}{\delta }, $$
(11)

where \( \delta \) is the allowed surface area change. Based on the agreement with our expert surgeons, we set this value to 0.5%, thus allowing only very small surface area changes as consequence of bone cuts and bending.

Both malformations and curvature discrepancies were calculated for the pre-operative and simulated post-surgical shapes obtained from our optimal plan, as well as for the normal cases of the atlas to have a normative reference. In addition, we calculated the bone surface in both pre- and post-operative cranial shapes, and the Von Mises stresses in the frontal bones (Young modulus 380 MPa, Poisson’s ratio 0.22) [10]. Table 1 shows the average results for the 15 patients included in this study. Moreover, Fig. 2 shows an example of the surgical plan for one representative patient with metopic craniosynostosis and trigonocephaly (triangular shape of the forehead).

Table 1. Cranial malformations and curvature discrepancies calculated for the 15 patients (pre-operative) and for their simulated outcome of our surgical plan (post-operative).
Full size table
Fig. 2.
figure 2

(a) Superior (left) and anterior (right) views of the malformations in the frontal bones for the same case shown in Fig. 1, both for its pre-operative cranial mesh representation (top) and its simulated post-operative cranial shape obtained with the optimal surgical plan (bottom). (b) Curvature discrepancies in the frontal bones for the same case. The rest of the cranium is shown in white. The black wireframe represents the closest normal shape from the multi-atlas.

Full size image

4 Results

The mean post-operative malformations obtained with the proposed method were within the ranges obtained for the normal cases in the multi-atlas (1.19 ± 0.87 vs. 1.48 ± 0.99 mm, for metopic craniosynostosis and normal cases, respectively), obtaining \( p = 0.13 \) using a Mann-Whitney test. The post-operative curvature discrepancies were similar to the curvature discrepancies obtained for the healthy cases (0.74 ± 0.58 vs. 072 ± 0.61 mm−1, for metopic craniosynostosis and normal cases, respectively), with \( p = 0.65 \). The maximum values of malformations and curvature discrepancies were reduced to normal ranges (3.65 ± 0.96 mm and 2.48 ± 0.72 mm−1 on craniosynostosis cases vs. 3.84 ± 2.10 mm and 2.79 ± 1.20 mm−1 on normal cases, for malformations and curvature discrepancies, respectively), obtaining p = 0.73 and p = 0.35 for malformations and curvature discrepancies, respectively.

Compared to the other state-of-the-art work in automatic surgical planning [4], our method obtained slightly larger reduction of mean malformations (51% vs. 49% in [4]). Moreover, we obtained a reduction of 35% in curvature discrepancies, which cannot be achieved by using the rigid framework proposed in [4].

We obtained a bone area surface reduction of 0.02% within the frontal bones, which is a much lower value than the clinically allowed surface area change as a result of bone cutting and bending (0.5%, see Experiments section). On the other hand, the total cranial surface was incremented by 0.58%. This reflects the goal of normalizing cranial shape and creating space for brain development by advancing the forehead bones in metopic craniosynostosis corrective [5, 6].

Finally, we evaluated the feasibility of the estimated optimal plan in terms of bone stress. We obtained an average Von Mises of 33.59 ± 19.18 MPa, which is lower than the maximum allowable stress of 87 MPa [10].

5 Conclusion

We presented a new surface-based diffeomorphic registration framework that allows matching surfaces including areas with different affine properties, while keeping the global consistency of the transformation. In addition, we demonstrated its applicability to the automatic creation of surgical plans for the treatment of metopic craniosynostosis. Our automated framework addresses the challenges in the state-of-the-art of computer assisted surgical planning for cranial vault remodeling by adding objective metrics and reproducible computational methods to calculate an optimal plan for surgery. Our method implements all the basic operations that are performed during surgical treatment (bone cutting, bending, repositioning) and our quantitative evaluation of the results demonstrate that the simulated post-surgical cranial shapes are within healthy ranges. In our future work, we will extend our bone cut template to allow for automatic surgical treatment planning of other types of craniosynostosis.