Medical image registration with object deviation estimation through motion vectors using octave and level sampling

Abstract

Medical image analysis presents a significant problem in the field of image registration. Recently, medical image registration has been recognized as a helpful tool for medical professionals. Current state-of-the-art approaches solely focus on source image registration and lack quantitative measurement for object deviation in terms of loosening, subsidence and anteversion related to surgery. In this article, we have provided motion vectors for recognizing the object deviation, in addition to detecting and selecting the feature points. Firstly, the feature points will be detected using Hessian matrix determinants and octave and level sampling. Then the strongest feature points are selected which will be utilized for identifying the object deviation with respect to the reference image through motion vectors. The objective of this work is to combine image registration and temporal differencing to achieve independent motion detection. In comparison to state-of-the-art approaches, the proposed methodology achieves higher Information Ratio (IR), Mutual Information Ratio (MIR) and their lower bounds for image registration.

KEYWORDS:

• Medical image registration
• feature point detection
• feature point selection
• motion estimation
• motion vectors
• differential operators

1. Introduction

The term “image registration” refers to the method through which several data sets are transformed into a single coordinate system. Several images taken from different sensors, at various periods, depths, distances, and/or vantage points may constitute data. Registration is required for comparing or integrating data from these different measurements. Due to the vast range of applications to which image registration is applied and various types of degradations in source images, it is impossible to create a universal approach that is suited for all applications. Image registration methods should consider not only the assumed type of geometric deformation between source images, but also noise corruption, desired registration accuracy, radiometric deformations, and the properties of the data that are specific to the targeted application.

Image registration has numerous applications, including but not limited to automatic target recognition in the military [1], medical imaging, computer vision, and the collection and analysis of satellite imagery [2]. Image registration is a pre-processing step in all fusion techniques [3–6] where multiple modalities of registered images will be combined into a single image, object detection [7], object recognition [8,9], and image classification [10–12].

Medical image registration involves matching of feature points in images of the same patient recorded at different times (multitemporal analysis) for more accurate diagnosis, prognosis, treatment, and follow-up [13]. Clinical applications requiring spatial alignment of anatomical features such as monitoring tumour growth and healing therapy, treatment verification, and comparison of the patient’s data with anatomical atlases, make registration an indispensable part of medical image processing. Image registration is useful when multiple imaging modalities are being utilized to record anatomical body structure such as Magnetic Resonance Spectroscopy (MRS). Positron Emission Tomography (PET), Magnetic Resonance Imaging (MRI) [14], Computed Tomography (CT) [15], Single Photon Emission Computed Tomography (SPECT), Ultrasound (US).

To discuss the utilization of research advances into practice and to provide a fair standard across competing approaches, the Learn2Reg [16] evaluates intra- and inter-patient, mono- and multimodal medical image registration for 3D images. This evaluation covers a broad spectrum of anatomies (brain, abdomen, and thorax), modalities (ultrasound, CT, MR), and availability of annotations

In 2D/3D registration, the notions of two weighted histograms of gradient directions and second-order spatial histograms were introduced and their benefits and drawbacks were assessed [17]. The gradient-weighted histogram is responsive to transformations in rotation and scaling but not translation. It causes a complex registration procedure and restricts the algorithm to a narrow convergence range for translation parameters. For statistical feature extraction, the approach relies on a spatially weighted histogram of gradient directions.

An exhaustive survey of applications of Transformers to medical image processing, including segmentation, detection, classification, restoration, synthesis, registration, and the generation of clinical reports was done in [18]. The authors have created a taxonomy for each of these uses, pinpointed the unique difficulties of each with suggestions for how to overcome them, and brought attention to the most current developments in each area.

In [19], the featured surface matching method was shown to be effective for registering lung images. To deform and align a 3D volume image with a 3D surface image, the technique was applied to volume-to-surface deformable image registration. The surface feature is extracted by calculating the curvature variation, and the correspondence between two surface points is established using the matching scheme.

To perform 2D/3D femur registration, [20] suggests using the core decompression technique as a surgical example, with the goal of registering preoperative MRI with intraoperative 2D fluoroscopic images. The registration's ensuing transformation is then utilized to generate overlays of preoperative MRI annotations and planning data for use by surgeons as a type of visual guidance.

There are some hurdles during the registration process, namely low resolution of images [21,22], deformation, organ motion, identification of motion vectors for detecting object deviation from source images, images with different modalities and the noise of interventional images [23]. Motivated by the above limitations and the demand for faster algorithms we proposed a novel method for image registration in this paper, that not only registers but also finds the amount of deviation using motion vectors. In this article we have taken datasets involving hip replacement diagnosis, brain, retina, and other human body scans where image registration plays an important role. After surgery and during follow-up visits, the doctor should be able to determine with certainty whether or not the replaced metal in the hip inserted through surgery is in the correct position even years later. In this scenario, our proposed method aids medical professionals. Similar object detection was done for retinal and brain scans.

2. Related work

Feature detection plays a crucial part in image processing. In order to compare identical objects in two images, one of which may be misaligned, feature detection is helpful. Many methods for processing medical images [24] rely on feature detection, which entails extracting the features of a certain image automatically based on the object's unique content. The method described above was employed in a wide variety of computer vision applications [25–28]. Examples of computer vision include object recognition and object tracing, which involve registering the feature points of an image.

In this section, we formulated how differential operators can be used as feature points and first-order convolutional filters to calculate Hessain determinants which in turn are used to find feature points.

2.1. Differential operators for feature points

If an image is focused using transformations like shifting with some scale or otherwise rotating at a certain angle, then it is feasible for object identification in an image. The following three stages are commonly used to determine the correspondence between the pixels in an image, while there are many others that can be utilized. The first stage is interest points selection from various locations like corners, globules, and intersections in an image. The second stage involves determining the feature vector from the surrounding pixels of each point of interest under various observational scenarios [29]. The last stage involves calculating the motion vectors by comparing the feature vectors from the two source images. In order to execute this kind of matching, the distance between the vectors must be determined. The correspondence between pixels can be determined with the help of standard geometrical transformations like the Euclidean distance.

2.2. Convolutional filters

One of the simplest linear filters in the spatial domain is the box filter. It replaces a pixel in an image with its average value of neighbourhood pixels [30], and this is a simple average method performed on pixels. Box filters are used to find the average value of neighbourhood 4-connectivity or 8- connectivity pixels.

Convolutional filters working with masks of size 3*3 or 9*9 sometimes are called convolution-symmetric kernels [31]. Let the coordinates of the image be (x,y) and “L” be the scale parameter which can characterize the size of the filter. The general property of the convolutional filter is given in Equation (1). (1) $\mathrm{\forall }\left(x,y\right)\in Z^2,\left(f\ast g\right)\left(x,y\right)=\sum _{\left(i,j\right)\in \mathrm{\Omega }}f\left(x-i,y-j\right)g\left(i,j\right)$(1) where $f$ represents an image in spatial domain, $g$ is a mask of size 3*3 used in our proposed method, $\left(x,y\right)$ denotes spatial coordinates.

Discrete differential operators are used in convolutional filters with several scales which are used to improve the filter efficiency and perform parameterized operation with the suitable kernel by above-mentioned variable “L”. To do this operation, first-order and second-order convolutional filters can be used.

2.2.1. First-order convolutional filter

Let us consider a digital image “u” and this can be operated with the true values of the operator at a scale $L$ accomplished with prime coordinate referred to as “D” multiplied with $D_x^L$. It creates the convolutional filter with a suitable size given by Equation (2). (2) $\ell \left(\text{L}\right)=0.8L\in \mathit{N}$(2)

Recollecting the information where the symbol [.] represents the round-off operation to the next higher value of integer. The first-order convolutional filter operator is given by Equation (3). (3) $\begin{array}{rl}{D}_x^{L}u:& =\left(g_{\left[-l,-1\right]\times \left[-l,l\right]}-g_{\left[1,l\right]\times \left[-l,l\right]}\right)\ast u\\ D_y^{L}u:& =\left(g_{\left[-l,l\right]\times \left[-l,-1\right]}-g_{\left[-l,l\right]\times \left[1,l\right])}\right)\ast u\end{array}$(3) where $D_x^{L}u$ and $D_y^{L}u$ are responses of first-order filter.

These filters use first-order partial derivatives as operators in a preference scale ℓ. We can get a good suitable formula for its impulse of “D” multiplied with $y^l$ operator. (4) $D_y^L\delta \left(x,y\right)=\{\begin{array}{ll}1& \text{if}\left(x,y\right)\in \left[-l,l\right]\times \left[-l,-1\right]\\ -1& \text{if}\left(x,y\right)\in \left[-l,l\right]\times \left[1,l\right]\\ 0& \text{otherwise}\end{array}$(4)

Computational complexity calculation using the integrated image as explained in Equation (2) we can pre-determine the integrated image say “U”, seven consecutive additions are enough to evaluate those operators irrespective of size of the parameter ℓ. Now let us consider and evaluate [32] the same Equation (2) with “b” and “a” equal to ℓ and c = -ℓ, d = −1.

We can get, the final filter coefficients using Equation (5),

(5) $\begin{array}{rl}{D}_y^{L}u\left(x,y\right)& =U\left(x+l,y+l\right)+U\left(x-l-1,y\right)\\ & -U\left(x+l,y\right)-U\left(x-l-1,y+l\right)\\ & -U\left(x+l,y-1\right)-U\left(x-l-1,y-l-1\right)\\ & +U\left(x+l,y-l-1\right)+U\left(x-l-1,y-1\right)\end{array}$(5)

3. Proposed method for medical image registration

The proposed image registration process attempts to discover matching feature points between two images and spatially align them. There are five major steps as depicted in Figure 1, that the proposed method has to go through for image registration. These could be listed as follows:

• Step 1: Estimating differential operators for the source images.

• Step 2: Feature points detection.

• Step 3: Feature points selection.

• Step 4: Feature points matching.

• Step 5: Object deviation estimation through motion vectors

Figure 1. Schematic of proposed method.

Each step is explained in detail in the following sub-sections.

3.1. Estimating differential operators for the source images

Let the two source images be $U_x$ and $U_y$. Differential operators will be calculated separately for each source image by using a convolutional filter Equation (3).

3.2. Feature points detection

In Section 2, the structures of first-order convolutional filters and how these filters are configured with their parameters were explained. In this section, feature points through multi-scale adaptive local feature detection are found. Convolutional filter [33] produces the response with scale-invariant property, i.e. the filter response should not vary with its operator’s scale. The term “multi-scale” is introduced here because, for multi-modal source images, multi-scale feature point detection can be employed.

Feature detection follows the steps shown below in Algorithm 1. Here, the inputs are $U_x$ and $U_y$ and it produces a list of feature points from the source images.

Algorithm 1 for detecting feature points

Step-1: Initialization of two differential operators $D_x^{L}u$ and $D_y^{L}u$

Step-2: For loop starts with L = 1 to 2

    Calculate Hessian matrix

   End For loop

Step-3: Select feature points

Step-4: For loop starts here to perform Octave Sampling

    For loop starts here to perform Level Sampling

      Calculate $2^{\mathrm{oi}}$

      List feature points by adding to Hessian matrix coefficients

    End for loop

  End for loop

Step-5: Return all feature points

In Algorithm 1, $2^{\mathrm{oi}}$ indicates level sampling, we can detect feature points and some important terms related to this algorithm like Hessian matrix coefficients, and active samplers as explained below. Second-order partial derivatives are arranged in a square matrix, which forms a Hessian matrix. Usually, the Hessian matrix uses the functional determinants. An example of a Hessian matrix applied on an MRI image is shown in Figure 2. For the leftmost image of Figure 2 [34] when the Hessian operator is applied, the edges will smoothen as in a centre image of Figure 2. Further smoothening of the image can be obtained as shown in the rightmost image by increasing the order of the Hessian operator.

Figure 2. An example of Hessian matrix applied on MRI image.

Hessian determinant calculation using convolution filters: The scale-independent Hessian matrix determinants, $D_O{H}^L\left(u\right)$ can be calculated from Equation (6). (6) $D_O{H}^L\left(u\right)=\frac{1}{L^4}\left(D_{\mathrm{xx}}^{L}u.D_{\mathrm{yy}}^{L}u-{\left(w{D}_{\mathrm{xy}}^{L}u\right)}^2\right)$(6) where $D_{\mathrm{xx}}^{L}u$ and $D_{\mathrm{yy}}^{L}u$ are calculated from Equation (3). As we mentioned in Algorithm 1, scaling can be achieved from the relation $L=2^{o\left(i+1\right)}$ and here the constant predefined weight “w” is 0.912. Normalization factors are required to get the scale invariance property, generally, we will take $\frac{1}{L^4}$ as a normalization factor. While calculating the Hessian determinant, the numerical tolerance will come into the picture to overcome this weighting factor “w” which was used.

3.3. Feature points selection

After detecting the feature points, we have to select only a few of them which are adaptively considered matched with the reference image. This can be achieved by considering the neighbourhood pixels. In this selection, the class of transformation which we used at the time of the feature detection procedure was important. By considering similarity parameters like location and shape, scale estimation is jointly performed with similarity transformation.

In this methodology, the common interest or desirable points are considered as the local maxima of the determinants of the Hessian matrix $D_O{H}^L$ of an image $\left(u\right)$. Now, the detection of maxima is performed [30] by considering P, Q, and R neighbourhood pixels, generally P = Q = R = 3, and comparing the box space with the nearest neighbourhood pixels. The above-said procedure is given in Algorithm 2. This algorithm consists of input arguments that are maxima of the determinants of the Hessian matrix $D_O{H}^L$ of an image $\left(u\right)$, octave sampling $\left(o\right)$, Level of interest $\left(i\right)$, $D_O{H}^L$ values at “$o$” and “$i$”. The output arguments are a list of key points.

Algorithm 2 for key points selection

Step-1: initialization of octave sampler, level of interest and hessian determinants

Step-2: Calculate $2^{o\left(i+1\right)}$

Step-3: For loop start with 0 to M-1

   Perform $2^{o-1}$ Correspondence.

   For loop start with 0 to N-1

     Perform $2^{o-1}$ Correspondence.

     Use Equation (7) to fix the threshold

     Maximum dominance calculation for $(D_O{H}^L{)}_{X,Y}$

     Interpolation parameters Calculation for $(D_O{H}^L{)}_{X,Y}$

   End for loop

  End for loop

Step-4: Return the list of key points

Procedure to set the threshold: Thresholding is an essential process after quantization with octaves. Here we use four octaves and two levels for analysis so that we can analyse 8 scales of representation. Using the method as in [35], it is not only a computationally efficient method of image representation but also provides noise immunity. We can get the key features by applying a threshold on Hessian determinants $\left(\mathrm{Do}{H}^L\right)$. (7) $\mathrm{Do}{H}^L\left(U_x\mathrm{and}{U}_y\right)>L$(7) We have to remember that this operator is scale-invariant and is constant. This normalization can be achieved using simple mathematical operations. By using the trial and error method too, we can set the threshold [36], let the threshold value be 103 and the reference image interval be [0,255]. Figure 3 describes setting the points of interest with reference to Hessian determinants $\left(\mathrm{Do}{H}^L\right)$ after the thresholding operation. The radii of the green circles are variable, in this application, the radii are set to 2.5 times the L (Box Scale).

Figure 3. Illustration of interest point detection.

3.4. Matching the key points

After selecting the feature points from two images, $U_x$ and $U_y$ match the feature points which are represented by their point of interest $\left(\{X_k{\}}_k\mathrm{and}\{Y_l{\}}_l\right)$ with their corresponding descriptors $\left(\mathrm{SURF}\mathrm{of}\right(\{X_k{\}}_k\mathrm{and}\mathrm{SURF}\mathrm{of}\{Y_l{\}}_l)$. We know that, for this type of descriptor, the interest points must satisfy the condition X ε [−1,1]⁶⁴ which indicates a vector with dimension 64 and its Euclidean norm is $||X||=1$. The comparison between the interest points was performed using Euclidean distance since these interest points are considered as vectors. Using a thresholding technique called Nearest-Neighbour Distance Ratio (NNDR), we can combine the features. Here, we calculate the Euclidean distance by comparing the descriptor $X_k$ with the image $U_x$ and $U_y$. (8) $d_{k,l}^2=X_k-{Y_l}_2^2=2\left(1-X_k^T.Y_l\right)$(8)

Let us consider one-megapixel image, almost millions of scalar products are computed from the second image “$U_y$”, this extends the computational complexity. Therefore, to lessen the computational complexity in the matching procedure, we need to compare the Laplacian [37] signs of key features. If the two descriptor’s signs are not equal, then those are considered as unlikely similar. In this case, there is no need to compare the key features and the comparison is discarded.

To match the key points, we need to find the most important and reliable correspondence from several number of reputed features, this leads to the best match of features. In the other way, if we consider the similarity distance $d_{k,l}$ only into account, it is necessary to separate the true match from the false match. The nearest neighbour matching method may or may not produce the best match because it considers only correspondence between $X_k$ and its most similar one, denoted by $Y_{l1}$ where $l_1$ is given in Equation (9), (9) $l_1\in {\arg }_l^{min}{d}_{k,l}$(9)

The Nearest Neighbour Distance Threshold is a proximity search technique which is used in pattern recognition, computer vision, etc., this is a straightforward method to find the correspondence by fixing the threshold on the similarity measure. Indeed, such a method of matching produces a negotiation between false positive numbers and false negative numbers. Therefore, instead of using this type of thresholding method [38], we used the Nearest Neighbour Distance ratio, which is a simple and reliable statistical method to discard mismatches.

In this method of matching, we need to find the ratio of the first and second nearest neighbour pixels and measure the correspondence between them. For doing this, first, we have to compute Y_ℓ1 from the desired X_k, such that (10) $l_2\in {\arg }_{l,s.t.d_{k,l}\ge d_{k,l_1}}^{min}{d}_{k,l}$(10)

Finally, compare the corresponding ratio to a suitable threshold “t” and follow the condition, if $\frac{d_{k,l1}}{d_{k,l2}}\le t^{\frac{d_{k,l1}}{d_{k,l2}}}$, then the correspondence was confirmed in this decision, and the threshold was set to t = 0.8, we may increase or decrease the threshold using trial and error methods. We are dealing with medical images hence the threshold was set to maximum value. After matching is complete, it is necessary to find the deviation of the matched points. This is explained in the below section motion vector calculation.

3.5. Motion vectors calculation

The motion estimation theory between two consecutive images using features is described in this section. The key idea of motion estimation is to determine [39] the displacement between two successive frames as shown in Figure 4.

Figure 4. Displacement between frames: (a) Reference image 1 and (b) Reference image 2.

In the simplest form, the problem can be described by Equation (11). This equation only permits the registration of translational movement in x and y and the intensity levels are expected to be the same all over the image. (11) $I\left(x,y,t\right)=I\left(x+\mathrm{\Delta x},y+\mathrm{\Delta y},t-\mathrm{\Delta t}\right)$(11)

Frame Translation: The general model of inconsistency between two objects in two respective frames is translation. Translation can identify the movement of any particular region in an image. Naturally, if an image is affected by translational movement, then it indicates objects [40] in the foreground have been moved. Suppose the translational movement is outside the reference image, then the new frame is created by the combination of both old coordinates plus translation coefficients. (12) $\left[\begin{array}{c}{x}_{\mathrm{new}}\\ y_{\mathrm{new}}\end{array}\right]=\left[\begin{array}{c}{x}_{\mathrm{old}}\\ y_{\mathrm{old}}\end{array}\right]+\left[\begin{array}{c}\mathrm{\Delta x}\\ \mathrm{\Delta y}\end{array}\right]$(12) where $x_{\mathrm{old}}$ and $y_{\mathrm{old}}$ are the old coordinates of the image in which the translation operation was performed, $x_{\mathrm{new}}$ and $y_{\mathrm{new}}$ are the new coordinates after performing the translation operation. The numerical values $\mathrm{\Delta x}$ and $\mathrm{\Delta y}$ are the small translational movement in X-axis and Y-axis direction, respectively.

Each pixel of the frame was translated and new coordinates were obtained from Equation (12). The entire frame was updated with new coordinates and the frame was neither damaged nor corrupted by noise. In the translation operation [41], all pixels of the original image are imposed in the output image, and the region of interest is moved to the desired position [38,42] using Equation (12). If all pixels of the original image are left unchanged, then the translated image is the same as the original image which means that there is no change in the output image. We can observe a slight difference between moving the pixel region and translating the pixel region if the original pixel region was filled with homogenous grey levels. An example of a translation operation is given in Figure 5.

Figure 5. Example for translational operation.

Since the pixel values of an image are integers, the translation procedure is straight forward and if we use sub-pixel values, we may use the bilinear interpolation method.

Rotation of an image: Spatial transformation techniques have a vital role in many image processing applications like CT, and MRI images. Rotational operation is one of the linear spatial transformation methods [43]. This type of operation depends on the direction of rotation and rotational angle.

Let us consider the camera movement to be an anti-clockwise rotation. Then the entire frame will be a clockwise rotation of all pixel values to a new location $(x_{\mathrm{new}},y_{\mathrm{new}})$ computed using Equation (13). (13) $\left[\begin{array}{c}{x}_{\mathrm{new}}\\ y_{\mathrm{new}}\end{array}\right]=\left[\begin{array}{c}\cos \theta -\sin \theta \\ \sin \mathrm{\theta cos}\theta \end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{old}}\\ y_{\mathrm{old}}\end{array}\right]$(13) where θ represents the rotational angle. For more reliable analysis, Equation (14) is written where the quantity $\overline{x}$ measures the average value. (14) $\left[\begin{array}{c}{\overline{x}}_{\mathrm{new}}\\ {\overline{y}}_{\mathrm{new}}\end{array}\right]=\left[\begin{array}{c}\cos \theta -\sin \theta \\ \sin \mathrm{\theta cos}\theta \end{array}\right]\left[\begin{array}{c}{\overline{x}}_{\mathrm{old}}\\ {\overline{y}}_{\mathrm{old}}\end{array}\right]$(14)

Equation (14) is rewritten as in Equation (15), (15) $\left[\begin{array}{c}{x}_{\mathrm{new}}-{\overline{x}}_{\mathrm{new}}\\ y_{\mathrm{new}}-{\overline{y}}_{\mathrm{new}}\end{array}\right]=\left[\begin{array}{c}\cos \theta -\sin \theta \\ \sin \mathrm{\theta cos}\theta \end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{old}}-{\overline{x}}_{\mathrm{old}}\\ y_{\mathrm{old}}-{\overline{y}}_{\mathrm{old}}\end{array}\right]$(15)

In numerous cases and more desirable cases involving medical applications [44], analysing the rotational movement of individual objects in the frame is needed. Equation (13) enables such type of mechanism. As mentioned in the above equations, we conclude the new coordinates for the rotated image will be as given in Equation (16). (16) $\begin{array}{rl}\left[\begin{array}{c}{x}_{\mathrm{new}}\\ y_{\mathrm{new}}\end{array}\right]& =\left[\begin{array}{c}\cos \theta -\sin \theta \\ \sin \mathrm{\theta cos}\theta \end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{old}}-{\overline{x}}_{\mathrm{old}}\\ y_{\mathrm{old}}-{\overline{y}}_{\mathrm{old}}\end{array}\right]\\ & +\left[\begin{array}{c}{\overline{x}}_{\mathrm{new}}-{\overline{x}}_{\mathrm{old}}\\ {\overline{y}}_{\mathrm{new}}-{\overline{y}}_{\mathrm{old}}\end{array}\right]\end{array}$(16)

4. Results and discussion

4.1. Experiment setup

Subjective and objective evaluations have validated the performance analysis of the proposed medical image registration using octave and level sampling. The simulation is run on a computer equipped with a 12th Gen Intel(R) Core(TM) i7-12700H processor, a RAM size of 16GB, and an NVIDIA GEFORCE RTX GPU. Figure 6 presents a total of 10 test samples [45–47] to be considered to evaluate the proposed method. The proposed medical image registration using octave and level sampling is compared with four state-of-the art methods: Scale Invariant Optical Flow [48] (SIOF), Constrained Least Squares [49] (CLS), Cost Function Minimization [50] (CFM), and Coherent Pixel Correspondence [51] (CPC). We have used the default values that the authors of these existing methods have provided. For quantitative analysis, we have used four quality metrics [52,53]: Information Ratio (Q_IR), Lower bound Information Ratio (Q_LIR), Mutual Information Ratio (Q_MIR), and Lower bound Mutual Information Ratio (Q_LMIR).

Figure 6. Ten medical test samples used for registration: (i) Test sample 1a (ii) Test sample 1b (iii) Test sample 2a (iv) Test sample 2b (v) Test sample 3a (vi) Test sample 3b (vii) Test sample 4a (viii) Test sample 4b (ix) Test sample 5a (x) Test sample 5b (xi) Test sample 6a (xii) Test sample 6b (xiii) Test sample 7a (xiv) Test sample 7b (xv) Test sample 8a (xvi) Test sample 8b (xvii) Test sample 9a (xviii) Test sample 9b (xix) Test sample 10a (xx) Test sample 10b.

4.2. Quality metrics

In this section, we used quantitative methods to analyse the effectiveness of the proposed algorithm. An image is described through a matrix of dimension $\mathrm{NXM}$ where each element is referred to as a pixel. Each pixel is itself a vector in $[2^D-1{]}^3$, where D is the colour depth of an image and each dimension represents grey scale intensity.

We assume that a set of $L$ images $X=\left\{X_1,X_{2,}\dots \dots .X_L\right\}$ is drawn from the distribution $P_X$ where $X_t\in \left[N\right]X\left[M\right]X[2^D-1{]}^3$. In the following, we refer to each $X_t$ as a frame. For a given $X_t$, $X_t\left(C\right)$ let $C\in \left\{\mathrm{graylevel}\right\}$ be the 2-D vector corresponding to the colour $C$; we refer to $X_t\left(C\right)$ as the component image of the colour $C$or the $C$ channel of the image $X_t$. Finally, $X_t\left(C,n,m\right)$ for $\left(n,m\right)\in \left[N\right]X\left[M\right]$ the pixel $\left(n,m\right)$ for the grey level $C$ at frame $t$

4.2.1. Information ratio (Q_IR)

Using Equation (17), $Q_{\mathrm{IR}}$can be calculated as the self-information of an image intensity relative to the information content of individual pixels with the same intensity. (17) $r\left(C,X\right)=Q_{\mathrm{IR}}=\sum _{i\in \left[2^D-\right]}{\overline{h}}_i\left(C\right)r_i\left(C,X\right)$(17) where ${\overline{h}}_i\left(C\right)\cong \sum _{n,m\in \left[N\right]X\left[M\right]}{1}_{\left\{X_t\left(C,n,m\right)=i\right\}},\mathrm{\forall i}\in \left[2^D-1\right]$, $r_i\left(C,X\right)=\{\begin{array}{l}\frac{-\log \left(P_i\left(C\right)\right)}{\text{log}\left({\overline{h}}_i\left(C\right)\right)},{\overline{h}}_i\left(C\right)>1\\ 0,\text{ otherwise}\end{array}$

$P_i\left(C\right)$ is self-information of the $i^{\mathrm{th}}$ level estimation as $-\log \left({\overline{h}}_i\left(C\right)/\mathrm{NM}\right)$

4.2.2. Mutual information ratio (Q_MIR)

Equation (18) is used to calculate Q_MIR which is analogous to that of the Q_IR, but in terms of a pair of intensity values over two images.

(18) $m\left(C,X_1;X_2\right)=Q_{\mathrm{MIR}}=\sum _{\left\{j,i\right\}\in \left[2^D-1\right]}{\overline{h}}_{i,j}\left(C\right)\frac{\log \left(\frac{P_{i,j}\left(C\right)}{P_i\left(C\right)P_j\left(C\right)}\right)}{\log \left({\overline{h}}_{i,j}\left(C\right)\right)}$(18)

4.2.3. Lower bound information ratio (Q_LIR)

Equation (19) gives the value of Q_LIR as a function of entropy and mutual information. The computational difficulty of calculating the feature count is greatly reduced by using this approximation. (19) $r\left(C,X\right)=Q_{\mathrm{LIR}}\ge \frac{\mathrm{NM}}{\log \left(\mathrm{NM}\right)}\hat{H}\left(C\right)\ge 0$(19) where $\hat{H}\left(C\right)$ is entropy.

4.2.4. Lower bound mutual information ratio (Q_LMIR)

Q_LMIR gives useful bounds between the entropy and the IR, and between the mutual information and the MIR. A lower bound on $m\left(C,X_1;X_2\right)$ is calculated as a function of $\hat{I}\left(C,X_1;X_2\right)$ as shown in Equation (20). (20) $m\left(C,X_1;X_2\right)={\mathit{Q}}_{\mathit{LMIR}}\ge \frac{\mathrm{NM}}{\log \left(\mathrm{NM}\right)}\hat{I}\left(C,X_1;X_2\right)\ge 0$(20) where $\hat{I}\left(C,X_1;X_2\right)$ is mutual information.

4.3. Results

In this section, we used ten test samples of medical images. The experiment with the proposed method of image registration not only performed feature detection but also identified how much the position of artificial bone inserted during surgery deviated from the reference image after surgery and during follow-ups. Figure 7(i) is a reference image of bone immediately after surgery and Figure 7(ii) is the image taken during follow-up visits say after a few months of surgery. From Equation (7), the strongest key points are selected as shown in Figure 7 (iii) and (iv). Here we can identify clearly that there was no deviation in the bone position, hence as seen in Figure 7(vi) no motion vectors exist. The same statement can be justified by the quantitative values of Figure 7(vii). Two polygons are shown in Figure 7(vii) for two reference images Figure 7(i) and (ii) have no slope which indicates that there is no deviation. In Figure 7(v) slope of the matched features is zero. The motion vector analysis performed on the test sample further justifies the proposed algorithm.

Figure 7. (i–vii) Hip replacement surgery without bone position deviation: (i) Test sample 1a (ii) Test sample 1b (iii) Feature points selection 1a (iv) Feature points selection 1b (v) Feature matching (vi) Motion vector analysis (vii) Object deviation estimation.

Consider Figure 8(i) which shows a reference image of bone immediately after surgery and Figure 8(ii) is the image taken during follow-up visits say after a few months of surgery. From Equation (7), the strongest key points are selected as shown in Figure 8(iii) and (iv). As it can be observed clearly in Figure 8(v) there exists a slope for each of the features matched, which indicates there is a mismatch in bone position between images of Figure 8(i) and (ii), hence motion vectors exist in Figure 8(vi). Two polygons shown in Figure 8(vii) for two reference images Figure 8(i) and (ii) have a slope which indicates that there is some deviation. If we compare Figure 7(vi) and Figure 8(vi), in the former there are no motion vectors because there is no deviation, but in the latter, we can observe a few motion vectors. This way we have three ways of identifying deviation namely feature matching, motion vector and polygon slope.

Figure 8. (i–vii) Hip replacement surgery with bone position deviation: (i) Test sample 2a (ii) Test sample 2b (iii) Feature points selection 2a (iv) Feature points selection 2b (v) Feature matching (vi) Motion vector analysis (vii) Object deviation estimation.

Our algorithm was tested not only on hip replacement surgery but also on other medical applications involving the retina, brain and a few other body scans. Similar analysis as in Figure 7 and Figure 8 was made for test samples Figure 6(v) to (xx) whose quantitative analysis is given in Figure 9 and Table 1. The proposed algorithm also identifies the amount of deviation like very small deviation/ no deviation/ small deviation/ high deviation as discussed in Figure 1. During the MATLAB simulation of the test samples in 1a and 1b of Figure 6, “no deviation” was observed while test samples 2a and 2b of Figure 6 showed “high deviation”.

Figure 9. Objective assessment of different image registration methods including proposed method: (i) Information Ratio, Q_IR (ii) Lower bound Information Ratio, Q_LIR (iii) Mutual Information Ratio, Q_MIR (iv) Lower bound Mutual Information Ratio, Q_LMIR.

Table 1. Registration metrics.

Metric	Test Sample	SIOF	CLS	CFM	CPC	Proposed
IR	1	143.217	161.19	174.546	154.066	187.561
	2	143.217	161.19	174.546	151.113	187.561
	3	115.606	129.585	140.471	124.631	152.262
	4	132.028	148.477	160.575	141.78	172.847
	5	124.796	140.42	151.94	134.198	163.321
	6	124.796	140.42	151.94	134.198	163.321
	7	144.765	140.56	148.76	138.91	155.54
	8	158.654	152.768	150.65	153.221	161.21
	9	132.343	137.434	135.121	134.461	141.54
	10	148.351	142.67	147.34	145.233	150.23
LIR	1	78.3473	90.0022	98.8611	85.623	107.59
	2	78.3473	90.0022	98.8611	84.549	107.59
	3	65.2882	74.8208	82.2869	71.4939	90.1629
	4	72.2711	83.0021	91.151	78.9119	99.4178
	5	68.9627	79.2008	87.0091	75.3464	94.7196
	6	68.9627	79.2008	87.0091	75.3464	94.7196
	7	79.768	79.451	81.331	83.328	98.889
	8	81.998	89.492	84.341	88.446	97.344
	9	89.431	88.402	87.318	88.985	97.987
	10	88.2241	87.331	89.734	90.351	96.435
MIR	1	143.145	161.096	174.572	46.9168	187.601
	2	48.3193	48.0257	47.0589	46.9168	45.8579
	3	83.3262	90.4525	94.8126	87.4357	97.2554
	4	57.9456	59.5488	59.2475	57.7371	58.7392
	5	58.0502	60.7047	62.3834	59.6138	64.5034
	6	56.8043	59.7738	60.5487	57.5558	62.2276
	7	70.564	69.982	72.854	73.442	76.989
	8	70.112	69.667	71.334	70.667	74.552
	9	72.665	73.4498	74.375	76.234	78.318
	10	61.098	61.132	62.003	63.224	65.739
LMIR	1	78.3473	90.0022	98.8604	161.002	107.59
	2	146.906	169.367	186.518	158.854	203.397
	3	103.573	118.964	130.929	113.414	144.14
	4	132.2	152.545	168.046	144.858	183.783
	5	124.397	143.197	157.52	136.095	171.625
	6	123.754	142.435	156.778	135.416	170.9
	7	148.897	147.628	151.546	149.981	195.887
	8	166.496	165.894	159.548	171.44	178.752
	9	180.568	181.564	179.552	168.519	189.562
	10	169.451	168.435	170.784	167.334	179.532

For quantitative analysis of the image registration methods, Information Ratio (Q_IR) and Mutual Information Ratio (Q_MIR) [54] are used. Here Q_IR and Q_MIR will be calculated including their lower bounds, i.e. Q_LIR and Q_LMIR are the functions of mutual information and entropy. These metrics are useful to minimize the computational complexity in identifying the feature count. These metrics will be useful to evaluate the registered image from the perspective of common features count that can be identified among given samples. The metrics for the ten test samples of medical images are given in Table 1.

To visualize the registration effect of the four existing and our proposed medical image registration techniques, we offered a graphical representation of objective assessment indices in Figure 9. In every metric, the proposed algorithm has obtained superior outcomes than the state-of-the-art approaches.

Figure 10 presents the average objective performance analysis of medical image registration methods. When compared to the four state-of-the-art approaches, our proposed method outperforms them all.

Figure 10. Average metric values of proposed and existing techniques.

Whether using subjective visual evaluation or objective indicators, the suggested strategy yields improved overall registration outcomes. Nevertheless, we may state that the proposed method produces satisfactory results.

5. Conclusion and future scope

Medical image registration with object deviation estimation through motion vectors using octave and level sampling has been proposed in this article. Many medical registration algorithms have the limitation of quantitative analysis while detecting abnormalities like loosening, subsidence and anteversion. The proposed algorithm will perform quantitative analysis in terms of Information Ratio (Q_IR), Mutual Information Ratio (Q_MIR) and their lower bounds. Our algorithm helps medical practitioners detect the amount of object deviation post-surgery using motion vector analysis. Both qualitative and quantitative assessments have found that our method outperforms the current state of art methods. In the proposed method, we concentrated on translational movement only. There is a possibility to identify the abnormalities even if it is affected by rotational and skew motion. This will be possible through 3-D image analysis.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

Any data and code used in this paper are available from the authors upon reasonable request.