An Image Registration Framework for Discontinuous Mappings Along Cracks

Abstract

A novel crack capable image registration framework is proposed. The approach is designed for registration problems suffering from cracks, gaps, or holes. The approach enables discontinuous transformation fields and also features an automatically computed crack indicator function and therefore does not require a pre-segmentation. The new approach is a generalization of the commonly used variational image registration approach. New contributions are an additional dissipation term in the overall energy, a proper balancing of different ingredients, and a joint optimization for both, the crack indicator function and the transformation. Results for histological serial sectioning of marmoset brain images demonstrate the potential of the approach and its superiority as compared to a standard registration.

1 Introduction

Image registration for smooth displacements is a very well explored topic in the literature; see, e.g. [9, 13], and references therein. There also exists work on discontinuous deformations for applications such as lung motion compensation (sliding motion) or partially constrained registration (local rigidity); see, e.g. [12, 14] and [5, 7]. However, to our best knowledge, the stage for images showing cracks has not yet been set although cracks frequently appear in applications with very thin and fragile objects such as histological serial sectioning [15].

Histological sections are still a very important foundation of the analysis of organs and diseases and are often considered as ultimate ground truth. Tissue is being sectioned into very thin (1–20 $\upmu$m), essentially 2D slices which are then analyzed using high-resolution microscopes; see, e.g., [15] for the histological background and Fig. 3 for an example. Our objective is to regain the 3D information of the findings. This is achieved by registering the slices to compensate non-linear tissue deformation introduced by the sectioning process; see also [15]. With this fragile material, crack formation related to distortions, movements, and in particular drying processes is common [15].

An example is an animal model for brain behavior. The interest in brain models for studying diseases such as degeneration is enormous. Typically, primate species, such as marmoset monkey, are used to build these models [11]. These models obtain mesoscale neuro-anatomical information through the study of histological images. To this end, sliced brain tissue is stained with different chemical dyes to highlight key biomarkers and digitized to an appropriate image format. Afterward, the image stack has to be registered to correct for the deformations introduced by the sectioning process and to gain 3D information. The alignment of these sections is difficult, as the tissue undergoes several deformations steps during sectioning and preparation. The tissues may suffer from non-linear deformations including cracking, tearing, or splitting; see Fig. 3 for an example. The registration of histological serial section is hence an important task and well explored in the literature; see, e.g., [13, 15], and references therein.

The picture changes when cracks appear. One of the obstacles is that a solid mathematical model is non-trivial and still missing. We consider a crack as an interface along which material breaks apart. The formation of a crack is a non-smooth process and the deformation map attains a jump discontinuity at the crack interface [4]. This process develops a gap or a hole in the material. Note that a crack might be viewed as an opening of the material and registration may aim to mimic this formation process (opening) or its inverse (closing).

Borrowing ideas from damage mechanics [8], we present a mathematical model for the registration with discontinuous deformations along cracks. A first attempt is proposed in [6], where the global regularization parameter is made spatially-dependent. This enables to relax the regularity along crack interfaces, but it also requires a pre-segmentation of the crack. Spatially-dependent regularization has also been used in several studies [12, 14], mainly to handle discontinuity due to sliding motion. However, these studies assume that a one-to-one correspondence exists almost everywhere for the given images. But, this assumption is not valid for the images with cracks.

Our work differs from the approach in [6]. First of all, we assume a crack to be a set of measure zero. We find this a more natural assumption, but it also provides challenges when dealing with integral measures. Next, our framework estimates both the deformation map and the location of the cracks and holes. Therefore, additional segmentation is not required. Note that the identification of the crack and in particular, its origin is non-trivial; see also our estimated crack interfaces in the result section. The new model utilizes the physical principles underlying crack propagation in material as studied in damage mechanics; see, e.g., [8]. These principles provide insight into the deformation processes. A damage mechanic model enables discontinuities in the deformation map by varying the stiffness of a material. This approach is thus similar to the spatially-varying regularization idea used in image registration. But the damage model also assumes that the material dissipates energy during the formation of a crack. This energy model the properties of a crack in terms of a crack indicator function. Finally and similarly to [3], we also introduce an adaptive spatially-varying dissimilarity weight to handle non-correspondences due to cracks and/or holes.

The proposed framework belongs to the class of joint segmentation and registration problems, where both an indicator function and the deformation map are unknown. Such frameworks are also explored for image denoising and fracture mechanics problems. In particular, the approach of Mumford and Shah [10] is common. We use ideas similar to the Ambrosio and Tortorelli [1] approximation of the Mumford Shah functional to avoid numerical difficulties arising from the low dimensional cracks. More specifically, we follow a damage mechanics framework [8], which is closely related to the Ambrosio and Tortorelli approximation.

In this paper, the focus is on the mathematical model and its properties and less on a particular application and parameterization. Therefore, we show results for 1D test case, where analytical solutions are known and can be used for objective comparison and parameter estimation. We also show preliminary results for marmoset brain images for which parameters are hand-tuned. The results confirm our expectations: the new approach is much better suited for the registration of images with cracks than a generic approach.

The paper is organized as follows. In Sect. 2, we outline the generic registration approach and describe our modifications to handle images with cracks. We introduce the new model and its parameterization. Sect. 3 comments on the numerical scheme to solve the new model. We demonstrate the superiority of our model on the 1D and 2D datasets in Sect. 4 and conclude our findings in Sect. 5.

2 Problem Statement

We begin by briefly outlining a generic image registration approach, which is based on a variational formulation of the correspondence problem; see [9] for details. As our focus is on images with cracks, we state the framework in its simplest form and even restrict to the 1D case. We stress that our concepts are generic and cover the general setting.

We next describe our modifications that enable us to handle images with cracks. The central idea is to decompose the image domain into a crack part $\varGamma$ and its complement $\varGamma ^c$. Following ideas from damage mechanic, we add a dissipation energy with is essentially based on the $(d-1)$-dimensional volume of the crack. Using ideas similar to the Ambrosio and Tortorelli approximation, we end up with our new model (3). Finally, we discuss our parameterization. We propose a particular parameter setting, such that only two parameters (regularization and crack weight) need to be specified.

Fig. 1.

A 1D function $R=\chi _{[-a,a]}$ (left) is mapped to T (center) by $\widehat{y}$ (right). Due to the discontinuity at the crack $\varGamma =\{0\}$ (in blue) a gap $\varSigma$ (in red) has been introduced. The forward mapping y with $x'=y(x)$ represents the opening of a gap, its inverse z with $x=z(x')$ represents the closing. Note that $z(\varSigma )=\varGamma$ but $y(\varGamma )$ is not well-defined. (Color figure online)

We remark that the inverse mapping z may provide additional insight; see Fig. 1. Due to page limitations, we do not explore this mapping further.

2.1 Generic Registration Problem

We only briefly outline the variational image registration approach; see e.g. [9] for details. In our setting $\varOmega \subset \mathbb {R}$ is a bounded domain, say $\varOmega =]-1,1[$, and two images $T,R:\varOmega \rightarrow \mathbb {R}$ are to be registered; see Fig. 1 for an example. More precisely, we want to determine a function $y:\mathbb {R}\rightarrow \mathbb {R}$ such that ideally $T[y]\approx R$, i.e. $T[y](x):=T\circ y(x):=T(y(x))\approx R(x)$ for almost all $x\in \varOmega$. For some applications, y might be assumed to be diffeomorphic, but in the presence of cracks this assumption does not hold.

The similarity of T[y] and R is enforced by minimizing a suitable distance D of the images. As the minimization of D is ill-posed, regularization is required. As examples serve an $L_2$ distance and a hyperelastic potential for 2D and $\phi =((v-1)/v)^2$; see [2] for details:

$$\begin{aligned} D(y):= & {} \textstyle \Vert T[y]-R\Vert _{L_2(\varOmega )}^2=\int _\varOmega \ [T(y(x))-R(x)]^2\ dx,\\ S(y):= & {} \textstyle \int _\varOmega \ \mathrm {hp}(y(x)) \ dx; \qquad \mathrm {hp}(y(x)):=\sum _{i,j}|\partial _iy_j(x)|^2+\phi (\det \nabla y(x)); \end{aligned}$$

Weighting the data fidelity and the regularization with parameters $\alpha$ and $\beta$, the generic approach (G) is to minimize the total energy

$$\begin{aligned} \textstyle J^{\mathrm {G}}_{\alpha ,\beta }(y) :=\int _\varOmega \ \alpha \>[T(y(x))-R(x)]^2+\beta \>\mathrm {hp}(y(x))\ dx; \end{aligned}$$

(1)

for details and generalizations see [9] and references therein. Note that the specific ingredients are irrelevant and the special choices are only for presentation reasons; our framework is general.

The parameters $\alpha$ and $\beta$ control the weights of the parts. Typical choices are $\alpha =1$ and $\beta =\tau$, $\tau$ being the smallest parameter for user-specific criterion such as the deformation being one-to-one is still fulfilled. Spatially dependent parameters may also be used [3, 6]. In [3], it is even suggested to automatically adjust the penalty for the handling of non-corresponding structures in images.

In [6], a spatially dependent regularization parameter $\beta =\beta (x)$ is studied. We remark that well-posedness requires $\beta (x)\ge \varepsilon >0$ for all $x\in \varOmega$, but it also guarantees a minimizing element y for $J^{\mathrm {G}}_{\alpha ,\beta }$. Furthermore, we remark that the minimizer y is globally smooth and has no discontinuities.

From a physical point of view, the parameter $\beta (x)$ represents the elastic stiffness of a material at a point $x\in \varOmega$, where stiffness reflects the rigidity of the material. The less stiff the material is, the easier to deform or stretch it is. A close to zero stiffness allows extreme stretching of material, and due to this, a crack may form, or a hole may close in the material. This interpretation motivates us to define the function $\beta$ depending on the crack indicator p.

2.2 Registration for Images with Cracks

As outlined above, the generic framework is incapable of registering images with cracks. We, therefore, suggest extending the generic model as follows. The basic idea is to partition the domain $\varOmega$ into a crack $\varGamma$ and its complement $\varGamma ^c:=\varOmega \backslash \varGamma$. For the registration on $\varGamma ^c$, we essentially follow the generic approach, whereas the crack deserves special treatment.

In material science, a crack is an interface along which material breaks apart. The breaking may create a hole in the material. For illustration, we consider the previous example; see also Fig. 1. Here, R is characteristic function of an interval $[-a,a]$ whereas T is the characteristic function of two intervals. With the transformation $\widehat{y}$ in Fig. 1, it holds $T(\widehat{y}(x))=R(x)$ on $\varGamma ^c$. Obviously, there is no continuous one-to-one mapping y that creates the gap $\varSigma$ between the two parts in T. Paying the price for a discontinuous transformation and introducing a crack at $\varGamma =\{0\}$, we see that $\widehat{y}$ is a perfect solution, as it is essentially a translation on $\varGamma ^c$ and thus does have minimal regularization energy.

We may describe a crack set $\varGamma \subset \varOmega$ by its indicator or phase field function

$$\begin{aligned} \textstyle p:\varOmega \rightarrow [0,1], \quad p(x)=\chi _{\varGamma }(x) \quad \text{(i.e. } \text{1 } \text{ on } \varGamma \text{, } \text{0 } \text{ on } \varGamma ^c\text{) } \end{aligned}$$

(2)

The closed subset $\varGamma$ has Lebesgue measure zero, defines the crack interface, and is a set of jump discontinuities for y; see [4] for details. Since $\varGamma$ has Lebesque measure zero, information for y on $\varGamma$ is not required.

If the crack and thus p is known a priori, we would simply pick the parameters $\alpha$ and $\beta$ in (1) depending on p, such that $\alpha (p)=\beta (p)=0$ for $p=1$, and we end up with a generic registration problem on $\varGamma ^c$. However, in our new approach, we aim to estimate the transformation and the crack simultaneously. To this end, we add an additional dissipation energy for the cracks in (1). The overall energy, which is now to be minimized with respect to (y, p) thus reads

$$\begin{aligned} J^{\mathrm {C}}_{\alpha ,\beta ,\gamma }(y,p) :=\int _\varOmega \ \alpha (p)\>[T(y(x))-R(x)]^2+\beta (p)\>\mathrm {hp}(y(x))+ \gamma \>e(p)\ dx. \end{aligned}$$

(3)

Ideally, $\int _\varOmega e\ dx=\int _\varGamma \ ds$ would measure the $(d-1)$-dimensional volume of the crack. This is similar to the edge set in the approach of Mumford and Shah [10] and known for its challenges. Ambrosio and Tortorelli [1] proposed to use a converging sequence of approximations, where p is differentiable. Roughly speaking,

$$ \textstyle \int _{\varGamma }\ ds \approx \lim _{h\rightarrow 0} \int _\varOmega \frac{1}{4h}p^2+h\>|\nabla p|^2\ dx. $$

In our approach, we modify this idea based on concepts from damage mechanics [8]. In damage mechanics, it is assumed that a material dissipates energy during the formation of a crack. This energy is proportional to the area of the crack region in the reference configuration of the domain. Generally, the dissipation energy is expressed as

$$\begin{aligned} \textstyle e(p) = p^2 + \ell ^2\>|\nabla p|^2, \end{aligned}$$

(4)

where $\ell$ is the internal length and relates to material properties; see [8] for details.

We remark that the concept of $\varGamma$-convergence or dissipation energy can be interpreted as a relaxation of the ideal setting. Rather than working with Dirac measures of the crack, we now have a differentiable phase field p. As a consequence, the proposed model (3) implicitly assumes that the deformation field is continuous and differential everywhere. This is good news as the iterates in our numerical schemes are smooth and differentiable. But it is also bad news, as the discontinuity can only be reached in the limit. In our implementation, we work with a fixed $\gamma$ and accept an over-smoothed solution.

Our approach is similar to ideas discussed in damage mechanics [8] to handle cracks and in image registration to handle sliding motion [14]. There, the non-smoothness is generally achieved by decreasing the regularization parameter $\beta$ in the generic case. We remark that our approach is conceptually capable of dealing with sets of measure zero, while the approaches [6, 14] are not.

2.3 Parameter Selection Strategy

We now discuss options for the parameters $\alpha$, $\beta$ and $\gamma$ in (3). The weights $\alpha$ and $\beta$ should ideally be related to indicator functions of the complement of the gap. For our numerical studies, we restrict to the basic choices,

$$\begin{aligned} \textstyle \alpha (p)= & {} (1-p)^2, \qquad \beta (p)=\varepsilon +\lambda (1-p)^n,\quad n\in \mathbb {N},\ \varepsilon ,\lambda>0,\\ e(p)= & {} p^2, \quad \text{ i.e. } \ell =0 \text{ in } \text{(4), } \text{ and } \gamma >0 \text{ constant }. \end{aligned}$$

These choices are convex and monotonic functions for $\alpha$ and $\beta$. The similarity parameter $\alpha$ is simply a smooth characteristic function for $\varGamma ^c$. The regularity parameter $\beta$ is along the same lines. However, we add a global and typically small constant $\varepsilon >0$, which ensures the existence of solutions for (3). Our parameter $\lambda$ plays a similar role as the constant parameter $\tau$ for the generic case (1) and in our numerical studies, we use the same strategy for picking $\lambda$ and $\tau$. The power n enables different distributions of the regularization weight. The higher n, the more concentrated is the regularization on the crack complement. As common in damage mechanics [8], we picked $n=2$ in our experiments.

For the dissipation energy e we again follow [8] and use (4) with the non-physical choice $\ell =0$, simply to keep the number of tunable parameters small. We remark that we, therefore, have no regularization on the first variation of p on the crack. Physically, the parameter $\gamma$ reflects the toughness of the material. The less tough the material is, the easier it is to break it. If $\gamma$ is very big, the energy converges to the generic case, i.e., no cracks at all. On the other side, if $\gamma =0$, a solution will be $p(x)=1$ (all is crack) and y is such that $S(y)=0$. The impact of the weight $\gamma$ for the dissipation is discussed in our results part.

We remark that for $p\equiv 0$ and proper choices of $\alpha$ and $\beta$, the new model coincides with the generic model and is thus a generalization. With this setting, the only degrees of freedom in parameterization is the choices of $\lambda$ and $\gamma$.

3 Discretization and Numerical Optimization

The minimization problem (3) is solved numerically, using the “discretize-then-optimize” paradigm as outlined in the FAIR toolbox [9]. More precisely, the integrals are approximated by quadrature rules, where for the images a cell-centered grid is used. Our discretization of the transformation y depends on the regularization; see [9] for details. In our experiments, we use an elastic regularizer [9] for 1D problems and the hyperelastic regularizer in [2] for the 2D problems, and therefore a nodal discretization for y. The crack indicator function p is discretized on a cell-centered grid.

As common for combined problems, we use an alternating minimization approach. Starting with $p\equiv 0$, we solve for y. This is equivalent to solve the generic problem (1). Using the latest y, we improve p and continue this iteration until our stopping criteria are met. Each sub-problem is solved with the Gauss-Newton scheme using an Armijo linesearch and typical stopping rules; see [9].

We remark that in contrast to the generic case, we do not take full advantage of the FEM approach to compute the hyperelastic regularization energy. The reason is that due to factor $\beta$, we do not have an analytic expression for the integral. We use a midpoint quadrature rule for the FEM based hyperelastic potential instead. We also remark that a joint optimization approach for (y, p) may converge to local minima, as typical for combined problems.

Following [9], we use a generic multilevel strategy to smooth out local minima for both 1D and 2D problems. We know that this simplistic strategy may not be optimal in the crack scenario but do not have space to elaborate on the issue.

4 Results

We demonstrate the performance of our new model. The experiment on a 1D example shows the convergence properties of our scheme and manifest that the new model can resolve discontinuous transformations perfectly along the cracks. This is further reaffirmed by our remarkable registration results on a pair of marmoset brain images; the data courtesy of Harald Möller, Max Planck Institute for Human Cognitive and Brain Sciences, Leipzig, Germany.

In all experiments, we use the $L_2$ or SSD based similarity measure [9]. We choose the elastic regularizer for the 1D studies and the hyperelastic regularizer in [2] for the 2D data. We compare results from the crack capable model (3) with the generic model (1) (denote with subscript ${\mathrm {C}}$ and $\mathrm {G}$, respectively).

Fig. 2.

Detailed results for the 1D performance study; cf. Fig. 1. All functions are symmetric, thus we only show the non-negative parts. Results for the forward pair (y, p) (top) and backward/inverse pair (z, q) (bottom) are displayed. On the left, we show the displacement $u(x)=y(x)-x$ of the solution $\widehat{y}$ (black dashed), the numerical solutions $y^\mathrm {G}$ and $y^{\mathrm {C}}$ for the generic (red dashed) and new approach (solid) for various discretization levels $h\approx 1/m$; $m=16$ (blue) $m=64$ (green), and $m=256$ (red). On the right we show the corresponding crack indicator functions p and q, respectively. Ideally, $p=\chi _{\{0\}}$, $q=\chi _{[-b,b]}$, and z as in Fig. 1. In our computations, we choose optimal parameter values, such that the estimates y and z have minimal error with respect to analytical solutions, i.e. $\alpha =1$, $\beta =0.06$ for the generic approach and $\alpha = (1 - p)^2$, $\beta =0.06(1-p)^2+10^{-9}$, $\gamma =0.02$ for the new approach. (Color figure online)

In our first study, we demonstrate that our results converge to the solution as the discretization error h reduces. To this end we use the pair of 1D functions as shown in Fig. 1, where an analytic solution $\widehat{y}$ is known. Figure 2 shows the displacement $u(x)=y(x)-x$ for the solution $\widehat{y}$ (and its inverse z) as well as numerical results for the generic and the crack capable approach on different discretization levels. It can be seen that the new approach converges to the true solution: for $m=256$, the gap is below the discretization error. Note that already on the coarse discretization ($m=16$), the new approach is superior to the generic approach on finest discretization ($m=256$).

Our second study demonstrates the superiority of the crack capable model on a dataset of marmoset brain images; see Fig. 3 and its caption for our parameterization. Although the schemes seem to be similar when looking at the closing of the crack induced gap in the transformed images T(y) and also in the difference images, in fact, they produce very different transformations. This is most pronounced on the gap (see top row). For the generic approach, $y^\mathrm {G}$ is globally smooth and distributes the gap to a rather thick layer; cf. $T(y^\mathrm {G})$.

On the other side, the crack capable model yields a solution precisely as expected. Its solution $y^{\mathrm {C}}$ shows a discontinuity which is numerically resolved within the discretization error. Moreover, $y^{\mathrm {C}}$ is very smooth on the crack complement. This is also manifested in the hyperelastic potential and the determinant of the Jacobian, where the maximal values of them are 350 and 8 times higher, respectively, with respect to the generic approach.

Also, the crack indicator function p matches our expectations and resolves the crack as a very thin line. We remark that a number of additional areas have automatically been indicated as cracks. In the absence of ground truth, we can neither confirm nor reject these findings. Note that it is generally not easy to determine a crack manually. Therefore, pre-segmentation based registration approaches such as the one in [6] may fail. We remark that the regularization potential and the determinant of the Jacobian from the generic model may also be used as a crack indicator; this is subject of ongoing investigations.

Fig. 3.

Results for marmoset brain images. The rows show the reference and the template image (top), result for the generic (middle) and results for the new approach (bottom). For the results, we show the transformed image T(y), the unweighted distance $|T(y)-R|$, the regularity $\mathrm {hp}(y)$, and the determinant of the Jacobian $\det \nabla y$ (from left to right). The top row also displays details of the transformations as well as the crack indicator function p in the new model. As parameters, we use $\alpha =1$ and $\beta =500$ for the generic case and $\alpha =(1-p)^2$, $\beta =500(1-p)^2+10^{-9}$, $\gamma =10^3$ for the crack model, respectively.

5 Conclusion

A novel crack capable image registration framework is proposed. The approach is designed for registration problems suffering from cracks, gaps, or holes. The approach enables discontinuous transformation fields and also features an automatically computed crack indicator function and therefore does not require a pre-segmentation. The new approach is a generalization of the commonly used variational image registration approach [9]. New contributions are an additional dissipation term in the overall energy, a proper balancing of different ingredients, and a joint optimization for both, the crack indicator function and the transformation. The approach is very general and flexible and can be combined with a huge variety of image similarity measures and regularization strategies. It is also not limited to a particular class of applications. For a particular setting, we also propose a parameterization strategy. Results for histological serial sectioning of marmoset brain images demonstrate the potential of the approach and its superiority as compared to a standard registration. Future work address a deeper mathematical analysis of the framework and the exploration of potential applications such as brain shift registration problem.