Divide-and-Conquer Higher-Precision Thoracic CT Registration Achieved using Varying Degrees of Welsch Sliding Penalty: Thoracic CT Registration Using Varying Sliding Penalty

The registration of thoracic CT images poses challenges due to varying degrees of sliding motion in different regions. The commonly used l₂-norm displacement field regularization penalizes sliding motion heavily, which imposes strict constraints on sliding at areas like the lung surface. Previous studies have addressed this sensitivity to sliding motion by regularizing the displacement field with lₚ-norm (p ≤ 2), allowing sliding motion at regions such as the lung surface. However, these methods are global and can compromise the accuracy of rigid motion regions, and make optimization difficult due to non-smoothness. In this work, we adopt a divided-and-conquer approach to regularize the displacement field for regions including the lung surface, left rib, right rib, spine, and the rest regions using smooth Welsch's functions. By controlling the sensitivity to sliding motion in different regions through the Welsch parameters, we align the priori assumptions of the displacement field more closely with reality. Experimental results on the DIR-Lab 4DCT dataset demonstrate the superiority of our method over global methods.


INTRODUCTION
The registration of thoracic CT images serves as the foundation for various clinical applications, such as surgical planning, surgical guidance, and image-guided radiation therapy [1,2].The problem of image registration can generally be formulated as an optimization problem to find the optimal displacement field that maximizes the similarity between the moving image and the fixed image while maintaining a physically plausible displacement field.However, it poses a challenge to ensure the displacement field aligns with reality due to different motion patterns within the chest cavity [3].The common priori assumption for the displacement field is that the displacement vectors of points in local regions should tend to be similar.Therefore, the squared 2 -norm regularization is often employed to penalize the gradient of the displacement field.Although this regularization is computationally simple, it imposes strict constraints on sliding motion at areas like the lung surface due to excessive penalization of gradients with larger magnitudes compared to those with smaller magnitudes, resulting in registration results that do not align with reality at sliding regions.Previous studies have introduced -norm (p ≤ 2) regularization, such as Total Variation (TV, p=1) [4] and Isotropic Total Variation (pTV-Reg, p=2) [5], to allow for more sliding since these norms promote sparsity and reduce the difference in penalization between gradients with larger and smaller magnitudes.However, these methods are global and lead to insufficient penalization of sliding in bone regions that tend to exhibit rigid motion, resulting in registration results that do not align with reality in the bone regions.Additionally, these methods suffer from non-smoothness, making optimization more challenging [6].
In this work, we propose a divided and conquering approach for thoracic CT image registration, utilizing Welsch's functions with varying sparsity to handle the displacement field with different degrees of sliding in different regions of the thorax.For rigid regions, we apply a higher penalty to the gradient of the displacement field, while for sliding regions, we apply a lower penalty.Compared to non-smooth -norm (p ≤ 2) regularization, the Welsch's function used for lower gradient penalization is smooth [6], facilitating optimization.We compare our proposed method with global squared 2norm and 2 -norm regularization on the DIR-Lab 4DCT dataset, and experimental results demonstrate superior performance in terms of Target Registration Error (TRE), lung Dice, and bone Dice.These results validate the effectiveness of our divided and conquering approach.

METHODOLOGY 2.1 Problem Formulation
Let Ω = (x, y, z)|0 ≤ x ≤ X, 0 ≤ y ≤ Y, 0 ≤ z ≤ Z denotes the domain of the CT image.The primary objective of the registration process is to compute a 3D displacement field d : Ω → R 3 that effectively aligns the moving image I M : Ω → R with the fixed image I F : Ω → R. The estimation of displacement field is typically addressed as an optimization problem, formulated as follows: where D represents the term for measuring image dissimilarity, and R denotes the term for regularization of displacements, with serving as the parameter governing the degree of regularization.

Dissimilarity Metric Term
The dissimilarity term in the objective function is used to penalize the dissimilarity between the registered moving image and the fixed image, guiding the iterative process towards reducing the dissimilarity.In this work, Mean-Square Error (MSE) is employed as the dissimilarity metric, defined as follows: where x ∈ R 3 is the voxel position in the image, and d(x) is the displacement vector at x.

Displacement Regularization Term
Regularization of displacement fields penalizes the gradients of the displacement field, forcing neighboring regions to be as similar as possible.This is consistent with the smooth motion observed in most areas of the thorax, such as the lungs and liver.However, there are still some regions that do not conform to this priori assumption, such as the sliding of the lung surface within the thorax during the breathing process [3].Additionally, due to the rigid nature of the spinal region, the displacement vectors of all points in this region should be kept as consistent as possible.In the case of the rib region, there may be partial non-rigid displacement associated with lung deflation, but there is a tendency for rigid motion on each side of the ribs.Therefore, different displacement field regularization should be designed for these regions.
A commonly used method is to apply global squared 2 -norm regularization to the gradients of the displacement field.However, this regularization leads to excessive penalties for sliding motion, hindering the correct alignment of sliding regions [5].To address sliding motion, previous research has proposed using -norm (p ≤ 2) displacement field regularization that is relatively insensitive to sliding motion, such as 1 -norm TV algorithm [4] and 2 -norm pTV-reg algorithm [5].These regularizations reduce the difference in loss between displacement field gradients with large values and those with small values, allowing for more sliding.However, their optimization is challenging due to non-smoothness [6].
In this work, we employ the smooth Welsch's function to penalize the gradients, utilizing different Welsch parameters to control the degree of penalty for sliding.The Welsch's function is expressed as: where > 0 is a parameter to control the degree of sliding penalty.Then, the regularization term for the displacement field can be expressed as follows: where N = {0, 1, 2, 3} represents the index set for different regions, namely the lung surface, the spinal region, the left rib region, the right rib region, and the rest region.corresponds to the weight assigned to the displacement field regularization in region .Higher weights can prioritize the regularization of the displacement field in specific regions.represents the Welsch parameter that controls the penalty for sliding in region .d denotes the displacement field in region i, and ∇d (x) represents the gradient vector of the displacement field at point x.
The function graph of the Welsch's function with different parameters exhibits distinct characteristics, as shown in Figure 1.
It can be observed that the function monotonically increases in the interval [0, +∞] and has an upper bound of 1.Additionally, as approaches 0, Equation 2 tends to the 0 -norm [6].Therefore, by controlling the value of , we can regulate the degree of penalty for sliding motion.In this work, we set different values for different  regions.For the spine, where rigid displacement occurs, we set to a larger value of 4. For the left and right ribs, which undergo slight smooth motion, is set to 2. For the lung surface, where sliding motion occurs, is set to a smaller value of 0.1.For the rest region with smooth motion, such as the lungs and liver interior, is set to 0.5.Although a low value can reduce sensitivity to larger loss values, it may lead to an imbalance in the weights of the terms in the objective function which consists of Welsch's functions with various values, penalizing terms with low loss values more than those with higher values.One solution is to perform registration with five different values separately, followed by segmentation and concatenation of the displacement fields.However, this approach can introduce overlapping problems at the boundaries of the segments [7].In this study, we address this problem by setting the overall weight w i of each region to i , aiming to balance the weights across different regions.The flowchart of the registration algorithm is illustrated in Figure 2.

Regions Segmentation
Figure 3 illustrates the segmentation process for each region.Firstly, the lung mask is obtained using a U-Net based deep learning model proposed by Hofmanninger et al [8].Then, the lung surface is extracted through edge detection and dilation operations.Subsequently, the bone region mask is easily obtained by thresholding, as bones typically have higher intensity values.Next, the dilated lung mask was subjected to a AND operation with the bone mask to obtain the left rib, spine, and right rib regions.Finally, the masks for the remaining regions are obtained by subtracting the mask of the rest region from the entire image.

Data Preparation
Our method was evaluated using the DIR-Lab 4DCT dataset [9] to assess its accuracy.This dataset consists of 10 cases, each containing 3DCT images from ten different phases.We performed registration between the end-expiration (EE) and end-inspiration (EI) time phases, which exhibit the maximum differences.All CT images were resampled to an isotropic voxel size of 1 × 1 × 1 mm 3 and redundant regions were removed.The intensity values were cropped to the range [0, 2000] and then scaled to the range [0, 1].

Implementation Details
The AMSGrad optimizer [10] was employed to perform gradient descent on the objective function.The initial learning rate of the optimizer was set to 0.005.The convergence threshold was set to 10 -5 .The maximum number of iterations was set to 800.

Comparision Methods
We compared our method with commonly used global squared 2norm regularization and 2 -norm regularization used by pTV-Reg.The evaluation metrics for comparison were the Target Registration Error (TRE) of the 300 landmarks provided by the DIR-Lab 4DCT dataset, as well as the Dice coefficients for the lung and bone regions.
The Dice coefficient is calculated using the following formula:

Results
Table 1, Table 2 and Table 3 presents a comparison between our method and the comparative methods in terms of Target Registration Error (TRE) and Dice coefficients for the lung and bone regions.It is evident that our method outperforms the comparative methods in all the evaluated metrics.

DISCUSSION
The Welsch's function has been widely used in robust image processing tasks due to its ability to adjust the sensitivity to noise and outliers through its parameter while maintaining function smoothness [6,11].In fact, the Welsch's function penalizes different losses in a regularized term with Gaussian weights, assigning weights to all loss terms with a mean of 0 and a standard deviation of [11].According to the 3-rule, terms with losses exceeding 3 are assigned small weight and thus excluded from optimization, making the optimization process insensitive to high-loss terms.
Since different regions of the thorax exhibit different motion patterns during respiration, the assumption of completely smooth motion (such as squared 2 -norm regularization) is not appropriate.Squared 2 -norm regularization imposes strict constraints on sliding  [5] is less sensitive to sliding motion, as it reduces the penalty difference between large and small loss terms through the square root operation.However, this may lead to insufficient penalization of rigid motion in the bone region.In this work, we employ the Welsch's functions with different parameters to penalize displacement field gradients in different regions, effectively handling regions with varying degrees of sliding motion and aligning the priori assumptions with reality.
The results demonstrate that our method outperforms 2 -norm regularization and squared 2 -norm regularization in terms of target registration error, lung Dice coefficient, and bone Dice coefficient, thus validating the effectiveness of our approach.Furthermore, as expected, the lung Dice coefficient with 2 -norm regularization is higher than that with squared 2 -norm regularization, indicating that reducing sensitivity to sliding motion is beneficial for lung registration since the lung surface undergoes sliding motion within the chest cavity.Conversely, the bone Dice coefficient with squared 2 -norm regularization is higher than that with 2 -norm regularization, suggesting that reducing sensitivity to sliding motion is detrimental to bone registration as the bone is rigid and experiences minimal sliding motion.

CONCLUSION
In this study, we propose a divide-and-conquer approach for displacement field regularization in thoracic CT image registration.We introduce a regularization scheme that utilizes smooth Welsch's functions with different parameters to control the penalty for sliding motion in different regions, aligning the priori assumptions with reality.The effectiveness of this approach is validated through experimental evaluation.

Figure 1 :
Figure 1: function with different values of the parameter .

Figure 2 :
Figure 2: The workflow diagram of the proposed algorithm.

Table 1 :
Comparison of Mean Target Registration Error (TRE) between Our Method and Comparative Methods

Table 2 :
Comparison of Lung Dice between Our Method and Comparative Methods

Table 3 :
Comparison of Bone Dice between Our Method and Comparative Methods