Topology-Preserving 3D Image Segmentation Based on Hyperelastic Regularization

Abstract

Image segmentation is to extract meaningful objects from a given image. For degraded images due to occlusions, obscurities or noises, the accuracy of the segmentation result can be severely affected. To alleviate this problem, prior information about the target object is usually introduced. In Chan et al. (J Math Imaging Vis 60(3):401–421, 2018), a topology-preserving registration-based segmentation model was proposed, which is restricted to segment 2D images only. In this paper, we propose a novel 3D topology-preserving registration-based segmentation model with the hyperelastic regularization, which can handle both 2D and 3D images. The existence of the solution of the proposed model is established. We also propose a converging iterative scheme to solve the proposed model. Numerical experiments have been carried out on the synthetic and real images, which demonstrate the effectiveness of our proposed model.

Appendices

Computation of A in (27)

$A = I_{3}\otimes (A_{1}^{T},A_{2}^{T},A_{3}^{T})^{T}$, $A_{1} = I_{(n_{3}+1)}\otimes I_{(n_{2}+1)}\otimes \partial _{n_{1}}^{1,h_{1}}$, $A_{2}=I_{(n_{3}+1)}\otimes \partial _{n_{2}}^{1,h_{2}}\otimes I_{(n_{1}+1)}$, $A_{3}=\partial _{n_{3}}^{1,h_{3}}\otimes I_{(n_{2}+1)}\otimes I_{(n_{1}+1)}$ and

$$\begin{aligned} \partial _{n_{l}}^{1,h_{l}} = \frac{1}{h_{l}}\begin{pmatrix} -1 &{} 1 &{} &{} \\ &{} \cdot &{} \cdot &{} \\ &{} &{} -1 &{} 1 \end{pmatrix}\in {\mathbb {R}}^{n_{l},n_{l}+1},\quad 1\le l \le 3. \end{aligned}$$

(47)

Here, $\otimes$ indicates Kronecker product.

Computation of $\varvec{s}(Y)$ and $\varvec{v}(Y)$ in (28)

In each tetrahedron $\varOmega ^{i,j,k,l}$, set $\mathbf{L} ^{i,j,k,l}(\varvec{x})= (L_{1}^{i,j,k,l}(\varvec{x}),L_{2}^{i,j,k,l}(\varvec{x}),L_{3}^{i,j,k,l}(\varvec{x}))= (a^{i,j,k,l}_{1}x_{1}+a^{i,j,k,l}_{2}x_{2}+a^{i,j,k,l}_{3}x_{3}+b_{1}^{i,j,k,l}, a^{i,j,k,l}_{4}x_{1}+a^{i,j,k,l}_{5}x_{2}+a^{i,j,k,l}_{6}x_{3}+b_{2}^{i,j,k,l},a^{i,j,k,l}_{7}x_{1}$$+a^{i,j,k,l}_{8}x_{2}+a^{i,j,k,l}_{9}x_{3}+b_{3}^{i,j,k,l})$, which is the linear interpolation for $\varvec{y}$ in the $\varOmega ^{i,j,k,l}$. Note that

$$\begin{aligned} \begin{aligned} \partial _{x_{1}} L^{i,j,k,l}_{1}&= a^{i,j,k,l}_{1}, \partial _{x_{2}} L^{i,j,k,l}_{1} = a^{i,j,k,l}_{2},\partial _{x_{3}} L^{i,j,k,l}_{1} = a^{i,j,k,l}_{3},\\ \partial _{x_{1}} L^{i,j,k,l}_{2}&= a^{i,j,k,l}_{4}, \partial _{x_{2}} L^{i,j,k,l}_{2} = a^{i,j,k,l}_{5},\partial _{x_{3}} L^{i,j,k,l}_{2} = a^{i,j,k,l}_{6},\\ \partial _{x_{1}} L^{i,j,k,l}_{3}&= a^{i,j,k,l}_{7}, \partial _{x_{2}} L^{i,j,k,l}_{3} = a^{i,j,k,l}_{8},\partial _{x_{3}} L^{i,j,k,l}_{3} = a^{i,j,k,l}_{9}. \end{aligned} \end{aligned}$$

(48)

Then the following approximation can be built:

$$\begin{aligned}&\int _{\varOmega } \alpha _{s}\phi _{w}(\mathrm {cof}\nabla \varvec{y})+\alpha _{v}\phi _{v}(\det \nabla \varvec{y})\mathrm {d}\varvec{x}\nonumber \\&\quad \approx \frac{h^{3}}{6}\sum _{i=1}^{n}\sum _{j=1}^{n}\sum _{k=1}^{n}\sum _{l=1}^{6}(\alpha _{s}\phi _{w}(s^{i,j,k,l})+\alpha _{v}\phi _{v}(v^{i,j,k,l})), \end{aligned}$$

(49)

where

$$\begin{aligned} s^{i,j,k,l} = \begin{pmatrix} a_{5}^{i,j,k,l}a_{9}^{i,j,k,l}-a_{6}^{i,j,k,l}a_{8}^{i,j,k,l} &{} a_{6}^{i,j,k,l}a_{7}^{i,j,k,l}-a_{4}^{i,j,k,l}a_{9}^{i,j,k,l} &{} a_{4}^{i,j,k,l}a_{8}^{i,j,k,l}-a_{5}^{i,j,k,l}a_{7}^{i,j,k,l} \\ a_{3}^{i,j,k,l}a_{8}^{i,j,k,l}-a_{2}^{i,j,k,l}a_{9}^{i,j,k,l} &{} a_{1}^{i,j,k,l}a_{9}^{i,j,k,l}-a_{3}^{i,j,k,l}a_{7}^{i,j,k,l} &{} a_{2}^{i,j,k,l}a_{7}^{i,j,k,l}-a_{1}^{i,j,k,l}a_{8}^{i,j,k,l} \\ a_{2}^{i,j,k,l}a_{6}^{i,j,k,l}-a_{3}^{i,j,k,l}a_{5}^{i,j,k,l} &{} a_{3}^{i,j,k,l}a_{4}^{i,j,k,l}-a_{1}^{i,j,k,l}a_{6}^{i,j,k,l} &{} a_{1}^{i,j,k,l}a_{5}^{i,j,k,l}-a_{2}^{i,j,k,l}a_{4}^{i,j,k,l} \end{pmatrix} \end{aligned}$$

(50)

and

$$\begin{aligned} \begin{aligned} v^{i,j,k,l}&=a^{i,j,k,l}_{1}a^{i,j,k,l}_{5}a^{i,j,k,l}_{9}+a^{i,j,k,l}_{2}a^{i,j,k,l}_{6}a^{i,j,k,l}_{7}+a^{i,j,k,l}_{4}a^{i,j,k,l}_{8}a^{i,j,k,l}_{3}\\&\quad -a^{i,j,k,l}_{2}a^{i,j,k,l}_{4}a^{i,j,k,l}_{9}-a^{i,j,k,l}_{1}a^{i,j,k,l}_{6}a^{i,j,k,l}_{8}-a^{i,j,k,l}_{3}a^{i,j,k,l}_{5}a^{i,j,k,l}_{7}. \end{aligned} \end{aligned}$$

(51)

In order to write (49) into a compact form, we construct $D_{l}, 1\le l\le 9$:

$$\begin{aligned} \begin{matrix} D_{1} = [M_{1},0,0], &{} D_{4} = [0,M_{1},0], &{} D_{7} = [0,0,M_{1}],\\ D_{2} = [M_{2},0,0], &{} D_{5} = [0,M_{2},0], &{} D_{8} = [0,0,M_{2}],\\ D_{3} = [M_{3},0,0], &{} D_{6} = [0,M_{3},0], &{} D_{9} = [0,0,M_{3}], \end{matrix} \end{aligned}$$

(52)

where $M_{1}$, $M_{2}$ and $M_{3}$ are the discrete operators of $\partial _{x_{1}}$, $\partial _{x_{2}}$ and $\partial _{x_{3}}$ respectively and how to construct them is shown in “Appendix C”. Then we define $\varvec{s}(Y)$ and $\varvec{v}(Y)$ as follows:

$$\begin{aligned} \begin{aligned} \varvec{s}(Y)&= \begin{pmatrix} D_{5}Y\odot D_{9}Y-D_{6}Y\odot D_{8}Y &{} D_{6}Y\odot D_{7}Y-D_{4}Y\odot D_{9}Y &{} D_{4}Y\odot D_{8}Y-D_{5}Y\odot D_{7}Y \\ D_{3}Y\odot D_{8}Y-D_{2}Y\odot D_{9}Y &{} D_{1}Y\odot D_{9}Y-D_{3}Y\odot D_{7}Y &{} D_{2}Y\odot D_{7}Y-D_{1}Y\odot D_{8}Y \\ D_{2}Y\odot D_{6}Y-D_{3}Y\odot D_{5}Y &{} D_{3}Y\odot D_{4}Y-D_{1}Y\odot D_{6}Y &{} D_{1}Y\odot D_{5}Y-D_{2}Y\odot D_{4}Y \end{pmatrix},\\ \varvec{v}(Y)&= D_{1}Y\odot D_{5}Y\odot D_{9}Y + D_{2}Y\odot D_{6}Y\odot D_{7}Y +D_{4}Y\odot D_{8}Y\odot D_{3}Y\\&\quad -D_{2}Y\odot D_{4}Y\odot D_{9}Y-D_{1}Y\odot D_{6}Y\odot D_{8}Y-D_{3}Y\odot D_{5}Y\odot D_{7}Y,\\ \end{aligned} \end{aligned}$$

(53)

where $\odot$ denotes the Hadamard product of two matrices. Furthermore, set the ith component of $\varvec{s}(Y)$ as

$$\begin{aligned} \varvec{s}(Y)_{i} = \begin{pmatrix} (D_{5}Y\odot D_{9}Y-D_{6}Y\odot D_{8}Y)_{i} &{} (D_{6}Y\odot D_{7}Y-D_{4}Y\odot D_{9}Y)_{i} &{} (D_{4}Y\odot D_{8}Y-D_{5}Y\odot D_{7}Y)_{i} \\ (D_{3}Y\odot D_{8}Y-D_{2}Y\odot D_{9}Y)_{i} &{} (D_{1}Y\odot D_{9}Y-D_{3}Y\odot D_{7}Y)_{i} &{} (D_{2}Y\odot D_{7}Y-D_{1}Y\odot D_{8}Y)_{i} \\ (D_{2}Y\odot D_{6}Y-D_{3}Y\odot D_{5}Y)_{i} &{} (D_{3}Y\odot D_{4}Y-D_{1}Y\odot D_{6}Y)_{i} &{} (D_{1}Y\odot D_{5}Y-D_{2}Y\odot D_{4}Y)_{i} \end{pmatrix}. \end{aligned}$$

(54)

Then we can see that $\varvec{s}(Y)$ and $\varvec{v}(Y)$ contain all approximated cofactors and determinants for all tetrahedrons.

Computation of $M_{1}$, $M_{2}$ and $M_{3}$ in (52)

We first investigate the linear approximation $L(x_{1},x_{2},x_{3}) = a_{1}x_{1}+a_{2}x_{2}+a_{3}x_{3}+b$ in the tetrahedron $V_{3}V_{4}V_{5}V_{7}$ (Fig. 1). Denote these 4 vertices of this tetrahedron by $V_{3} = \varvec{x}^{1,1,1}$, $V_{4} = \varvec{x}^{2,2,2}$, $V_{5} = \varvec{x}^{3,3,3}$ and $V_{7} = \varvec{x}^{4,4,4}$. Set $L(\varvec{x}^{1,1,1}) = y^{1,1,1}$, $L(\varvec{x}^{2,2,2}) = y^{2,2,2}$, $L(\varvec{x}^{3,3,3}) = y^{3,3,3}$ and $L(\varvec{x}^{4,4,4}) = y^{4,4,4}$. Substituting $V_{3},V_{4}$, $V_{5}$ and $V_{7}$ into L, we get

$$\begin{aligned} \begin{pmatrix} x_{1}^{1} &{} x_{2}^{1} &{} x_{3}^{1} &{} 1\\ x_{1}^{2} &{} x_{2}^{2} &{} x_{3}^{2} &{} 1\\ x_{1}^{3} &{} x_{2}^{3} &{} x_{3}^{3} &{} 1\\ x_{1}^{4} &{} x_{2}^{4} &{} x_{3}^{4} &{} 1 \end{pmatrix} \begin{pmatrix} a_{1} \\ a_{2} \\ a_{3} \\ b \end{pmatrix} = \begin{pmatrix} y^{1,1,1}\\ y^{2,2,2}\\ y^{3,3,3}\\ y^{4,4,4} \end{pmatrix}. \end{aligned}$$

(55)

Then eliminating b, we obtain

$$\begin{aligned} \begin{pmatrix} x_{1}^{1}-x_{1}^{4} &{} x_{2}^{1}-x_{2}^{4} &{} x_{3}^{1}-x_{1}^{4} \\ x_{1}^{2}-x_{2}^{4} &{} x_{2}^{2}-x_{2}^{4} &{} x_{3}^{2}-x_{2}^{4} \\ x_{1}^{3}-x_{3}^{4} &{} x_{2}^{3}-x_{2}^{4} &{} x_{3}^{3}-x_{3}^{4} \end{pmatrix} \begin{pmatrix} a_{1} \\ a_{2} \\ a_{3} \end{pmatrix} = \begin{pmatrix} y^{1,1,1}-y^{4,4,4}\\ y^{2,2,2}-y^{4,4,4}\\ y^{3,3,3}-y^{4,4,4} \end{pmatrix}. \end{aligned}$$

(56)

Set

$$\begin{aligned} C = \begin{pmatrix} x_{1}^{1}-x_{1}^{4} &{} x_{2}^{1}-x_{2}^{4} &{} x_{3}^{1}-x_{1}^{4} \\ x_{1}^{2}-x_{2}^{4} &{} x_{2}^{2}-x_{2}^{4} &{} x_{3}^{2}-x_{2}^{4} \\ x_{1}^{3}-x_{3}^{4} &{} x_{2}^{3}-x_{2}^{4} &{} x_{3}^{3}-x_{3}^{4} \end{pmatrix}. \end{aligned}$$

(57)

Then we have

$$\begin{aligned} \begin{pmatrix} a_{1} \\ a_{2} \\ a_{3} \end{pmatrix} = \frac{1}{\det }\begin{pmatrix} C_{11} &{} C_{21} &{} C_{31} \\ C_{12} &{} C_{22} &{} C_{32} \\ C_{13} &{} C_{23} &{} C_{33} \end{pmatrix} \begin{pmatrix} y^{1,1,1}-y^{4,4,4}\\ y^{2,2,2}-y^{4,4,4}\\ y^{3,3,3}-y^{4,4,4} \end{pmatrix}, \end{aligned}$$

(58)

where $\det$ is the determinant of C and $C_{ij}$ is the (i, j) cofactor of C. Since the domain $\varOmega$ has been divided into N voxels, in order to find all $a_{1}$ in the tetrahedron with the same position of each voxel, we can make it as follows:

$$\begin{aligned} \begin{pmatrix} a_{1}^{1} \\ \vdots \\ a_{1}^{N} \end{pmatrix} = \frac{1}{\det }(C_{11}(E_{3}Y-E_{7}Y)+C_{21}(E_{4}Y-E_{7}Y)+C_{31}(E_{5}Y-E_{7}Y)), \end{aligned}$$

(59)

where $E_{l},l\in \{3,4,5,7\}$ is a matrix which extracts the corresponding positions of the vertices. Set $G_{1} = \frac{1}{\det }(C_{11}(E_{3}-E_{7})+C_{21}(E_{4}-E_{7})+C_{31}(E_{5}-E_{7}))$. For other 5 tetrahedrons, we can also build $G_{l}, l\in \{2,...,6\}$. Then we get

$$\begin{aligned} M_{1} = \begin{pmatrix} G_{1}\\ \vdots \\ G_{6} \end{pmatrix}. \end{aligned}$$

(60)

Similarly, we can obtain $M_{2}$ and $M_{3}$.