A Deformable 3D-3D Registration Framework Using Discrete Periodic Spline Wavelet and Edge Position Difference

Neck pain is one of the most common symptoms of cervical spine disease, and segmenting neck muscles to create volumetric measurements may assist clinical diagnosis. While image registration is used to segment medical images, registration is highly challenging for neck muscles due to their tight proximity, shape and size variations among subjects, and similar appearance. These challenges cause conventional multi resolution-based registration methods to be trapped in local minima due to their low degree of freedom geometrical transforms. This article presents a novel object-constrained hierarchical registration framework for aligning inter-subject neck muscles. First, to handle large scale local minima, the proposed framework uses a coarse registration technique, which optimizes the new edge position difference (EPD) similarity measure, to align large mismatches. Also, a new transformation based on the discrete periodic spline wavelet (DPSW), affine and free-form-deformation (FFD) transformations are exploited. Second, to avoid monotonous nature of using transformations in multiple stages, a fine registration technique is designed for aligning small mismatches. This technique uses a double-pushing system by changing edges in the EPD and switching transformation resolutions. The EPD helps in both coarse and fine techniques to implement object-constrained registration via controlling edges, which is not possible when using traditional similarity measures. Experiments are performed on clinical 3D magnetic resonance imaging (MRI) scans of the neck, with the results showing that the EPD is more effective than the mutual information (MI) and sum of squared difference (SSD) measure in terms of volumetric dice similarity coefficient (DSC). Additionally, the proposed method is compared with the diffeomorphic Demons and SyN state-of-the-art approaches with ablation studies in inter-subject deformable registration. The proposed method achieves better accuracy, robustness and consistency than the reference methods, with an average volumetric DSC of 0.7029 compared to 0.6654 and 0.6606 for the Demons and SyN methods, respectively.


I. INTRODUCTION
Neck pain is a common musculoskeletal condition affecting the cervical spine. The estimated prevalence of significant episodes of neck pain is 40 to 70 per cent of an The associate editor coordinating the review of this manuscript and approving it for publication was Trivikram Rao Molugu. affected person's lifetime [1]. Further, neck pain results in significant socio-economic costs related to health care expenses, work absences, reduced productivity and insurance claims [2]- [4]. Patients may present with either acute pain, particularly as a result of trauma such as that experienced in a motor vehicle accident, or more chronic pain. In some instances of neck pain, the cervical muscles, in magnetic resonance (MR) images, appear to undergo pseudo hypertrophy due to fat infiltration or atrophy. However, these changes are inconsistent between the muscles and vertebral levels [5]- [7]. These inconsistencies may be caused by the measurement method rather than the marker [5]. First, quantification of the sizes of muscle and fat infiltrations using cross-sectional measures from two-dimensional (2D) images at selected vertebral levels is not representative of 3D muscle measurements [8]. Second, when using 2D quantification methods, the data may suffer from errors associated with the effects of a partial volume or alignment of the axial slice relative to the cervical spine [9]. Therefore, volumetric studies of neck muscles are required before the use of their sizes and fat infiltrations as reliable markers can be verified or refuted.
3D-3D cervical MRI registration can provide volumetric muscle segmentation for neck pain treatments. However, registration between the neck images of two individuals is difficult due to the potential of large anatomical variability (see Figure 1). Within the small and narrow region of the neck, the 27 muscles that control the movements of the cervical spine are compactly arranged (see Figures 1,5,6 and 8) with approximately similar composition resulting in similar intensity. This intricate relationships make them difficult to identify accurately. Further, the variable quantities of fat present within (see Figure 1) and between the muscles make their boundaries challenging to demarcate. This may be further complicated by unconscious movements of the patient (e.g., swallowing and breathing) during the MR scanning process, which can deteriorate the image quality. Since, to the best of our knowledge, there is no registration algorithm that has been developed for 3D-3D registration between neck volumes, it is essential to develop such an algorithm. Registration methods generally attempt to find a global minimum and avoid local ones, but cannot avoid all local minimums due to the many that are classified as small-scale dips and large-scale basins [10]. Escaping from basins is vital as they yield large mis-registrations and usually occur in large initial mismatch scenarios. Conversely, dips yield small mis-registrations occurring in small mismatch scenarios and occur more often. Neck medical data contain both scenarios discussed above. Conventional multi resolution-based registration methods [11]- [14] cannot avoid basins and dips completely due to the lower available degree of freedom of geometrical transformation and the monotonous types of transformations in multiple stages. Information theoretic-based methods have been widely used for 3D-3D medical image registration [15]- [17]. However, global information theoretic-based methods are less sensitive to local deformation and often encounter the mis-correspondences problem due to the lack of distinctiveness of the similarity measure (SM). Due to the compactness and similar appearance of neck muscles, the problem of mis-correspondences is more acute in neck data. Further, local information theoretic-based methods are computationally expensive and face statistical instability problems. In addition, information theoretic-based methods ignore anatomical information, which is crucial for guiding correspondence detection and registration. The feature-based registration methods establish correspondence through high-order anatomical information [18]. However, some features are often partially invariant in images with a different appearance, which is common in clinical applications. Some descriptor-type features show distinctive characteristics; nevertheless, its performance faces challenges in scenarios involving large anatomical variation.
Our proposed work develops a coarse-to-fine 3D-3D registration approach for dealing with the neck's high level of variability. It uses a hybrid registration framework divided into coarse and fine mismatch correction sections to handle the basins and dips, respectively, and exploits the diffeomorphic Demons algorithm in its last stage to boost alignment. The coarse and fine sections reduce the chances of becoming stuck in basins and dips, respectively, the latter through its double-pushing system. This method makes the following four key contributions: • In our multistage framework, multiple 3D transformations at the coarse level are used in multiple stages for large mismatches, which are usually of different types. Therefore, they can push basins from different directions with greater angle differences, as opposed to similar transformations in conventional multi resolution methods, thereby reducing the chance of becoming stuck in basins.
• A novel discrete periodic spline wavelet (DPSW)-based 3D transformation, which requires fewer parameters than the free-form-deformation (FFD) one but has similar benefits, is developed in the coarse section. It reduces the burden of optimization and provides variations at the global level with no stretching or shrinking effect, unlike a discrete cosine (DC), which is popular in video coding.
• The framework incorporates a new SM called edge position difference (EPD), which uses our modified 3D Chamfer distance transform algorithm. Since it uses the edges of an object, it provides an opportunity to tune through multiple stages using different sets of edges. VOLUME 8, 2020 It aligns the neck's trunk using the same set of strong edges in the multistage coarse section and, gradually, incorporates the weaker ones for the muscles and other small objects in the fine section. It can do object-wise alignment through multiple transformations, whereas traditional multi resolution methods use SM methods that are incapable of object-wise alignment.
• A double-pushing system is designed for the fine section to reduce the chance of becoming stuck in dips, which occur more often than basins, whereas a single-pushing system is used in most multi resolution methods. The double-pushing system yields small deformations formed through changing the number of edges of the EPD and the resolution levels of the transformation, rather than changing only the latter, as in traditional multi resolution techniques. The EPD, mutual information (MI) and sum of squared difference (SSD) SMs are compared using the affine transformation and the EPD achieves good accuracy for our clinical dataset. Additionally, our proposed method is compared with the diffeomorphic Demons [19] and SyN [20] algorithms, which are state-of-the-art registration approaches, and outperforms various other non-rigid registration algorithms [21]. We compute the volumetric dice similarity coefficient (DSC) in our real clinical 3D MRI dataset using the proposed, diffeomorphic Demons and SyN methods, with the proposed method achieving a substantial improvement in accuracy.
The rest of this article is organized as follows. Related work is discussed in Section II, we describe details of our deformable 3D-3D registration method in Section III, Section IV presents the experimental procedure and results, and Section V and Section VI provide a discussion and conclusion, respectively.

II. RELATED WORK
Image registration [22]- [24] is a basic image-processing technique whereby two or more images are aligned by keeping one stationary (called a fixed image) and moving another (called a moving image) towards it [18]. It comprises a geometrical transformation, similarity measure and optimization, with 3D-2D registration currently being developed commercially and rigid registration [25] practically available [26]. Inter-subject registration (ISR) in 3D, a kind of deformable image registration [27], is a key challenge due to anatomical variability [26] preventing uniformity. Although much work on 3D-3D ISR has been conducted in the last two decades, mainly on brain images, its accuracy is not clinically acceptable, with specialists considering that further research must yet be undertaken [26], [28]. Moreover, 3D-3D ISR often faces some problems compared to 2D-2D registration. First, the optimization becomes difficult as more parameters are required for geometrical transformation. Second, more dips and basins are required to handle higher-dimensional space. Third, the computational cost becomes expensive. Finally, the SM faces discontinuity in intensity problems along out of plane, since most medical imaging modalities keep spacing between the slices. Most studies on 3D-3D ISR have focused on optimization [13], [21], [29], local regularization [30], [31], multi resolution FFD [11], [12] and the application of the diffeomorphic log-demons algorithm [32], [33] to 3D-3D ISR. The most recently developed competitive registration algorithms are the multi resolution FFD and diffeomorphic log-demons approaches. Hua et al. [11] proposed a 3D-3D deformable registration for handling discontinuities by adding extra degrees of freedom to a multi resolution framework using a parameter up-sampling method that required segmenting a target image to determine discontinuities and allowing more time to optimize additional parameters. Sun et al. [34] proposed a random perturbation technique for a multi resolution nonlinear registration framework for 3D-3D and 2D-2D applications using a lower-order B-spline, retaining the same smoothness as a higher-order spline to reduce the execution time. However, in terms of accuracy, it could not perform well in other clinical applications. Sun et al. [35] proposed another 3D-3D simultaneous multi resolution strategy in which different resolutions of the spline and data were used together to improve performance, but their method was dependent on a parameter whose value varies for different applications.
Overall, conventional multi resolution methods use only the coarse resolution levels of a transformation to resolve large mismatches. These are not sufficient to achieve proper correction, since, when the optimization process is stuck in a basin, to escape easily, it must be pushed in a specific direction depending on the particular basin. However, consecutive resolution levels in multiple stages of multi resolution methods have almost the same characteristics with different directions and small variations in angle but are considered different transformations. Therefore, consecutive coarse levels may not thrust in the required direction for some basins because of their small angle variations. Large angle variations among multiple transformations in multiple stages may be required so that a transformation can push basins in the appropriate directions. The optimization process will not face any difficulty to determine the optimal direction because the multiple transformations are not applied simultaneously rather they are applied consecutively as each stage of the algorithm is performed. Therefore, as traditional multi resolution methods cannot push in the required direction for some basins, they are unable to avoid all the basins. Further, these methods try to tackle small mismatches, which cause dips more frequently, by changing only the fine resolution levels of a transformation in multiple stages with the same SM and optimization, and are thus incapable of eliminating all of them.
Feature-based registration methods have great effects in 3D-3D ISR due to their use of anatomical information, which helps to find correspondence detection effectively. There are many types of features in the registration literature, such as histogram of oriented gradient (HOG) [36], gradient location and orientation histograms (GLOH) [37], scale-invariant feature transform (SIFT) [38] and speeded up robust features (SURF) [39]. The gradient information-based features are partially invariant [18], which causes misregistration. The SIFT and SURF methods require the same features to be detectable in both fixed and moving images, which is not possible in neck MRI data as the muscles are very compact and images are obtained with different acquisition protocols. In fact, SIFT and SURF are more suitable for natural image analysis than medical image analysis. The RANdom SAmple Consensus (RANSAC) algorithm is often used in conjunction with feature-based registration methods to filter out excellent matches. For example, Kahaki et al. [40] proposed a local intensity maxima feature-based registration method for in vivo time lapse microscopy images. They used a two-step feature matching procedure in which features are initially matched coarsely and then the matching features are refined through RANSAC. The iterative closest point (ICP) algorithm is a popular method in shape registration [41]- [43], which has high accuracy for point set registration. A 3D Canny edge-based objective function is used in medical image registration for pose estimation and shape reconstruction [44], multi-modal geometric matching [45] and respiratory motion correction [46]. A self-similaritiesbased feature called a modality-independent neighborhood descriptor (MIND) [47] was proposed to provide distinctive correspondences in objective function by incorporating neighborhood pixels' information. It showed better results in cases of similar local structural patterns in small regions than other SMs. However, its performance can be challenged in cases of large local anatomical variation. Further, it cannot hide the influence of contrast enhancement and embeds unwanted information.

III. DEFORMABLE 3D-3D REGISTRATION MODEL A. OVERVIEW
A diagram of the operational flow of the proposed hybrid registration framework is displayed in Figure 2. As a pre-processing step, all the original fixed and moving MRI volumes are trimmed and interpolated to volumes of 128 × 128 × 128 voxels, since they contain some unwanted information. The volume of interest is selected as the volume between the C1 and C7 vertebral levels, which represent the top and bottom vertebrae of the neck and are the landmarks most commonly used to assess muscle morphometry [7]. The processed volumes are then manually delineated to obtain the ground truths, as discussed in detail in Section IV.
The framework is divided into coarse and fine sections. All the geometrical transformations at the different stages in Figure 2 are separable, since each stage uses the registered moving volume of the prior stage as its moving volume and the same fixed volume used by the prior stage as its fixed volume. Thus, the transformations are obtained separately from the framework. To describe our experimental results in Section IV, we have assigned the following labels for different stages of the algorithm: Affine-EPD, DPSW-EPD, Coarse-EPD and Fine-EPD.

1) COARSE SECTION
All the stages in this section use the same strong edges of the MRI volumes to align the neck's trunk and the boundaries of other large objects. In each stage, the geometrical transformation is changed to avoid basins and, importantly, obtain a good alignment. Other elements, such as the SM or optimization method, could also be changed in each stage. A local minimum is considered for a specific combination of the transformation, SM and optimization. We use the affine, the DPSW and the coarsest level of the FFD as transformations to combat basins. Therefore, multiple transformations can attack the optimization from different directions using large angle variations to pull out from basins.

2) FINE SECTION
In this section, five stages are used to obtain fine deformations. This section helps to reduce the chance of the optimization algorithm converging to a dip. As dips occur more frequently than basins and could cause the optimization to be stuck at any stage, we have designed a double pushing system to combat dips. The double pushing system is implemented by changing the transformation and SM simultaneously at every stage. The transformation change is performed by using different levels of the spline in the FFD. The SM change is accomplished by using the attributes of the EPD which allow different sets of edges to indicate different ranges of values. Different sets of edges for specific volume pairs are used for different stages. These sets of edges are changed gradually from strong to weak. The strong edges are a subset of the set of weaker edges. The first four stages in this section use coarse to fine levels of the spline for the FFD with corresponding different sets of strong to weak edges, respectively. The fine deformations are achieved by using the VOLUME 8, 2020 weak sets of edges and the fine levels of the spline. The strong sets of edges and the coarse levels of the spline are used to correct coarse deformations. The gradual change protects the framework against mis-correspondences. Finally, we apply the Demons algorithm [19], [48] in the fifth stage of the fine section to correct more fine mismatches. Actually, the stages before Demons bring the moving volume closer to the fixed volume which helps Demons to align more effectively than when only using the Demons method.

B. GEOMETRICAL TRANSFORMATION
The choice of transformation has a large effect on the registration process [49], with the most appropriate one not known as a prior [50]. In the registration process, the transformation parameters are estimated using an optimization technique, with the number of them referring to the deformation's degrees of freedom. In our application, the registration needs to be performed between the MRI volumes of two different individuals' necks; this cannot be achieved using only an affine or rigid transformation because both have a limited number of parameters. Therefore, an elastic transformation with a higher degree of freedom is required to tackle the morphological complexity and variability of the population. However, to deal with the neck's variability, we use a mixture of affine, DPSW and FFD transformations to align the neck's trunk first by exploiting the advantages of the EPD. The same strong edges that correspond mainly to the neck's trunk are used to correct a coarse mismatch with a different transformation. This is because a good deformation cannot be achieved through a single transformation.
In this study, considering F(x, y, z) and I (x , y , z ) as fixed and moving volumes, respectively, their coordinates are involved in elastic registration as follows: where m k are the motion parameters, k is the parameter index, P is the total number of motion parameters and ϕ k is the basis functions for the complex mapping, given as: There are many types of basis functions in the literature, including polynomial, Fourier, radial, B-spline, DC and wavelet. Of these, the Fourier, B-spline and wavelet functions support a multi resolution decomposition that provides a coarse-to-fine representation of the displacement field. Hence, these basis functions are usually used in medical image registration. However, wavelet basis functions can achieve a local deformation more effectively than Fourier basis functions due to its localization in both the frequency and spatial domains [51]. The DPSW and FFD transformations will be described in Section III-B1 and Section III-B2.

1) DISCRETE PERIODIC SPLINE WAVELET
We use the Cai-Wang [52] wavelet, which is compactly supported, in our DPSW-based basis functions and a fourth-order B-spline as a scaling function. The wavelet is arranged in a periodical form with the center of the main lobe translated to a coordinate origin.
The fourth-order B-spline is defined as: where, for any integer (n): Then, the spline wavelet is: The DPSW-based basis functions are: where k = 2su + sv + w + 1, u, v, w = 0, 1, 2, · · · , s − 1, s = 3 P 3 and ψ λ () is the DPSW. The supports are 4 and 3 for the spline and spline wavelet respectively, as shown in Figure 3, with the periodic spline wavelet generated by considering the wavelet as one period. Our 3D basis functions are generated from the fifth resolution of the wavelet using a point-to-point multiplication of the 1D functions, whereas the FFD uses a tensor product of the 1D non-periodic third-order spline to generate basis functions. The DPSW-based transformation requires fewer parameters than the spline or B-spline-based wavelets to represent local deformations. Specifically, one basis function in the DPSW can represent local deformations over an entire image, whereas the B-spline or spline wavelet cannot, since they require more parameters to represent the same level of local deformation for a whole image, which places a burden on the optimization process. In particular, for a 128 × 128 × 128 image, a B-spline-based transformation requires 15, 27 and 51 parameters for the fifth, fourth and third resolutions, respectively, while the DPSW uses only 24 parameters for all resolutions. Although DC-based basis functions can obtain local deformations over an entire image and are widely popular in video and various image-processing applications [53], they cause shrinking and stretching in several parts of the image. There are two reasons for these effects. First, the parameterizations in the available DC-based basis functions do not use the full cycle of a cosine wave within a cubical image support, whereas we use multiple cycles of the periodical spline wavelet in the DPSW-based basis functions, which have greater variations in values, as shown in Figure 4. We also perform registrations using the DPSW and DC on the 3D MRI volumes shown in Figure 5 and obtain stretching effects in the latter's results. Second, there is a lower span in the negative lobe of the spline wavelet with regularization when compared to the cosine one.

2) FREE FORM DEFORMATION
We use the fourth-order B-spline defined in (3) in the FFD-based basis functions to obtain smoother local deformations than those in the traditional FFD using the third-order B-spline, with the basis functions: where t is the translation index, j is the resolution level, R is the starting resolution level, L is the maximum resolution level and () is the basis functions. As L depends on the volume's size, the number of parameters for a resolution level is (2 −j N + 1) × 3 if N = N x = N y = N z and the 3D basis functions use the tensor product of the 1D spline, given by:

C. EPD SIMILARITY MEASURE
As the SM is another significant part of the registration process, in our method, a new measure called the EPD is leveraged in our coarse-to-fine registration framework to accomplish alignments in both the coarse and fine levels. This measure is based on the hierarchical Chamfer matching algorithm [54] and determines the distance between the corresponding edges in two images. The EPD is not exactly the same as root mean square (RMS), Euclidean distance or closest distance which are commonly used in shape matching. First, the EPD uses the arithmetic mean of Chamfer values of the moving image at the position of edges in the fixed image. Second, the Chamfer distance uses an approximation of the Euclidean distance to calculate distance from a pixel to the nearest edge. Third, the edge points of the fixed image may map to Chamfer values which may correspond to different edges other than the edge in the fixed image. One of the most popular intensity-based SMs in the literature is the MI measure. However, it is not suitable for neck MRI datasets in which multiple muscles are near each other and have similar compositions, with large deformations between subjects. Therefore, using a MI-based SM causes mis-registrations between MRI volumes. Conversely, the EPD is a feature-based SM that uses the edges of the muscles and the neck's trunk. In this registration technique, every edge pixel in a moving image contributes to a registration error with a value proportional to the distance to its closest edge in the reference image. Therefore, the edges of overlapping images are attracted to each other, which leads the EPD SM to have better registration accuracy than the MI measure, as justified in Section IV-C. Moreover, the EPD supports coarse-to-fine tuning because it uses edges that are VOLUME 8, 2020 controllable by selecting their detection thresholds. Consequently, since it is more suitable than the MI, this new SM is used in our proposed framework.
To calculate the EPD, first, edge volumes of the moving I and fixed F volumes (E I and E F , respectively) are calculated using a canny edge detector, with a value of 1 corresponding to an edge; otherwise, the value is 0. We can select strong or weak edges by choosing thresholds in the canny edge detector, which assists in designing a double-pushing system for combating dips in the fine section of the registration framework. Then, the locations of the edges (β) in E F are calculated and the distance transform volume (C) from E I is determined using the 3D Chamfer distance transform method. This transform is an approximation of the Euclidean distance transform, in which each value of a voxel represents a distance to the closest edge in the edge volume. Finally, the EPD S, which is the arithmetic mean of the values of the voxels in the β positions, is calculated as: where T is the total number of edge voxels in β, 3 is used to compensate local distance in the Chamfer distance. Figure 6 presents an example of a MRI slice with its corresponding edge image and the latter's Chamfer distance image, obtained from a 3D MRI volume using the 3D Chamfer distance algorithm with local distances of d 1 = 3, d 2 = 4 and d 3 = 5 [55]. However, if there are no edges in some consecutive slices of a 3D volume or some local 3D regions of it, infinity remains in the border voxels of the cubical support. Therefore, we modify the original 3D Chamfer distance transform algorithm, as described below.

1) MODIFIED THREE-DIMENSIONAL CHAMFER DISTANCE TRANSFORM ALGORITHM
To tackle the problem mentioned above regarding the boundary voxels, we use two forward and two backward passes rather than only two passes as in the original algorithm. The new algorithm handles the other voxels similarly to those in the original algorithm. We process the boundary voxels of the last and first slices separately in the first forward and backward passes, respectively. The masks of the original algorithm are changed according to the positions of the boundary voxels and are used in our algorithm, with Figure 7 (a) showing their modified positions and Figure 7 (b) showing the modified mask for special condition 3. There are five special conditions in the first forward and backward passes, and four in the second ones. The framework is summarized in Algorithm 1.

D. TRANSFORMATION PARAMETERS OPTIMIZATION
Optimizing a registration is considered ill posed and is actually a multidimensional problem that maximizes or minimizes the SM with respect to the transformation parameters. Generally, optimization methods for medical image registration are classified into three categories: continuous, discrete and miscellaneous. The first category includes GD, conjugate gradient, Powell's conjugate directions, quasi-Newton, Gauss-Newton (GN), Levenberg-Marquardt and stochastic GD approaches, the second includes graph-based, belief propagation and linear programming techniques and the third includes greedy and evolutionary algorithms. The Powell's conjugate directions and stochastic GD methods have been applied for transformations with low degrees of freedom [49], while evolutionary algorithms, which are used mainly in linear registration, have shown slow convergence. Hence, GN and GD approaches have been used in many medical image registration approaches. In our application of mono modal registration, the GNGD method is used because the EPD involves summing the function values.
The optimization procedure required to find m k is described as follows: for q = 1 : N do 5: Update C by Q where the value of E I is zero, and by zero otherwise. 6: end for 7: end for 8: end for 9: for r = 2 : O do 10: for p = 2 : M do 11: for q = 2 : N do 12: if p < M ∧ q < N then 13: Update C according to forward mask in Figure 17 in [55]. 14: else if r = O then 15: for p = 1 : M do 16: for q = 1 : N do 17: if p = 1 ∧ q = 1 then 18: Update C according to the first special condition in Figure 7. 19: else if p = 1 ∧ 1 < q < N then 20: Update C according to the second special condition in Figure 7. 21: else if p = 1 ∧ q = N then 22: Update C according to the third special condition in Figure 7. Update C as first forward pass using backward mask in Figure 17 in [55] and special conditions 6, 7, 8, 9 and 10 in Figure 7 for the first slice. 38: end for 39: end for 40: end for 41: Return C To minimize S, it is necessary to estimate its values in a small neighborhood of m k , as: where m is a vector added to m k in each iteration. The value of m is determined using the first-order Taylor series approximation, as: The ∇S(m k ) is calculated as: where ∂C ∂x , ∂C ∂y and ∂C ∂z are the spatial gradients of C, while ∂x ∂m k , ∂y ∂m k and ∂z ∂m k are obtained from corresponding transformation basis functions (ϕ k ) in equations (6) and (7).
The Hessian matrix ∇ 2 S(m k ) in equation (12) is defined in the GN optimization technique as: where J (m k ) is the Jacobian of C(x , y , z ), calculated as: where b is the parameter indices and c is the position index of β. Finally, we obtain the vector m from equation (12), which is used to update the transformation parameters as follows: As most of the elements in the fine resolution's basis functions are zero in the FFD transformation, most of those in the Hessian matrix become zero, since this matrix is calculated by multiplying the spatial gradients of C and basis functions in the β positions, making it impossible to invert the matrix. Therefore, we omit this matrix in the FFD transformation, which changes the optimization to a GD method.

IV. EXPERIMENTAL PROCEDURE AND RESULT ANALYSIS A. DATA AND ANNOTATIONS
Our experiments were performed on neck MR images as part of a larger study of neck pain undertaken at the Australian National University (ANU) Medical School. The study had ethics approval from the Human Research Ethics Committees of the ANU and Australian Capital Territory Health, with written informed consent from all participants. T1 SE axial spin echo MR images with voxel spacings of 0.8594 mm × 0.8594 mm×4 mm and sizes of 256×256×45 were acquired using a 3 Tesla Skyra scanner (Siemens, Erlangen, Germany) VOLUME 8, 2020 and then cropped and interpolated to 128 × 128 × 128. Data from 11 participants collected at different times using different MRI machine settings were used in our experiments. The participants' demographics and the MR sequence parameters for each scan are shown in Table 1. An ANU graduate-entry medical student with a degree in anatomy manually delineated the right and left sternocleidomastoid, semispinalis capitis and splenius capitis muscles between the C1 and C7 vertebral levels. These segmentations were validated and edited by two medical experts from the ANU Medical School and Canberra Hospital. A MATLAB graphical user interface developed by our team for segmentation could use as many vertices as necessary to capture small details of the contours of these muscles. The contours of a slice from an MRI volume are shown in Figure 8. Different colors are used for each separate muscle due to annotation convenience; however, the same color is used for symmetric muscles.

B. EVALUATION METRIC
To evaluate our registration results, we used the DSC due to its popularity in medical image research. The DSC is expressed as: where g is the annotation contour and s is the transformed contour, with the volumetric DSC for each muscle calculated separately.

C. NUMERICAL RESULTS ANALYSIS
We performed experiments on a HP z230 tower workstation with a 16 GB RAM and 3.40 GHz Intel(R) Core(TM) i7-4770 processor running the Windows 7 operating system using MATLAB, and C and C++ MEX programming. A canny edge detector with sigma 1.5 was used to calculate the edge image for the EPD similarity measure. The lower and higher thresholds pairs 0.1, 0.9; 0.08, 0.7; 0.04, 0.4; 0.01, 0.2 and 0.001, 0.1 were selected in the 5 th , 4 th , 3 rd , 2 nd and 1 st resolution levels of the FFD registration respectively. The threshold pair 0.1, 0.9 was also used with affine and DPSW registrations.

1) COMPARISON OF EPD, MI AND SSD SIMILARITY MEASURES
We performed affine EPD, affine MI and affine SSD registrations on our dataset to compare the EPD, MI and SSD similarity measures. The affine MI registrations were performed using advanced normalization tools (ANTs). For each fixed patient, the other 10 patients were considered moving, giving a total of 11 × 10 = 110 cases (registrations). Table 2 shows the experimental results obtained for affine registrations-that is, the DSC values for 10 of the 110 cases for the six muscles in which PT-1's volume was considered fixed and the others moving. Similar tables were constructed for the other patients. The bold values, which represent the highest DSC in each column, indicate that the EPD has the highest DSC in eight cases and the MI and SSD each in only one. It is also clear that the PT-2 case has the lowest best value, with best DSCs for the EPD, MI and SSD of 0.8179, 0.7845, and 0.6735 respectively. Overall, the EPD achieved better accuracy than the MI and SSD. Table 3 shows the mean and overall mean DSCs for all the fixed volumes of the affine registration. The last column represents the final means of the 110 cases and the last three rows display the means of 10 cases when considering the other 10 patients as moving. The table indicates that the EPD has better registration accuracy than the MI and SSD for the nine patients, and SSD has better accuracy than EPD and MI for the two patients, EPD achieving a final DSC of 0.5477, significantly better than the MI's 0.2609 and the SSD's 0.5188. This large difference in DSC values proves that the MI similarity measure is not suitable for our dataset and, therefore, for our coarse-to-fine framework. Although the accuracy of SSD is close to that of EPD, it can not provide object-wise alignment. Therefore, we have selected EPD as the similarity measure for our proposed method. Although the volumetric DSCs for EPD are better than for the other methods, its values are not in the excellent range because strong edges were used to align the volumes coarsely. The coarse alignment is further refined through the latter multiple stages in our proposed method. These results are obtained TABLE 2. Results obtained for the edge proposition difference (EPD), mutual information (MI) and sum of squared difference (SSD) measures from the affine registration experiments in terms of volumetric dice similarity coefficient (DSC) (DSC values for patient-1 [PT-1] fixed and others considered moving; left and right sternocleidomastoid, left and right semispinalis capitis, and left and right splenius capitis muscles denoted as Muscles 1, 2, 3, 4, 5 and 6, respectively, where higher values indicate better performances, with maximum value 1, and best DSC in each column marked in bold). using only the single stage "Affine-EPD" from our proposed framework. The EPD yields mediocre results in some circumstances such as PT-11. There could be two possible reasons. First, the patient is demographically and anatomically more different than the other patients. Second, there may be some noise or intensity inhomogeneity problems. These problems can be eliminated by using our proposed complete framework.
To investigate further, we used a box plot to statistically analyze our experiments. Figure 9 exhibits the registration accuracies obtained by the EPD, MI and SSD for each muscle separately (see Figure 9(a)) and combined (see Figure 9(b)). Figure 9(a) indicates that the EPD has better median DSCs than the MI and SSD for four muscles, and the SSD is better for the other two muscles. In addition, the maximum DSCs for EPD are higher than for the MI and SSD measures for all muscles except for Muscle-5 when using MI. Similarly, Figure 9(b) reveals that the EPD has better overall median and maximum DSCs for all the muscles combined (0.5866 and 0.8627, respectively) than the MI (0.2430 and 0.8441, respectively) and SSD (0.5484 and 0.8220, respectively).

2) PERFORMANCE COMPARISON
We used strong and weak sets of edges to obtain coarse and fine level deformation. Figure 10 shows edge maps and the corresponding deformation fields for an inter-subject registration case using strong and weak sets of edges. These maps  illustrate a coarse displacement field for the strong set and a fine displacement field for the weak set. To evaluate the effectiveness of every stage of our proposed method, we calculated the volumetric DSC for one of our inter-subject registrations at every stage. Figure 11 shows the 2D visual results at different stages of our proposed registration method with a corresponding volumetric DSC value. It can be seen from Figure 11(c) that the required deformations are very large across subjects. Every stage gradually aligns the two images as indicated by the alignment of the vertebra, muscles and neck trunk with the superimposed edges. The effectiveness of some stages is more clearly visible in some other inter-subject registration cases. The rising trend of volumetric DSC proves the effectiveness of every stage in our proposed framework. Figure 12 shows the corresponding coronal views for the patient shown in Figure 11. We have not shown the sagittal views due to anterior information loss caused by inhomogeneity problems with the MRI scanner. It exhibits similar alignment improvement with every stage as in the axial views. Figure 13 shows the results in terms of muscles contours of different methods compared with the ground truths. It can be seen that the proposed method's contours are more fairly matched than for the SyN and Demons algorithms.
We compared our full proposed method for the 110 inter-subject cases with Coarse-EPD and Fine-EPD using volumetric DSC. We excluded other stages for presentation   [20] and Demons [19] registrations in our neck dataset. The SyN and Demons registrations were conducted through ANTs and MATLAB, respectively. Multistage registrations were used in the SyN method, in which the rigid, affine and SyN geometrical transformations used the same smoothing sigmas, shrink factors, and convergence, except that SyN used different convergence. The metric CC means cross-correlation.
convenience. We also performed ISR using the diffeomorphic Demons [19], [48] and SyN [20] methods on our dataset for the 110 cases. The two methods are powerful and popular registration algorithms considered as a gold standard in the deformable registration field. We used ANTs and MATLAB for implementing the SyN and diffeomorphic Demons, respectively. The parameter information for the SyN and diffeomorphic Demons is shown in Table 4. Table 5 shows the mean and overall mean volumetric DSCs for all the fixed volumes of the proposed, diffeomorphic Demons [19], SyN [20], Coarse-EPD and Fine-EPD registrations. Tables 3 and 5 were generated in a similar way. The overall mean values in the last column represent the final means for our 110 cases and those in the last five rows represent the means for the 10 cases considering the other 10 patients as moving. The proposed method clearly yields better registration accuracy than the others for nine of the 11 participants, whereas Demons [19] and SyN [20] yield good accuracies for only one each of the 11 participants. The proposed method achieves a final DSC of 0.7029, significantly better than those of the Demons (0.6654) and SyN (0.6606), with the overall mean volumetric DSC values obtained for five muscles higher by the proposed method than others, except for Muscle 2 by SyN. However, the rising trend of final DSC of Coarse-EPD and Fine-EPD indicate that the stages before Demons stage in our proposed method bring the moving volumes closer to the fixed volume which helps Demons to align more effectively than the only Demons method. Figure 14 shows the results obtained from the statistical analysis of the proposed, diffeomorphic Demons, SyN, Coarse-EPD and Fine-EPD experiments for each muscle separately (see Figure 14(a)) and combined (see Figure 14(b)). Figure 14(a) indicates that the proposed method has better median and maximum DSCs than the others for Muscles 1, 4, 5 and 6, whereas the SyN algorithm performs better for Muscles 2 and 3.
Similarly, Figure 14(b) reveals that the proposed method has better overall median and maximum volumetric DSCs for all the muscles combined (0.7385 and 0.9075, respectively) than the diffeomorphic Demons (0.7215 and 0.8680, respectively) and SyN (0.7137 and 0.8940, respectively) methods. Moreover, the diffeomorphic Demons method has more outliers (53) than the proposed method (20) and SyN (32) for all the muscles combined.
To perform a complete analysis, we also used the Hausdorff distance (HD) as a distance error metric. Figure 15 shows the HD results for all muscles combined. It reveals that the proposed method has better median HD (5.7446 mm) than the SyN (6.0000 mm) and diffeomorphic Demons (5.9161 mm) methods.
Finally, the experimental results and analyses reveal that the proposed method outperforms the diffeomorphic Demons and SyN algorithms in terms of registration accuracy and consistency.

V. DISCUSSION
A deformable 3D-3D fully automatic registration framework using a novel DPSW transformation and modified 3D Chamfer distance transform of the EPD was developed. The empirical outcomes demonstrated that our method outperformed the well-established diffeomorphic Demons [19], [48] and SyN [20] algorithms in terms of accuracy and consistency. Also, interestingly, the proposed framework was robust in terms of input volumes because it worked on multiple MRI scanner settings. In particular, the same thresholds were used in the EPD SM for the input volumes with different repetition and echo times, which is common when acquiring MR images in clinical practice.
The computational time required to register two neck volumes was calculated for the EPD, SSD, MI, proposed, Demons [19] and SyN [20] methods. These times are shown in Table 6. Although, some parts of the proposed method were implemented using C/C++ MEX coding, the computational time of the proposed method could be reduced if the full code was implemented in C/C++. The proposed method is slower than the Demons algorithm but has better accuracy. The low computational time for Demons was due to the professional implementation of this algorithm in MATLAB.
The experiments in our work were conducted on real clinical neck MRI data. Our proposed registration framework was thoroughly evaluated against well-recognized deformable registration methods [19], [20], [48] applied to our dataset, and obtained a better overall mean DSC value, confirming its effectiveness. The reasons for the consistently enhanced accuracy of our method are as follows:  • A conventional multi resolution registration method attempts to avoid large-scale basins in the function to be optimized, using only the coarse resolution levels of the spline. This approach has limited success because of the monotonous SM and optimization. In our work, the framework was constructed based on the novel notion that correct deformations cannot be obtained through a single transformation due to the constraints placed on the parameters of the geometrical transformation. This approach brings the floating volume nearer to the reference one by applying multiple transformations, optimizations and SMs. However, if an exact deformation is near a transformation, it cannot be obtained if the same stage is used repeatedly because it becomes trapped in a local optimum of the SM, to escape from which a different SM and optimization approach is required. In the past, FFD-based registrations attempted to escape large-scale basins by exploiting multi resolution versions of the B-spline with the same SM and optimization. However, they did not succeed due to the similar characteristics of the different B-spline resolutions. In this work, completely different combinations were used in each stage of the coarse section to align the neck's trunk first, which avoids large-scale basins in the optimization function and provides good alignment. • A registration process can frequently become stuck due to dips in the optimization function, rather than basins [14]. Although all the multi resolution techniques using the B-spline have tried to avoid these, they could not mitigate all their effects because they considered only variations of the spline resolution and the nature of their cost functions was non-convex [11], [34], [35]. For our framework, a double-pushing system was designed in the fine section to prevent dips yielding small deformations and to reduce the chance of becoming stuck in dips. The first process in this system involved changing the edges of the EPD from the strongest to weakest, while the second involved altering the resolution of the spline. Therefore, our framework obtained better accuracy by avoiding more dips.
• The negative lobe of the spline wavelet provides an implicit regularization in the transformation and, therefore, a smooth deformation. The DPSW has fewer parameters than the FFD, which facilitates the optimization process to avoid large-scale basins and improves matching.
• The EPD provides an opportunity to achieve an alignment through multiple stages for the same objects or edges. Therefore, tuning an alignment can be performed object-wise, a capability not available using any other SM. Although a limitation of our work is the lack of explicit regularization, which normally introduces a parameter that imposes a burden on the user, in our framework, we preferred the DPSW, which provides an implicit regularization and obtains smooth transformations through the fourth-order spline.
In our future work, the proposed framework will be validated using some public brain and lung datasets with a more theoretical insight into DPSW. Also, some other optimization methods, such as stochastic ones, will be explored, and other evaluation metrics for the proposed framework will be   [19] and SyN [20] algorithms on our computer for a single case. The affine transformation is used for the EPD, SSD and MI registration. The MI and SyN [20] are implemented through advanced normalization tools (ANTs).
investigated. Further, the neck dataset will be increased for experimental analysis.

VI. CONCLUSION
In this research, we proposed a registration framework with different combinations of transformations, optimizations and SMs with a new transformation DPSW and a modified 3D Chamfer distance transform algorithm. The DPSW requires fewer parameters than the FFD to represent similar local deformations, which decreases the burden on the optimization process. Our framework enhances the probability of finding the global minimum by reducing the effects of basins and dips, with outcomes obtained from experiments on real clinical datasets confirming its effectiveness.