## 1 Introduction

*et al.*, 2008; Oliveira and Tavares, 2012). When evaluating the efficiency of treatment or progress of the disease, pre- and post-treatment, CT scans must be made (for the same patient) and compared by aligning (registering) these two (or more) scans or particular slices from the scans.

*et al.*(2011), the problem of registering CT scans in a body atlas has been considered. It is also required for navigating automatically to certain regions of a scan or if sub-volumes have to be identified automatically. Various methods are developed for problems of such type (Emrich

*et al.*, 2010; Feulner

*et al.*, 2009; Fernández

*et al.*, 2014). An automated method is developed in Shi

*et al.*(2007) in order to identify the corresponding nodules in serial thoracic CT scans for interval change analysis. The method uses rib centrelines as a reference for the initial nodule registration. The rib anatomy is used as a reference point in the CT scan analysis. In Kindig and Kent (2013), a model is introduced to describe the centroidal path of a rib (i.e. the sequence of centroids connecting adjacent cross-sections) in terms of several physically-meaningful and intuitive geometric parameters in CT scans. This model addresses a critical need for the accurate characterization of rib geometry in the biomechanics literature. A six-parameter shape model of the human rib centroidal path using logarithmic spirals is proposed in Holcombe

*et al.*(2016).

*et al.*(2016) cannot be applied here. The problem arises in selecting a proper function defining a contour, bounded by the fragments of the rib bone in the slice. These fragments are a result of the cross-section of a bone with the transverse plane. The slices may be compared using the bone tissue areas from the cross-section of bones in the slice. The authors in Bilinskas

*et al.*(2015) and Bilinskas

*et al.*(2017) offer a cardioid-type curve defining the rib-bounded contour on the slice. However, it is not the only one possible. A snake-type curve may serve as an alternative (Kass

*et al.*, 1988), but computing of such a curve will face problems in the spine area.

*et al.*(2017), where a method for analysing transversal plane images, obtained by computer tomography scans, is presented. Such a mathematical model was created and the problem of approximation is solved by finding out the optimal parameters of the model in the least-squares sense. The authors of paper (Bilinskas

*et al.*, 2017) disclose the problems that appear in building the proper model. Only slices, where ribs are visible, are considered. The methods of analysis of this part of the body are important because many internal organs are located here: liver, heart, stomach, pancreas, lungs, etc. The model is flexible and describes the rib-bounded contour independently of the patient age, sex, and disease.

*et al.*(2002) report that several entries in the DICOM header are often imprecise or even completely wrong.

*et al.*, 1995). In our case, parameters of the mathematical model and of distribution of the bone tissue on the CT scan slice form a set of features describing a particular slice.

## 2 The Model

*et al.*(2017), a method is proposed for bone tissue segmentation and further development of the mathematical model of the tissue configuration in a particular slice. Let us denote a set of bone tissue pixels of the slice by $B=\{({b_{1i}},{b_{2i}}),\hspace{2.5pt}i=\overline{1,m}\}$. Bilinskas

*et al.*(2017) is not the only possibility to extract the bone tissue from the slice, see e.g. Banik

*et al.*(2010), Zhang

*et al.*(2012). The bone tissue is approximated by a mathematical model. The model consists of two parts. The first part is a modification of cardioid, the second one is a supplement to reduce the spine influence. The first part of the model is described as follows:

##### (1)

\[ x(\varphi )={x_{0}}+a\rho (\varphi )\cos \varphi \cos \theta -b\rho (\varphi )\sin \varphi \sin \varphi ,\]*s*defines the spine cave ‘strength’,

*c*is the ‘strength’ of subtrahend for breastbone,

*l*is the steepness of subtrahend for breastbone,

*a*and

*b*are horizontal and vertical zoom respectively,

*θ*is the rotation of the human body, and $({x_{0}},{y_{0}})$ is the starting point (the ‘spike’ point of the model, see Fig. 1).

##### (4)

\[ \begin{array}{l}\text{a})\hspace{1em}({x_{0}},{y_{0}}),\\ {} \text{b})\hspace{1em}\big({x_{0}}+({y_{0}}-{\min _{y}})\sin \theta ;\hspace{2.5pt}{y_{0}}-({y_{0}}-{\min _{y}})\cos \theta \big),\end{array}\]## 3 The Problem of Registration in CT Image Analysis

##### (5)

\[ k=\arg \underset{{A^{\prime\prime }_{j}}\in {\boldsymbol{A}^{\prime\prime }}}{\min }\mathit{dist}\big({A^{\prime }};{A^{\prime\prime }_{j}}\big).\]## 4 Data for Slice Registration

*φ*runs through the interval $[-\pi /2;3\pi /2)$ with a step $2\pi /n$, we get the sequence $\boldsymbol{C}=({C_{i}}=({x_{i}},{y_{i}}),\hspace{2.5pt}i=\overline{0,n-1})$ of points of the curve (1)–(2), where ${x_{i}}$, ${y_{i}}$ are defined by Eq. (1) and (2) respectively, ${x_{i}}=x(\frac{2\pi }{n}i)$ and ${y_{i}}=y(\frac{2\pi }{n}i)$. The length of the sequence $\boldsymbol{C}$ is

*n*because the second part of the model (a line segment) is not used here.

*n*groups of bone tissue points are formed. Model points $({x_{i}},{y_{i}})$, $i=\overline{0,n-1}$ have weights $W=({w_{0}},{w_{1}},\dots ,{w_{n-1}})$, where ${w_{i}}$ is the number of bone tissue points in the

*i*th group; the

*i*th group contains points closer to the model point $({x_{i}},{y_{i}})$ than $({x_{j}},{y_{j}})$, $\forall j\ne i$. Without loss of generality, further we will use $W=({w_{0}},{w_{1}},\dots ,{w_{n-1}})$ as normalized weights, where ${\textstyle\sum _{i=0}^{n-1}}{w_{i}}=1$.

*B*, the sequence $\boldsymbol{C}$ of discrete points of the curve of model (1)–(3), array

*M*of 8 parameters describing the mathematical model, and weights

*W*of model curve points. Registration is applied to two slices. Denote

*B*, $\boldsymbol{C}$,

*M*, and

*W*of the source slice by ${B^{\prime }}$, ${\boldsymbol{C}^{\prime }}$, ${M^{\prime }}$, and ${W^{\prime }}$, and that of target slice by ${B^{\prime\prime }}$, ${\boldsymbol{C}^{\prime\prime }}$, ${M^{\prime\prime }}$, and ${W^{\prime\prime }}$, e.g. ${W^{\prime }}=({w^{\prime }_{0}},{w^{\prime }_{1}},\dots ,{w^{\prime }_{n-1}})$ are the weights of the source slice model, and ${W^{\prime\prime }}=({w^{\prime\prime }_{0}},{w^{\prime\prime }_{1}},\dots ,{w^{\prime\prime }_{n-1}})$ are the weights of the target slice model.

## 5 Registration

### 5.1 Rotation Invariance

*θ*, describing the rotation of a patient with respect to the bed. This parameter indicates the rotation of the model curve about the point $({x_{0}},{y_{0}})$ as well. Rotation invariance is realized rotating the model curve about the point $({x_{0}},{y_{0}})$ by the angle $-\theta $. This procedure should be applied both to source and target slices. The revised parameters of the models become ${M^{\prime }}=\langle {s^{\prime }},{c^{\prime }},{l^{\prime }},{a^{\prime }},{b^{\prime }},0,{x^{\prime }_{0}},{y^{\prime }_{0}}\rangle $ and ${M^{\prime\prime }}=\langle {s^{\prime\prime }},{c^{\prime\prime }},{l^{\prime\prime }},{a^{\prime\prime }},{b^{\prime\prime }},0,{x^{\prime\prime }_{0}},{y^{\prime\prime }_{0}}\rangle $. Without loss of generality and seeking for simplicity of notation, we redefine $\boldsymbol{C}$ by the sequence of points after the rotation described above, where entire points from $\boldsymbol{C}$ are rotated by $-\theta $ about $({x_{0}},{y_{0}})$. Therefore, in the further text, the sequences after rotations of source and target slices are denoted as ${\boldsymbol{C}^{\prime }}=(({x^{\prime }_{i}},{y^{\prime }_{i}}),\hspace{2.5pt}i=\overline{0,n-1})$ and ${\boldsymbol{C}^{\prime\prime }}=(({x^{\prime\prime }_{i}},{y^{\prime\prime }_{i}}),\hspace{2.5pt}i=\overline{0,n-1})$, respectively.

### 5.2 Scale Invariance

*a*and

*b*(see Eqs. (1) and (2)). For comparison of two slices, the parameters

*a*and

*b*of these slices should be scaled. Large scale invariance may be attained varying

*a*and

*b*. Therefore, the pyramid technique (Burt, 1981) has no use here. Three options O1, O2, and O3 of the unification of scales are considered below. Options O2 and O3 are used when DICOM metadata tags are not precise or are even lost.

##### (6)

\[\begin{array}{l}\displaystyle {M^{\prime }}=\big\langle {s^{\prime }},{c^{\prime }},{l^{\prime }},{a^{\prime }}\cdot z,{b^{\prime }}\cdot z,0,{x^{\prime }_{0}},{y^{\prime }_{0}}\big\rangle ,\\ {} \displaystyle z=\mathit{width}\big({\boldsymbol{C}^{\prime\prime }}\big)/\mathit{width}\big({\boldsymbol{C}^{\prime }}\big),\\ {} \displaystyle \mathit{width}\big({\boldsymbol{C}^{\prime }}\big)=\underset{i}{\max }{x^{\prime }_{i}}-\underset{i}{\min }{x^{\prime }_{i}},\hspace{2em}\mathit{width}\big({\boldsymbol{C}^{\prime\prime }}\big)=\underset{i}{\max }{x^{\prime\prime }_{i}}-\underset{i}{\min }{x^{\prime\prime }_{i}}.\end{array}\]##### (7)

\[\begin{aligned}{}\mathit{area}\big({\boldsymbol{C}^{\prime }}\big)=& \frac{1}{2}\big({x^{\prime }_{0}}{y^{\prime }_{1}}-{x^{\prime }_{1}}{y^{\prime }_{0}}+{x^{\prime }_{1}}{y^{\prime }_{2}}-{x^{\prime }_{2}}{y^{\prime }_{1}}+\cdots +{x^{\prime }_{n-2}}{y^{\prime }_{n-1}}-{x^{\prime }_{n-1}}{y^{\prime }_{n-2}}\\ {} & +{x^{\prime }_{n-1}}{y^{\prime }_{0}}-{x^{\prime }_{0}}{y^{\prime }_{n-1}}\big),\\ {} \mathit{area}\big({\boldsymbol{C}^{\prime\prime }}\big)=& \frac{1}{2}\big({x^{\prime\prime }_{0}}{y^{\prime\prime }_{1}}-{x^{\prime\prime }_{1}}{y^{\prime\prime }_{0}}+{x^{\prime\prime }_{1}}{y^{\prime\prime }_{2}}-{x^{\prime\prime }_{2}}{y^{\prime\prime }_{1}}+\cdots +{x^{\prime\prime }_{n-2}}{y^{\prime\prime }_{n-1}}-{x^{\prime\prime }_{n-1}}{y^{\prime\prime }_{n-2}}\\ {} & +{x^{\prime\prime }_{n-1}}{y^{\prime\prime }_{0}}-{x^{\prime\prime }_{0}}{y^{\prime\prime }_{n-1}}\big).\end{aligned}\]### 5.3 Translation Invariance

*s*may slightly compensate it.

#### 5.3.1 Pointwise Comparison (PW)

*i*th source slice model point and the

*i*th target slice model point were as minimal as possible. The problem may be formulated as a least-squares one:

*ϕ*is

##### (10)

\[ \frac{d\phi }{d\Delta y}={\sum \limits_{i=0}^{n-1}}2\big({y^{\prime\prime }_{i}}-\big({y^{\prime }_{i}}+\Delta y\big)\big)\]##### (11)

\[ \Delta y=\frac{1}{n}{\sum \limits_{i=0}^{n-1}}\big({y^{\prime\prime }_{i}}-{y^{\prime }_{i}}\big),\]##### (12)

\[ \mathit{dist}\big({A^{\prime }};{A^{\prime\prime }}\big)={\sum \limits_{i=0}^{n-1}}\big(\big({x^{\prime\prime }_{i}}-{\big({x^{\prime }_{i}}+\Delta x\big)\big)^{2}}+\big({y^{\prime\prime }_{i}}-{\big({y^{\prime }_{i}}+\Delta y\big)\big)^{2}}\big)\]##### Fig. 2

##### Fig. 3

#### 5.3.2 Total Least Squares (TLS)

*i*th point of the source slice model should be compared not with the

*i*th point of the target slice model, but with the nearest point on this model. Therefore, the problem to search for optimal $\Delta y$ may be formulated as follows:

##### (13)

\[ \underset{\Delta y}{\min }\phi (\Delta y)={\sum \limits_{i=0}^{n-1}}\big({y^{\prime\prime }_{i}}-\big({\overline{y}^{\prime }}{\big({x^{\prime\prime }_{i}},{y^{\prime\prime }_{i}},\Delta x,\Delta y\big)+\Delta y\big)\big)^{2}}\]*n*sampled points and problem (13) was solved using one-dimensional search.

##### (14)

\[\begin{aligned}{}\mathit{dist}\big({A^{\prime }};{A^{\prime\prime }}\big)=& {\sum \limits_{i=0}^{n-1}}\big(\big({x^{\prime\prime }_{i}}-({\overline{x}^{\prime }}{\big({x^{\prime\prime }_{i}},{y^{\prime\prime }_{i}},\Delta x,\Delta y\big)+\Delta x)\big)^{2}}\\ {} & +\big({y^{\prime\prime }_{i}}-({\overline{y}^{\prime }}{\big({x^{\prime\prime }_{i}},{y^{\prime\prime }_{i}},\Delta x,\Delta y\big)+\Delta y)\big)^{2}}),\end{aligned}\]#### 5.3.3 Weighted Total Least Squares (WTLS)

##### (15)

\[\begin{aligned}{}\underset{\Delta y}{\min }\phi (\Delta y)=& {\sum \limits_{i=0}^{n-1}}\big(\big({y^{\prime\prime }_{i}}-({\overline{y}^{\prime }}{\big({x^{\prime\prime }_{i}},{y^{\prime\prime }_{i}},\Delta x,\Delta y\big)+\Delta y)\big)^{2}}\\ {} & \times \big({w^{\prime\prime }_{i}}-{\overline{w}^{\prime }}{\big({x^{\prime\prime }_{i}},{y^{\prime\prime }_{i}},\Delta x,\Delta y\big)\big)^{2}}\big),\end{aligned}\]##### (16)

\[\begin{aligned}{}\mathit{dist}\big({A^{\prime }};{A^{\prime\prime }}\big)=& {\sum \limits_{i=0}^{n-1}}\big(\big(\big({x^{\prime\prime }_{i}}-({\overline{x}^{\prime }}{\big({x^{\prime\prime }_{i}},{y^{\prime\prime }_{i}},\Delta x,\Delta y\big)+\Delta x)\big)^{2}}\\ {} & +\big({y^{\prime\prime }_{i}}-\big({\overline{y}^{\prime }}{\big({x^{\prime\prime }_{i}},{y^{\prime\prime }_{i}},\Delta x,\Delta y\big)+\Delta y\big)\big)^{2}}\big)\\ {} & \times \big({w^{\prime\prime }_{i}}-{\overline{w}^{\prime }}{\big({x^{\prime\prime }_{i}},{y^{\prime\prime }_{i}},\Delta x,\Delta y\big)\big)^{2}}\big).\end{aligned}\]#### 5.3.4 Weighted Ordinary Least Squares (WOLS)

##### (17)

\[ \underset{\Delta y}{\min }\phi (\Delta y)={\sum \limits_{i=0}^{n-1}}\big({y^{\prime\prime }_{i}}-\big(\overline{\overline{y}}{\hspace{0.1667em}^{\prime }}{\big({x^{\prime\prime }_{i}},\Delta x\big)+\Delta y\big)\big)^{2}}{w^{\prime\prime }_{i}}\cdot \overline{\overline{w}}{\hspace{0.1667em}^{\prime }}\big({x^{\prime\prime }_{i}},\Delta x\big).\]*i*th point $({x^{\prime\prime }_{i}},{y^{\prime\prime }_{i}})$ of the target model is on top or bottom of the model curve. $\Delta x={x^{\prime\prime }_{0}}-{x^{\prime }_{0}}$ as above. ${w^{\prime\prime }_{i}}$ is the weight of the point $({x^{\prime\prime }_{i}},{y^{\prime\prime }_{i}})$ with respect to the density of the bone tissue spread near to this point. ${\overline{\overline{w}}^{\prime }}({x^{\prime\prime }_{i}},\Delta x)$ is the weight of the point of the source model curve (shifted by $\Delta x$) at the abscissa point ${x^{\prime\prime }_{i}}$, dependently whether the

*i*th point $({x^{\prime\prime }_{i}},{y^{\prime\prime }_{i}})$ of the target model is on top or bottom of the model curve.

*i*runs from 0 to $n-1$, the points $({x_{i}},{y_{i}})$ lie on the top of the slice model curve starting from such smallest

*i*(denote it by ${i^{\ast }}$), where ${x_{i+1}}<{x_{i}}$, and ending with such

*i*(denote it by ${i^{\ast \ast }}$), where ${x_{i+1}}>{x_{i}}$. The bottom curve consists of all the remaining points and has two points $({x_{{i^{\ast }}}},{y_{{i^{\ast }}}})$ and $({x_{{i^{\ast \ast }}}},{y_{{i^{\ast \ast }}}})$ common to the top curve. The criterion whether the point $({x_{i}},{y_{i}})$ lies on the bottom or top of the slice model is illustrated in Fig. 5.

##### (18)

\[ \Delta y=\frac{{\textstyle\textstyle\sum _{i=0}^{n-1}}({y^{\prime\prime }_{i}}-\overline{\overline{y}}{\hspace{0.1667em}^{\prime }}({x^{\prime\prime }_{i}},\Delta x)){w^{\prime\prime }_{i}}\cdot \overline{\overline{w}}{\hspace{0.1667em}^{\prime }}({x^{\prime\prime }_{i}},\Delta x)}{{\textstyle\textstyle\sum _{i=0}^{n-1}}{w^{\prime\prime }_{i}}\cdot \overline{\overline{w}}{\hspace{0.1667em}^{\prime }}({x^{\prime\prime }_{i}},\Delta x)}.\]##### (19)

\[ \mathit{dist}\big({A^{\prime }};{A^{\prime\prime }}\big)={\sum \limits_{i=0}^{n-1}}\big({w^{\prime\prime }_{i}}-\overline{w}{\hspace{0.1667em}^{\prime }}{\big({x^{\prime\prime }_{i}},{y^{\prime\prime }_{i}},\Delta x,\Delta y\big)\big)^{2}}.\]## 6 Experiments

##### Table 1

${\varepsilon _{\mathit{mean}}}$ | σ |
${\varepsilon _{\max }}$ | |

O1 | 6.575 | 8.405 | 36.25 |

O2 | 9.102 | 11.627 | 47.50 |

O3 | 9.648 | 11.976 | 47.50 |

##### Table 2

${\varepsilon _{\mathit{mean}}}$ | σ |
${\varepsilon _{\max }}$ | |

O1 | 9.974 | 9.226 | 38.75 |

O2 | 10.052 | 14.607 | 60.00 |

O3 | 8.737 | 13.507 | 60.00 |

*σ*, and the maximum error ${\varepsilon _{\max }}$. The results are presented in Tables 1–4 for various strategies of translation invariance of the first pair of CT scans. The results of weighted ordinary least-squares with invariance option O3 on the second pair of CT scans with 2.5 mm slice thickness are presented in Table 5.

*I*is the intersection of ${B^{\prime }}$ and ${B^{\ast }}$, $|\hspace{0.1667em}\cdot \hspace{0.1667em}|$ is the cardinality of the set. Here the intersection is assumed as the pixel-wise logical AND operator of two binary images ${B^{\prime }}$ and ${B^{\ast }}$. The function $\mathit{dist}({B^{\prime }};{B^{\ast }})$ returns values from interval $[0;1]$. The experiments have shown that the examined combination of Pyramidal Implementation of the Lucas Kanade Feature Tracker and RANSAC method is not very accurate, it leads to very poor registration results, as shown in Table 6.

##### Table 3

${\varepsilon _{\mathit{mean}}}$ | σ |
${\varepsilon _{\max }}$ | |

O1 | 0.977 | 0.604 | 2.50 |

O2 | 0.469 | 0.657 | 2.50 |

O3 | 0.703 | 0.740 | 3.75 |

##### Table 4

${\varepsilon _{\mathit{mean}}}$ | σ |
${\varepsilon _{\max }}$ | |

O1 | 0.508 | 0.614 | 1.25 |

O2 | 0.326 | 0.549 | 1.25 |

O3 | 0.339 | 0.555 | 1.25 |