Stepwise local stitching ultrasound image algorithms based on adaptive iterative threshold Harris corner features.

Herein, a Harris corner detection algorithm is proposed based on the concepts of iterated threshold segmentation and adaptive iterative threshold (AIT-Harris), and a stepwise local stitching algorithm is used to obtain wide-field ultrasound (US) images.Cone-beam computer tomography (CBCT) and US images from 9 cervical cancer patients and 1 prostate cancer patient were examined. In the experiment, corner features were extracted based on the AIT-Harris, Harris, and Morave algorithms. Accordingly, wide-field ultrasonic images were obtained based on the extracted features after local stitching, and the corner matching rates of all tested algorithms were compared. The accuracies of the drawn contours of organs at risk (OARs) were compared based on the stitched ultrasonic images and CBCT.The corner matching rate of the Morave algorithm was compared with those obtained by the Harris and AIT-Harris algorithms, and paired sample t tests were conducted (t = 6.142, t = 31.859, P < .05). The results showed that the differences were statistically significant. The average Dice similarity coefficient between the automatically delineated bladder region based on wide-field US images and the manually delineated bladder region based on ground truth CBCT images was 0.924, and the average Jaccard coefficient was 0.894.The proposed algorithm improved the accuracy of corner detection, and the stitched wide-field US image could modify the delineation range of OARs in the pelvic cavity.

Introduction

Cone beam computer tomography (CBCT) is a common image guiding method and has been extensively used in adaptive radiotherapy (ART) of pelvic tumors.[ However, the soft tissue resolution of CBCT is low,[ and the artifacts caused by the hardening effect and electron scattering worsen the quality of the reconstructed image.[ Ultrasound (US) possesses numerous advantages, including portability, use of nonionizing radiation, and real-time performance.[ It is commonly used in the monitoring of ventricular mechanical performance, cardiac motion estimation, and cardiovascular boundary detection, among others.[ In the radiotherapy of pelvic tumors, US images can be combined with CBCT images to assist in image-guided treatment.[ Arcadipane et al[ demonstrated that 3D US images could clearly show the contours of the anterior rectal wall and bladder neck and reduce the error of delineation of the prostate, bladder, and other organs in prostate cancer radiotherapy. However, for patients with pelvic tumors, in the case of high bladder filling, it is difficult to include the gross tumor region (GTV) and complete organs in the imaging range with narrow beam US images, which is not conducive to clinicians’ judgment of tumor location or organs at risk (OARs). At the same time, incomplete organ structure imaging also has a certain degree of impact on the set up accuracy of radiotherapy. To solve the above problems, narrow-beam US images are stitched into wide-field images through image stitching technology, so that the GTV and the complete tissues or organs are completely included in the imaging range. Wide-field US images can better present the contour edge of organs and relative position information between different organs, so that doctors can make more accurate delineation of target areas and OARs, and also help to improve the accuracy of the set up. Therefore, US combined with CBCT image guidance method can be better applied in the treatment of pelvic tumors.

Image stitching refers to the combination of a group of overlapping images based on image registration and image fusion and the generation of a new wide-field image, which includes all the pixel information of the images before stitching.[ One of the key steps in image stitching is the extraction of corner features. Corner features refer to the points of the maximum curvature on the edge curve or the pixels of 2D image with drastic gray change. These features can effectively reduce the amount of image information, effectively improve the speed of calculation, and facilitate the accuracy of images matching. Accordingly, corner feature detection methods can be classified into the following 2 categories: corner detection based on image edges[ and corner detection based on image grayscale values.[ The first type mainly determines corner points according to the detected edge intersection points, which mainly relies on image segmentation and edge extraction. It involves a large amount of calculations and a small fault tolerance rate, and its application scope is limited. Edge detection operators mainly include the Laplacian of Gaussian, Canny, and Sobel operators.[ The second type searches for pixels that correspond to the maximum values of local grayscale changes and sets them as corner points. Gray-level detection operators mainly include the Moravec, Harris, and Susan operators.[ The Harris operator is a common method based on image grayscale detection corners. Harris corner detection algorithm is a signal-based point feature extraction method that is based on the Moravec algorithm. It uses a fixed window to move in any direction in the image and compares the degree of grayscale change in the window before and after sliding to determine the corner points.[ A schematic of the Harris corner detection method is shown in Figure 1. If the window moves in any direction without grayscale change, it is the flat region. However, if the window moves in a certain direction, the grayscale of the image in the window changes greatly, and the grayscale in other directions does not change, it is the edge region. If the window moves in any direction and the gray scale changes obviously, it is the corner region. The Harris corner features are insensitive to image rotation and grayscale and noise changes and is relatively simple to calculate. It can extract a large number of feature points from complex human tissue images.[ However, the traditional Harris corner detection algorithm needs to manually set the threshold when corner points are extracted, which is not scale invariant. If the threshold setting is too large, it will lead to the loss of some corners, and if it is too small, it will lead to a large number of pseudo corners.[ Previous studies have improved the threshold selection in the case of the traditional Harris corner detection method, but they are not ideal in terms of robustness and self-adaptation, and the calculation amount is also large.[

Figure 1

Harris corner detection schematic. (A) Flat region; (B) edge region; (C) corner region.

In this study, the segmentation method based on iterative threshold was proposed to extract the corner features of abdominal US images. The selected features were more accurate and adaptive for identifying the threshold of corner detection. According to the precision of the extracted corner features, the methods of stepwise local brute-force (BF) and random sample consensus (RANSAC) matching were used to stitch the abdominal ultrasonic images obtained from different scanning angles. In this respect, the stitched ultrasonic images provided a wider field of vision than the scope of the probe. Finally, the delineated bladder region was compared with that in the CBCT image to prove its applicability in the delineation of OARs in abdominal tumor radiotherapy. The flow chart of wide-field US acquisition is shown in Figure 2.

Figure 2

Flow chart of wide-field US acquisition.

Methods

Data

The data included 9 patients with cervical cancer who did not undergo surgery and 1 patient with prostate cancer. The protocol of this study was approved by the IRB of Second People's Hospital of Changzhou, Nanjing Medical University, Changzhou, China (Ethical no.:2019KY056-01). All experiments were performed in accordance with relevant guidelines, and the informed consent was obtained from all participants.

CBCT and US image acquisition

Before image acquisition, the indoor laser lamp coordinate system and US probe were calibrated with a calibration phantom. This operation assured the same coordinate positions for each 3D ultrasonic voxel point. After calibration, the CBCT of Infinity Linear Accelerator (Elekta Company, Sweden) was used to scan the patients. The scanning conditions included 120 kV tube voltage, 20 mA tube current, and 1 × 1 × 3 mm3 voxels. US images were acquired using the Clarity US-guided device (Elekta Company, Sweden). Scanning was conducted with a C5–2/60 US probe and a central frequency of 3.5 MHz. The voxel size was 1 × 1 × 3 mm3. US images were scanned immediately after acquisition of CBCT images to reduce errors caused by the movement of OARs. The acquisition interval between CBCT and US images was <2 minutes. During ultrasonic scanning, the ultrasonic probe was attached to a robotic arm that performs constant pressure scans over a predetermined orbit and kept in a perpendicular orientation to the body surface of patients. The patients were scanned in 2 directions (75 and 105 degree relative to the coach) in succession to avoid bone tissues.

Harris corner detection algorithm based on an iterative threshold method

Two groups of 3D US images for stitching were Gaussian filtered and resampled. They were then registered with CBCT images acquired during the same treatment stage with automatic affine registration. The algorithm found the most likely correspondences between the bladder contour of both images in a noniterative manner based on a modified Viterbi algorithm.[ Two sets of 2D US slices after registration were aligned in sequence. This avoided the problem that arises when directly acquired 2D image pairs are not at the same scanning level and have specific angular differences. Based on iterative threshold segmentation, we proposed an adaptive corner response threshold method to select corner features[ and improve the traditional Harris corner feature detection algorithm. The are 4 steps. First, a 5 × 5 test window was set in the 2D US image and centered on specific pixels within the image. The test window was then moved slightly in all spatial directions. If the window was too large, the convolution operation process could lead to corner point offsets during the calculation of the image pixel adaptive matrix. Conversely, if it was too small, a large number of pseudo corner points would be extracted from the US image owing to the influence of noise points. Second, the angular point response value within the window was calculated as follows:

where R is the response value of the corner point; r and c represent the horizontal and vertical coordinates of the image pixel points, respectively; and μ1 and μ2 are the singular values of the autocorrelation matrix. Each pixel in the US image had a specific R value. Third, the initial corner response threshold T was set, and the pixels in the image were divided into groups G1 and G2. G1 was composed of pixels with R values greater than T, and G2 was composed of pixels with R values less than T. Fourth, the average R values m1 and m2 of G1 and G2 classes were calculated, and the new threshold Tk = (m1 + m2)/2 was obtained and compared with the initial T value. Lastly, the 4th step was repeated until the difference ΔT between the Tk values in the successive iterations became smaller than the error value, and the optimal threshold T was obtained. The process of obtaining the optimal threshold based on the iterative threshold segmentation method is depicted in Figure 3.

Figure 3

Flow chart of optimal corner response threshold based on iterative threshold segmentation method.

Stepwise local US image stitching algorithms

After obtaining the corner features, the US images were stitched based on the 8 steps. First, to improve the efficiency of corner feature matching and avoid obvious corner mismatching, a 60 × 150 local corner matching window was set. Second, a multiscale-oriented patch descriptor was constructed with each corner as the center. The descriptor was invariant and could be used to record the regional feature information around the corner. Third, the BF method was used to match the 2 US images for stitching. The Euclidean distance between the descriptors was calculated, and the matched corner pairs were obtained. Fourth, the RANSAC algorithm was used to purify the angular pairs obtained by BF, and point pairs with false matches were eliminated. The calculated mean distance between the angular matching corner pairs could be used as the optimal stitching range for the local US images.[ Fifth, the steps listed above were repeated, and the corresponding weight value was set based on the ratio of the corner pairs in the windows to the total corner matching pairs of the entire image. Then, the optimal matching range of the entire image was calculated according to the following equation:

where i represents the i-th corner matching window, the range is (0, 1), ω is the weight, n is the number of matching points among the corner points, d is the optimal matching range of the local US images, and d represents the optimal matching range of the entire image. Sixth, the progressively evolving fusion algorithm was used to stitch the US images to form a 2D wide-field US image. Seventh, the stitched 2D images were reconstructed into a 3D wide-field US image. Eighth, the reconstructed 3D US and CBCT images were registered by global mutual information and local deformation registration methods. The bladder region in the US image was extracted using regional saliency and used for local deformation registration. The whole registration process would alternate between global and local registrations until it converged.[ Local corner matching results of US images are shown in Figure 4A–C show the pairs of US images to be matched after corner detection. The red points are corners, and the red rectangular window is a local corner matching window, which limits the matching range of corner features. Figure 4D–F are the corresponding results of the local matching of corner points of the US images. Every 2 matched corner points in the 2 US images are connected by a color line. It can be seen that in each local region, most of the corner points can be matched correctly, which can improve the accuracy of image stitching.

Figure 4

Local corner matching of US images. (A–C) The US images used for stitching. (D–F) Local corner matching results of the US images.

Evaluation

The accuracy of the US image stitching process can be measured using the matching rate, which can not only measure the reliability of the corner feature detection algorithm but also demonstrate the accuracy of the matching algorithm. In this experiment, the corner matching rate of the ultrasonic image was defined as follows:

where N is the number of correctly matched corner pairs in all slices, which includes the bladder region of the tested patients, and represents the sum of the corner matching pairs that were not eliminated after purification by the RANSAC algorithm. N is the number of mismatched corner point pairs, and N + N represents the sum of the unpurified pair of supporting parts after BF stitching.

In this study, the Dice similarity coefficient (DSC)[ and Jaccard similarity coefficient[ were used to evaluate the segmentation accuracy of the bladder region based on the stitched wide-field US images. The ground truth was defined by the contour of the bladder, which was manually delineated by medical physicists in CBCT images. By calculating the overlap of the bladder volume, an accurate segmentation result outcome should yield a high organ volume overlap rate. The DSC and Jaccard coefficients were defined as follows:

where L and L represent the results of binary segmentation of the bladder in CBCT and US images, respectively. As the DSC and Jaccard coefficient values approach unity, the automatically segmented bladder region (from the stitched US image) becomes progressively more similar to that outlined in the CBCT image.

Results

All images used in this study were obtained from 9 cervical cancer patients and 1 prostate cancer patient. The wide-field US images obtained after local stitching based on corner features extracted by the AIT–Harris, Harris, and Morave algorithms are shown in Figure 5. The corner matching rates obtained by the three algorithms are listed in Table 1. The average corner matching rate based on the AIT–Harris algorithm was the highest. In addition, the corner matching rate obtained by the Morave algorithm was defined as the ground truth and was compared with that obtained by the Harris and AIT–Harris algorithms with estimated paired sample t test values of t = 6.142 and t = 31.859, respectively (P < .05). Compared with the Harris algorithm, the corner matching rate obtained by the improved algorithm was significantly different from that obtained by the ground truth.

Figure 5

Results of local corner point matching in US images. Symbols a1 and a2 refer to the US images for stitching, whereas b1,2, c1,2, and d1,2 are the wide-field ultrasonic images spliced based on corner features extracted by the AIT–Harris, Harris, and Morave algorithms, respectively.

Table 1

Matching results after corner feature extraction based on 3 corner detection algorithms.

In this study, CBCT was used as a registration reference image. Compared with the affine registration results of CBCT and unstitched US images, the same registration method was used to quantitatively evaluate the improvement of the registration results of wide-field US and CBCT. The registration errors between US and CBCT are shown in Table 2. US (75) and US (105) refer to the US images obtained from 2 different angles (75 and 105 degree, respectively). US(matching) refers to the stitched US images. The mean registration errors of the stitched US images and CBCT are 0.33, 0.38, and 0.32 in the lateral, longitudinal, and vertical directions, respectively. Compared with the registration results of unstitched US images and CBCT, the improvement effect is significant.

Table 2

The registration errors between US (US(75), US(105), and US(matching)) and CBCT images.

Figure 6 shows a comparison of bladder contours. Figure 6A shows a comparison of bladder contouring results obtained from the patient with prostate cancer. Figs. 6B and C show another comparison of bladder contouring results obtained from cervical cancer patients. The yellow line is the outline of the bladder based on the CBCT image and serves as the ground truth. The white line is the bladder contour extracted by using the concept of region saliency based on the stitched US images. Two quantitative indices, namely, the DSC and Jaccard coefficients, were used to evaluate the similarity between the automatically sketched bladder contour based on US image and CBCT ground truth contour. The experimental results are listed in Table 3. The mean DSC measured between the segmented bladder regions on the US image obtained based on the improved algorithm splicing and the bladder region manually outlined on the CBCT image was 0.924, and the mean value of the Jaccard coefficient was 0.894.

Figure 6

Comparison of bladder contours. The yellow line is the outline of the bladder based on the CBCT image, and the white line is the bladder contour extracted after segmentation based on the stitched US images. (A) Comparison of bladder contouring results obtained from the prostate cancer patient. (B and C) Comparison of bladder contouring results obtained from cervical cancer patients.

Table 3

Evaluation results of the DSC and the Jaccard coefficient based on the automatic segmentation of bladder contours on ultrasound images and the manual delineation of ground truth contours on CBCT.

Discussion

In this study, the Harris corner feature detection algorithm was improved using the iterative threshold segmentation method. Additionally, the stepwise local RANSAC corner matching method was used to stitch the US images without using the wide-field imaging function. US images acquired at one angle did not show larger GTV and OARs within the imaging range. The stitched wide-field US images also facilitated the delineation of OARs in SBRT and ART in CBCT images. The experimental results showed that the corner feature obtained by the AIT–Harris corner detection method is more accurate. As can be observed from the comparison of the contours of the bladder area, the US image has high contrast in the soft tissue area and clearly distinguishes the contours of the bladder and surrounding tissues. According to the results of the DSC and Jaccard coefficients, the bladder region in the wide-field US image is similar to that delineated by CBCT; however, their contours at the junction of the bladder and surrounding tissues and organs are significantly different. The contour of the bladder region based on wide-field US could be corrected in the CBCT images.

The corner detection threshold in the traditional Harris corner detection algorithm is an empirically derived value and can only be applied to the detection of corner features in some specific images. Accordingly, it cannot be used universally.[ The structures of tissues and organs in the lower abdomen are quite complex. As such, many corner points and edge features can be observed in the US image of this region. A single, fixed threshold T would yield a large number of invalid points while detecting corner points. To improve the accuracy, the threshold of manually repeated debugging is needed, which reduces the efficiency of corner detection. In this study, the adaptive iterative threshold Harris corner detection algorithm was improved. Based on iterative threshold segmentation, the corner response R values of all pixels in the image were substituted for the pixel value. Additionally, repeated iterative calculations were performed to obtain the optimal corner response threshold by replacing the threshold value of the traditional algorithm—which was determined based on experience—so as to provide accurate corner features for the subsequent stitching process. This algorithm avoided the corner response threshold, which was set based on experience, and extracted corner features more accurately. Meanwhile, it reduced the extraction of pseudo corners in US images and improved the adaptability of traditional Harris corner detection methods. Each slice has its own optimal threshold to extract corner features, and the adaptive effect is obvious. In global image matching, the mapping function needs to traverse all pixel points in the image. Regarding redundant corner features in the overlapped areas outside the 2 US images for stitching, the mapping function interfered with the normal matching of the corner features inside the overlapped areas of the US images, resulting in mismatching.[ Based on the overlapping area of 2 US images for stitching, a local matching window is set, and the entire overlapping area is traversed in a stepwise manner to limit the stitching range. This operation can reduce the erroneous matching of corner pairs in the matching algorithm. In this experiment, the accuracy of corner feature extraction and image stitching algorithms based on US images was preliminarily studied in 9 patients with postoperative cervical cancer and 1 prostate cancer. The results proved the application prospect of US images in the auxiliary delineation of OARs in the RT of clinical abdominal tumors. In future research, we will continue to generate new data and verify the application prospect of the algorithm. Considering the abundance of corner feature information in 3D ultrasonic images and the large values of the calculated corner response R, the operation time of corner detection based on 3D images is very long, ultimately affecting the efficiency of image stitching. Therefore, this study adopted the 2D ultrasonic image stitching algorithm. In future experiments, we will conduct an in-depth study on the algorithm based on a deep learning method to extract the corner features of US images so as to improve the accuracy and real-time performance of the stitching process.

Conclusions

In this study, a corner detection method was proposed based on iterative threshold segmentation, and the accuracy of the new method was verified experimentally. Based on the preliminary study on the accuracy of the stepwise local US stitching algorithm, narrow beam US images are stitched into wide-field images, so that the GTV and organs are completely included in the imaging range, which is convenient for doctors to make more accurate delineation of target areas and OARs, and also helps to improve the accuracy of set up. This makes US combined with CBCT image guidance method better applied in the treatment of pelvic tumors.

Author contributions

Hongfei Sun is a graduate student at Northwestern Polytechnical University. He currently pursuing the Ph.D. degree with the School of Automation, Northwestern Polytechnical University, China. In this study, he performed the experimental characterizations and analyzed the results.

Jianhua Yang received the Ph.D. degree from Northwestern Polytechnical University. She is currently a Professor with the School of Automation, Northwestern Polytechnical University. Her research interests include biomedical image processing, detection, and control technology. In this study, she revised the manuscript.

Rongbo Fan is a graduate student at Northwestern Polytechnical University. He currently pursuing the master's degree with the School of Automation, Northwestern Polytechnical University, China. In this study, he preprocessed the images and assisted to perform the experimental characterizations.

Kai Xie is an assistant research fellow at the Affiliated Changzhou No.2 People's Hospital of Nanjing Medical University. In this study, he focused on the delineation of the organ at risk.

Conghui Wang is a graduate student at Northwestern Polytechnical University. She currently pursuing the master's degree with the School of Automation, Northwestern Polytechnical University, China. For this study, she conceived the experiments.

Xinye Ni received the Ph.D. degree from Nanjing University of Aeronautics and Astronautics, Nanjing, with an expertise in radiotherapy and biomaterials. In this study, he performed the ultrasound images segmentation and revised the manuscript.