CN104299209A - Lung 4D-CT image super-resolution reconstruction method based on fast sub-pixel motion estimation - Google Patents

Abstract

The invention discloses a lung 4D-CT image super-resolution reconstruction method based on fast sub-pixel motion estimation, including step (1) reading the initial lung 4D-CT image, the lung 4D-CT image consists of multiple different phases 3D-CT images of the lungs, and randomly select a phase of the 3D-CT images of the lungs as the 3D-CT images of the lungs to be reconstructed; (2) remove the 3D-CT images of the lungs to be reconstructed from the 4D-CT images of the lungs The lung 3D-CT image after the CT image is compared with the lung 3D-CT image to be reconstructed for initial motion estimation, and the initial motion vector field between them is obtained; (3) The precision optimization of the obtained initial motion vector field is obtained to obtain Sub-pixel motion vector field; (4) Based on the sub-pixel motion vector field, the 3D-CT image of the lung to be reconstructed is reconstructed to obtain a reconstructed high-resolution image with the same phase as the 3D-CT image of the lung to be reconstructed Rate lung 4D-CT images. This method can improve the resolution of lung 4D-CT images.

Description

Super-resolution reconstruction method of lung 4D-CT images based on fast sub-pixel motion estimation

technical field

The invention relates to the technical field of medical image processing, in particular to a method for super-resolution reconstruction of lung 4D-CT images based on fast sub-pixel motion estimation.

Background technique

Lung 4D-CT images play an important role in lung cancer radiotherapy, it can provide a comprehensive high-precision radiotherapy respiratory motion characterization, help to track tumor movement, implement precise radiotherapy, and reduce damage to normal tissues . However, due to the inherent high-dose radiation of CT, it is often only possible to reduce the sampling along the longitudinal direction (Z-axis direction) to reduce the scanning time of lung 4D-CT in order to reduce the radiation dose, resulting in a poor inter-layer resolution of lung 4D-CT images. Lower than the intra-layer resolution, resulting in significant anisotropy of the data. This makes it necessary to perform an interpolation operation to obtain a correct display when multi-plane observation (coronal sagittal plane, etc.) is performed on the data, and this operation easily leads to blurring of the image.

During super-resolution reconstruction, the estimation of the motion field between images is the main factor affecting the reconstruction accuracy and speed. The super-resolution technology based on motion estimation previously proposed by the applicant uses a full-search motion estimation method to estimate the motion field between different frame images. This method can reconstruct a clearer coronal (sagittal) image of the lung than the traditional interpolation method, but its main disadvantages are: slow speed and only a fixed search step. This not only affects the reconstruction speed, but also affects the reconstruction accuracy. Therefore, it is necessary to provide a super-resolution reconstruction method of lung 4D-CT images based on fast sub-pixel motion estimation to overcome the shortcomings of the existing technology.

Contents of the invention

The purpose of the present invention is to provide a lung 4D-CT image super-resolution reconstruction method based on fast sub-pixel motion estimation, which can improve the resolution of the lung 4D-CT image.

The object of the present invention can be achieved through the following technical measures: a lung 4D-CT image super-resolution reconstruction method based on fast sub-pixel motion estimation, the method comprising the following steps:

(1) Read the initial 4D-CT image of the lung, which is composed of multiple 3D-CT images of the lung with different phases, and arbitrarily select a 3D-CT image of the lung with a certain phase as the Reconstruction of 3D-CT images of the lungs;

(2) Perform initial motion estimation on multiple lung 3D-CT images with different phases after removing the lung 3D-CT image to be reconstructed in the lung 4D-CT image relative to the lung 3D-CT image to be reconstructed, and obtain the lung The initial motion vector field between multiple different phase lung 3D-CT images after removing the lung 3D-CT image to be reconstructed and the lung 3D-CT image to be reconstructed in the first 4D-CT image, the initial motion vector field The precision is an integer;

(3) Carry out precision optimization to the obtained initial motion vector field, make it accurate to sub-pixel, obtain sub-pixel motion vector field;

(4) Based on the sub-pixel motion vector field obtained in step (3), the 3D-CT image of the lung to be reconstructed is reconstructed to obtain a reconstructed high-resolution image with the same phase as the 3D-CT image of the lung to be reconstructed. Rate lung 4D-CT images.

In the present invention, in the step (2), a three-step search method is used for initial motion estimation, and multiple lung 3D-CT images of different phases after removing the lung 3D-CT image to be reconstructed are searched for in the lung 4D-CT image. The integer pixel displacement between the image and the 3D-CT image of the lung to be reconstructed. The algorithm uses a coarse-to-fine search pattern, starting from the center point of the search window, taking 8 surrounding points at a certain step to form a point group for each search, and then performing matching calculations according to the matching criteria to find the center of the matching block with the smallest error point. The specific process includes:

(2.1) Determine a center point, set the maximum search length, and use 1/2 of the maximum search length as the step size, and find the minimum block error point according to the matching criteria for the center point and 8 detection points with the same step length around the distance, If the minimum block error point is located at the original center point, the algorithm ends, otherwise, proceed to step (2.2);

(2.2) The step size is halved, and the center point of the minimum error matching block is found in the minimum block error point determined in the previous step and 8 detection points of the same step length around it;

(2.3) Repeat (2.1) and (2.2) until the step size meets the search accuracy requirement, that is, the best matching point is obtained.

What the present invention adopts is that the minimum absolute error matching criterion commonly used is defined as:

$S(v_x,v_{\mathrm{the\; y}})={\mathrm{\Σ}}_{i=0}^{N-1}{\mathrm{\Σ}}_{j=0}^{N-1}|I_t(p+i,q+j)-I_{t-{\mathrm{\Δ}}_t}(p+v_x+i,q+v_{\mathrm{the\; y}}+j)|$ ......Formula 1)

In the formula, the size of the image block is N×N, the coordinates of the upper left corner are (p, q), and the block size selected by the present invention is 16×16; the motion vector is (v _x , v _y ); I _t (i, j) and are the values at pixel (i, j) of the current frame and the reference frame, respectively. The motion vector with the smallest value of S(v _x , v _y ) in the search area is the optimal motion vector of the current block.

In the present invention, in the step (3), the optical flow method is used to optimize the accuracy of the initial motion vector, which is accurate to sub-pixels. The specific process includes:

(3.1) Assuming that the two consecutive images before and after motion are f(x,y) and g(x,y) respectively, according to the simplified model of classical optical flow, we can get

$g(x,\mathrm{the\; y})=f(x+\mathrm{\Δ}{x}_{\mathrm{the\; s}},\mathrm{the\; y}+\mathrm{\Δ}{\mathrm{the\; y}}_{\mathrm{the\; s}})\≈f(x,\mathrm{the\; y})+\mathrm{\Δ}{x}_{\mathrm{the\; s}}\frac{\∂}{\∂x}f(x,\mathrm{the\; y})+\mathrm{\Δ}{\mathrm{the\; y}}_{\mathrm{the\; s}}\mathrm{\Δ}{x}_{\mathrm{the\; s}}\frac{\∂}{\∂\mathrm{the\; y}}f(x,\mathrm{the\; y})$ ...Formula (2)

This formula is approximately a first-order Taylor series expansion.

(3.2) Minimize and solve formula (2) to obtain the optimal displacement vector:

$\underset{\mathrm{\Δ}{x}_{\mathrm{the\; s}},\mathrm{\Δ}{\mathrm{the\; y}}_{\mathrm{the\; s}}}{\min \mathrm{imize}}\mathrm{\Φ}(\mathrm{\Δ}{x}_{\mathrm{the\; s}},\mathrm{\Δ}{\mathrm{the\; y}}_{\mathrm{the\; s}})$ ... Formula (3);

in:

$\mathrm{\Φ}(\mathrm{\Δ}{x}_{\mathrm{the\; s}},\mathrm{\Δ}{\mathrm{the\; y}}_{\mathrm{the\; s}})=\underset{x,\mathrm{the\; y}}{\mathrm{\Σ}}{(g(x,\mathrm{the\; y})-f(x,\mathrm{the\; y})-\mathrm{\Δ}{x}_{\mathrm{the\; s}}\frac{\∂}{\∂x}f(x,\mathrm{the\; y})-\mathrm{\Δ}{\mathrm{the\; y}}_{\mathrm{the\; s}}\frac{\∂}{\∂\mathrm{the\; y}}f(x,\mathrm{the\; y}))}^2$ ...Formula (4);

(3.3) Equation (4) can be regarded as a linear least squares solution problem, and the partial derivatives of the objective function can be set to 0 to obtain the optimal values Δx _s , Δy _s . Therefore, the following equation can be obtained:

$\frac{\∂\mathrm{\Φ}}{\∂\mathrm{\Δ}{x}_{\mathrm{the\; s}}}=0$ and $\frac{\∂\mathrm{\Φ}}{\∂\mathrm{\Δ}{\mathrm{the\; y}}_{\mathrm{the\; s}}}=0$ ...Formula (5);

(3.4) Solving equation (5), the optimal solution Δx _s , Δy _s can be obtained.

(3.5) Suppose the motion vector obtained by the three-step search method is Δx _I , Δy _I , then the final total motion vector Δx, Δy is:

Δx = Δx _I + Δx _s , Δy = Δy _I + Δy _s ...... formula (6);

In the present invention, the step (4) uses the iterative back projection method to reconstruct a high-resolution lung 4D-CT image, the abbreviation of the iterative back projection method is IBP, and the specific process includes:

(4.1) interpolating and enlarging the low-resolution lung 3D-CT image to be reconstructed as the initial high-resolution image H ⁽⁰⁾ ;

(4.2) Simulate the imaging process according to the degradation model to obtain a collection of low-resolution images corresponding to the original low-resolution image sequence K represents the number of images in the original low-resolution image sequence, and n is the number of iterations.

Specifically, in the nth iteration process, H ⁽ⁿ⁾ simulation degenerates to The process is as follows:

$L_k^{\left(\mathrm{no}\right)}=\left(T_k\left(h^{\left(\mathrm{no}\right)}\right)h\right)\↓\mathrm{the\; s}$ ...Formula (7);

Among them, T _k represents the two-dimensional geometric transformation from H to L _k , which is the motion deformation field obtained in step (3); h is the Gaussian blur operator; ↓ _s is the downsampling operator;

(4.3) Judgment error Whether it reaches the minimum value, if reached, then stop the iteration, the previously estimated H ⁽ⁿ⁾ is the final super-resolution image; if not reached, then enter step (4.4);

Judgment error in the above step (4.3) Whether the minimum value is reached can be determined by judging whether the error function e ⁽ⁿ⁾ is smaller than the set threshold ε. The specific calculation formula of the error function is:

$e^{\left(\mathrm{no}\right)}=\sqrt{\frac{1}{K}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\left|\right|L_k-L_k^{\left(\mathrm{no}\right)}\left|\right|}_2^2}$ ...Formula (8);

(4.4) Update the current high-resolution image according to the error, and the update process is as follows:

$h^{(\mathrm{no}+1)}=h^{\left(\mathrm{no}\right)}+\frac{1}{K}{T}_k^{-1}\left(((L_k-L_k^{\left(\mathrm{no}\right)})\&Up\; Arrow;\mathrm{the\; s})p\right)$ ... Formula (9);

In the formula, ↑s represents the up-sampling operator; p represents the back-projection operator, which depends on h and T _k ;

(4.5) Use the updated high-resolution image as the initial high-resolution image, and enter step (4.2).

Compared with the prior art, the present invention has the following beneficial effects:

(1) Due to the optimization of motion vector estimation speed and precision, the calculation time of the method of the present invention is reduced by nearly 20 times compared with the full search algorithm, which is a significant improvement.

(2) In the high-resolution lung 4D-CT image reconstructed by the present invention, the brightness and clarity of blood vessels in the lung parenchyma and surrounding tissues are significantly enhanced.

Description of drawings

The present invention will be described in detail below in conjunction with the accompanying drawings and specific embodiments.

Fig. 1 is the flow chart of the lung 4D-CT image super-resolution reconstruction method based on fast sub-pixel motion estimation of the present invention;

Fig. 2 is the reconstruction result diagram of data 2 phase 0 coronal plane using different methods in Embodiment 1 of the present invention. From left to right, the bilinear interpolation method is used, and the reconstruction result diagram based on the full search motion estimation algorithm and the method of the present invention is used;

Fig. 3 is an enlarged schematic diagram corresponding to the box part in Fig. 2 according to Embodiment 1 of the present invention. From left to right, it is a partial enlarged view of the result map reconstructed based on the full search motion estimation algorithm and the method of the present invention using the bilinear interpolation method. ;

Fig. 4 is a reconstruction result diagram of the data 2 phase 0 sagittal plane using different methods according to Embodiment 1 of the present invention. From left to right, it is the reconstruction result diagram based on the full search motion estimation algorithm and the method of the present invention using the bilinear interpolation method;

Fig. 5 is an enlarged schematic diagram corresponding to the box part in Fig. 4 of Embodiment 1 of the present invention. From left to right, it is a partial enlarged view of the result map reconstructed based on the full search motion estimation algorithm and the method of the present invention using bilinear interpolation method .

Detailed ways

Embodiment one

Fig. 1 shows the specific flow of the method of the present invention. The processing process of the method of the present invention is described in detail below in conjunction with a set of publicly available lung 4D-CT data sets. The dataset consists of 10 groups of lung 4D-CT data, each containing 10 phase images. The present invention is based on fast sub-pixel motion estimation lung 4D-CT image super-resolution reconstruction method, the specific steps are as follows:

(1) Read the initial 4D-CT image of the lungs. This data is selected from the second group of the public data set. The image size is 256*256*112, the resolution within the image layer is 1.16mm, and the resolution between layers is 2.5mm , the lung 4D-CT image is composed of 10 different phase lung 3D-CT images, and the phase 0 coronal sagittal plane low-resolution lung 3D-CT image is selected as the lung 3D-CT image to be reconstructed, and the selected The lung 3D-CT image is randomly selected, which phase of the lung 3D-CT image is selected, then the high-resolution lung 4D-CT image reconstructed in the subsequent step (4) is exactly the same as the lung 3D-CT image to be reconstructed. CT images correspond to lung 4D-CT images with the same phase;

(2) Perform initial motion estimation on multiple lung 3D-CT images with different phases after removing the lung 3D-CT image to be reconstructed in the lung 4D-CT image relative to the lung 3D-CT image to be reconstructed, and obtain the lung The initial motion vector field between multiple different phase lung 3D-CT images after removing the lung 3D-CT image to be reconstructed and the lung 3D-CT image to be reconstructed in the first 4D-CT image, the initial motion vector field The precision is an integer;

(3) Carry out precision optimization to the obtained initial motion vector field, make it accurate to sub-pixel, obtain sub-pixel motion vector field;

(4) Based on the sub-pixel motion vector field obtained in step (3), the 3D-CT image of the lung to be reconstructed is reconstructed to obtain a reconstructed high-resolution image with the same phase as the 3D-CT image of the lung to be reconstructed. Rate lung 4D-CT images.

In the present invention, in the step (2), a three-step search method is used for initial motion estimation, and multiple lung 3D-CT images of different phases after removing the lung 3D-CT image to be reconstructed are searched for in the lung 4D-CT image. The integer pixel displacement between the image and the 3D-CT image of the lung to be reconstructed. The algorithm uses a coarse-to-fine search pattern, starting from the center point of the search window, taking 8 surrounding points at a certain step to form a point group for each search, and then performing matching calculations according to the matching criteria to find the center of the matching block with the smallest error point. The specific process includes:

(2.1) Determine a center point, set the maximum search length, and use 1/2 of the maximum search length as the step size, and find the minimum block error point according to the matching criteria for the center point and 8 detection points with the same step length around the distance, If the minimum block error point is located at the original center point, the algorithm ends, otherwise, proceed to step (2.2);

(2.2) The step size is halved, and the center point of the minimum error matching block is found in the minimum block error point determined in the previous step and 8 detection points of the same step length around it;

(2.3) Repeat (2.1) and (2.2) until the step size meets the search accuracy requirement, that is, the best matching point is obtained.

What the present invention adopts is that the minimum absolute error matching criterion commonly used is defined as:

$S(v_x,v_{\mathrm{the\; y}})={\mathrm{\Σ}}_{i=0}^{N-1}{\mathrm{\Σ}}_{j=0}^{N-1}|I_t(p+i,q+j)-I_{t-{\mathrm{\Δ}}_t}(p+v_x+i,q+v_{\mathrm{the\; y}}+j)|$ ......Formula 1)

In the formula, the size of the image block is N×N, the coordinates of the upper left corner are (p, q), and the block size selected by the present invention is 16×16; the motion vector is (v _x , v _y ); I _t (i, j) and are the values at pixel (i, j) of the current frame and the reference frame, respectively. The motion vector with the smallest value of S(v _x , v _y ) in the search area is the optimal motion vector of the current block.

In the present invention, in the step (3), the optical flow method is used to optimize the accuracy of the initial motion vector to make it accurate to sub-pixels. The specific process includes:

(3.1) Assuming that the two consecutive images before and after motion are f(x,y) and g(x,y) respectively, according to the simplified model of classical optical flow, we can get

$g(x,\mathrm{the\; y})=f(x+\mathrm{\Δ}{x}_{\mathrm{the\; s}},\mathrm{the\; y}+\mathrm{\Δ}{\mathrm{the\; y}}_{\mathrm{the\; s}})\≈f(x,\mathrm{the\; y})+\mathrm{\Δ}{x}_{\mathrm{the\; s}}\frac{\∂}{\∂x}f(x,\mathrm{the\; y})+\mathrm{\Δ}{\mathrm{the\; y}}_{\mathrm{the\; s}}\mathrm{\Δ}{x}_{\mathrm{the\; s}}\frac{\∂}{\∂\mathrm{the\; y}}f(x,\mathrm{the\; y})$ ...Formula (2)

This formula is approximately a first-order Taylor series expansion.

(3.2) Minimize and solve formula (2) to obtain the optimal displacement vector:

$\underset{\mathrm{\Δ}{x}_{\mathrm{the\; s}},\mathrm{\Δ}{\mathrm{the\; y}}_{\mathrm{the\; s}}}{\min \mathrm{imize}}\mathrm{\Φ}(\mathrm{\Δ}{x}_{\mathrm{the\; s}},\mathrm{\Δ}{\mathrm{the\; y}}_{\mathrm{the\; s}})$ ... Formula (3);

in:

$\mathrm{\Φ}(\mathrm{\Δ}{x}_{\mathrm{the\; s}},\mathrm{\Δ}{\mathrm{the\; y}}_{\mathrm{the\; s}})=\underset{x,\mathrm{the\; y}}{\mathrm{\Σ}}{(g(x,\mathrm{the\; y})-f(x,\mathrm{the\; y})-\mathrm{\Δ}{x}_{\mathrm{the\; s}}\frac{\∂}{\∂x}f(x,\mathrm{the\; y})-\mathrm{\Δ}{\mathrm{the\; y}}_{\mathrm{the\; s}}\frac{\∂}{\∂\mathrm{the\; y}}f(x,\mathrm{the\; y}))}^2$ ...Formula (4);

(3.3) Equation (4) can be regarded as a linear least squares solution problem, and the partial derivatives of the objective function can be set to 0 to obtain the optimal values Δx _s , Δy _s . Therefore, the following equation can be obtained:

$\frac{\∂\mathrm{\Φ}}{\∂\mathrm{\Δ}{x}_{\mathrm{the\; s}}}=0$ and $\frac{\∂\mathrm{\Φ}}{\∂\mathrm{\Δ}{\mathrm{the\; y}}_{\mathrm{the\; s}}}=0$ ...Formula (5);

(3.4) Solving equation (5), the optimal solution Δx _s , Δy _s can be obtained.

(3.5) Suppose the motion vector obtained by the three-step search method is Δx _I , Δy _I , then the final total motion vector Δx, Δy is:

Δx = Δx _I + Δx _s , Δy = Δy _I + Δy _s ...... formula (6);

In this example, the step (4) uses an iterative back projection method to reconstruct a high-resolution lung 4D-CT image. The abbreviation of the iterative back projection method is IBP. The specific process includes:

(4.1) interpolating and enlarging the low-resolution lung 3D-CT image to be reconstructed as the initial high-resolution image H ⁽⁰⁾ ;

(4.2) Simulate the imaging process according to the degradation model to obtain a collection of low-resolution images corresponding to the original low-resolution image sequence K represents the number of images in the original low-resolution image sequence, and n is the number of iterations.

Specifically, in the nth iteration process, H ⁽ⁿ⁾ simulation degenerates to The process is as follows:

$L_k^{\left(\mathrm{no}\right)}=\left(T_k\left(h^{\left(\mathrm{no}\right)}\right)h\right)\↓\mathrm{the\; s}$ ...Formula (7);

Among them, T _k represents the two-dimensional geometric transformation from H to L _k , which is the motion deformation field obtained in step (3); h is the Gaussian blur operator; ↓s is the downsampling operator;

(4.3) Judgment error Whether it reaches the minimum value, if reached, then stop the iteration, the previously estimated H ⁽ⁿ⁾ is the final super-resolution image; if not reached, then enter step (4.4);

Judgment error in the above step (4.3) Whether the minimum value is reached can be determined by judging whether the error function e ⁽ⁿ⁾ is smaller than the set threshold ε. The specific calculation formula of the error function is:

$e^{\left(\mathrm{no}\right)}=\sqrt{\frac{1}{K}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\left|\right|L_k-L_k^{\left(\mathrm{no}\right)}\left|\right|}_2^2}$ ...Formula (8);

(4.4) Update the current high-resolution image according to the error, and the update process is as follows:

$h^{(\mathrm{no}+1)}=h^{\left(\mathrm{no}\right)}+\frac{1}{K}{T}_k^{-1}\left(((L_k-L_k^{\left(\mathrm{no}\right)})\&Up\; Arrow;\mathrm{the\; s})p\right)$ ... Formula (9);

In the formula, ↑s represents the up-sampling operator; p represents the back-projection operator, which depends on h and T _k ;

(4.5) Use the updated high-resolution image as the initial high-resolution image, and enter step (4.2).

Phase 0 coronal low-resolution image reconstruction results are shown in Figure 2 and Figure 3. From left to right are the results reconstructed using bilinear interpolation method, based on the full search motion estimation algorithm and the method of the present invention; FIG. 3 is an enlarged schematic diagram corresponding to the box part in FIG. 2 .

Phase 0 sagittal low-resolution image reconstruction results are shown in Figure 4 and Figure 5. From left to right are the reconstruction results based on the bilinear interpolation method, the full search motion estimation algorithm and the method of the present invention; FIG. 5 is an enlarged schematic diagram corresponding to the box part in FIG. 4 .

It is worth noting that the reconstruction result of the method of the present invention is better than that of the reconstruction algorithm based on full search, which not only enhances the edges and details of the image, but especially reduces the error caused by motion estimation. As indicated by the white arrows in the second column of Fig. 2 and Fig. 4, since the full search motion estimation algorithm uses a fixed step size, the matching error of some image blocks will be large, and the error will be superimposed during the reconstruction process. The results are shown as highlighted areas. However, the motion estimation algorithm adopted by the method of the present invention can realize more accurate sub-pixel image block matching, so there will be no highlighted areas caused by error iterations.

In addition, this example also objectively evaluates the effectiveness of the present invention through quantitative indicators. In this example, the image average gradient is used to evaluate the phase 0 coronal sagittal plane reconstruction results.

The average gradient not only reflects the micro-detail contrast and texture change characteristics in the image, but also reflects the sharpness of the image. The larger the average gradient, the clearer the image, the better the contrast, and the better the details are preserved. It is defined as:

$\▿\stackrel{\‾}{f}=\frac{1}{(m-1)(N-1)}\×{\mathrm{\Σ}}_{i=1}^{m-1}{\mathrm{\Σ}}_{j=1}^{N-1}\sqrt{\frac{{\▿}_i^{2}f(i,j)+{\▿}_j^2(i,j)}{2}}$ ...Formula (10)

In the formula, f(i, j), ▽ _i f(i, j) and ▽ f _j (i, j) are the pixel gray level and its gradient in the row and column direction respectively; M and N are the image number of rows and columns.

According to Equation (10), the average gradient is calculated for the phase 0 coronal sagittal high-resolution images reconstructed based on the full search motion estimation algorithm and the method in this paper, respectively, using bilinear interpolation. The results are shown in Table 1. Table 1 is a comparison table of the average gradient of the coronal sagittal plane reconstruction results shown in Figure 2 and Figure 4 at phase 0. It can be seen that the average gradient value of the algorithm in this paper is significantly larger than that of the bilinear interpolation method, and compared with the reconstruction results based on full search, it also has a certain improvement, the image is clearer, and the details are better preserved.

Table 1: Comparison table of average gradient of coronal sagittal plane reconstruction results in phase 0 shown in Figure 2 and Figure 4

data bilinear interpolation full search algorithm The method of the invention phase 0 coronal plane 5.67 8.97 9.14

Phase 0 sagittal plane 5.18 7.68 8.45 average value 5.43 8.33 8.80

Finally, in this example, the time-consuming reconstruction results obtained by using different algorithms for the same set of data were counted, and the results are shown in Table 2. Table 2 is a time-consuming comparison table of reconstruction results obtained by different methods for the same data. It can be seen that due to the optimization of motion vector estimation speed and precision, the time consumption of the method of the present invention is reduced by nearly 20 times compared with the full search algorithm, which is a significant improvement.

Table 2: Time-consuming comparison table of reconstruction results obtained by different methods for the same data

algorithm full search algorithm The method of the invention time consuming 40min 2min

Integrating all the contents in Embodiment 1 of the present invention: Fig. 2, Fig. 3, Fig. 4, Fig. 5 show from the visual effect that the brightness and clarity of the blood vessels and surrounding tissues in the lung parenchyma of the high-resolution image reconstructed by the method of the present invention are uniform Significantly enhanced without highlighting areas caused by erroneous iterations. Table 1 further objectively evaluates the effectiveness of the method of the present invention from the quantitative indicators. Most importantly, Table 2 shows that the method of the present invention greatly reduces the time consumption, which is a significant improvement.

It can be seen that compared with the prior art, the present invention can obtain a clearer lung image, and the image structure is obviously enhanced. At the same time, compared with the full search algorithm, the operation time is significantly reduced.

Finally, it should be noted that the embodiments of the present invention are not limited thereto, and may be modified according to actual needs to meet different actual needs. Therefore, under the premise of the above-mentioned basic technical idea of the present invention, other modifications, replacements or changes made to the content of the present invention in accordance with ordinary technical knowledge and conventional means in this field all fall within the protection scope of the present invention.

Claims (4)

1., based on the lung 4D-CT image super-resolution rebuilding method of fast sub-picture element estimation, the method comprises the following steps:

(1) read initial lung 4D-CT image, this lung 4D-CT image is made up of lung's 3D-CT image of multiple out of phase, selects arbitrarily the lung 3D-CT image of wherein a certain phase place as lung 3D-CT image to be reconstructed;

(2) the lung 3D-CT image of the multiple outs of phase after removing lung 3D-CT image to be reconstructed in lung 4D-CT image is carried out initial motion estimation relative to lung 3D-CT image to be reconstructed, obtain the initial motion vectors field between the lung 3D-CT image of the multiple outs of phase after removing lung 3D-CT image to be reconstructed in lung 4D-CT image and lung 3D-CT image to be reconstructed, the precision of this initial motion vectors field is integer;

(3) precision optimizing is carried out to the initial motion vectors field obtained, make it be accurate to sub-pix, obtain sub-pel motion vector field;

(4), based on the sub-pel motion vector field obtained by step (3), lung 3D-CT image to be reconstructed is rebuild, obtains the high resolving power lung 4D-CT image after the identical reconstruction of the phase place corresponding with lung 3D-CT image to be reconstructed.

2. the lung 4D-CT image super-resolution rebuilding method based on fast sub-picture element estimation according to claim 1, it is characterized in that: in described step (2), adopt three step search algorithm to carry out initial motion estimation, the integer pixel displacement between the lung 3D-CT image of the multiple outs of phase after lung 3D-CT image to be reconstructed and lung 3D-CT image to be reconstructed is removed in search lung 4D-CT image, this algorithm passes through by the thick search pattern to essence, from search window centre point, get 8 points around by a fixed step size and form each point group searched for, then matching primitives is carried out according to matching criterior, the match block central point finding error minimum, detailed process comprises:

(2.1) central point is determined, setting maximum search length, using 1/2 of maximum search length as step-length, by 8 check points of central point and the around identical step-length of distance according to matching criterior, find smallest blocks error point, if smallest blocks error point is positioned at former central point, then algorithm terminates, otherwise, carry out step (2.2);

(2.2) step-length reduces by half, the smallest blocks error point determined in previous step and around identical step-length 8 check points in find the central point of least error match block;

(2.3) repeat (2.1) and (2.2) until step-length reaches search precision requirement, namely obtain optimal match point;

What the present invention adopted is conventional least absolute error matching criterior, is defined as:

$S(v_x,v_y)={\mathrm{\Σ}}_{i=0}^{N-1}{\mathrm{\Σ}}_{j=0}^{N-1}|I_t(p+i,q+j)-I_{t-{\mathrm{\Δ}}_t}(p+v_x+i,q+v_y+j)|$ ... formula (1)

In formula, the size of image block is N × N, and the coordinate in the upper left corner is (p, q), and the block size that the present invention chooses is 16 × 16; Motion vector is (v _x, v _y); I _t(i, j) and be respectively present frame and the value of reference frame at pixel (i, j) place, in region of search, make S (v _x, v _y) the minimum motion vector of value is the optimal motion vector of current block.

3. the lung 4D-CT image super-resolution rebuilding method based on fast sub-picture element estimation according to claim 1, it is characterized in that: in described step (3), utilize optical flow method to carry out precision optimizing to initial motion vectors, make it be accurate to sub-pix, detailed process comprises:

(3.1) suppose that motion front and back two frame consecutive images are respectively f (x, y) and g (x, y), according to classical light stream simplified model, can obtain

$g(x,y)=f(x+\mathrm{\Δ}{x}_s,y+\mathrm{\Δ}{y}_s)\≈f(x,y)+\mathrm{\Δ}{x}_s\frac{\∂}{\∂x}f(x,y)+\mathrm{\Δ}{y}_s\mathrm{\Δ}{x}_s\frac{\∂}{\∂y}f(x,y)$ ... formula (2)

This formula is approximately first order Taylor series expansion;

(3.2) minimize to solve to formula (2) optimum displacement vector can be obtained:

$\underset{\mathrm{\Δ}{x}_s,\mathrm{\Δ}{y}_s}{\min \mathrm{imize}}\mathrm{\Φ}(\mathrm{\Δ}{x}_s,\mathrm{\Δ}{y}_s)$ ... formula (3);

Wherein:

$\mathrm{\Φ}(\mathrm{\Δ}{x}_s,\mathrm{\Δ}{y}_s)=\underset{x,y}{\mathrm{\Σ}}{(g(x,y)-f(x,y)-\mathrm{\Δ}{x}_s\frac{\∂}{\∂x}f(x,y)-\mathrm{\Δ}{y}_s\frac{\∂}{\∂y}f(x,y))}^2$ ... formula (4);

(3.3) formula (4) can regard linear least-squares Solve problems as, the partial derivative of objective function can be set to 0 and obtain optimal value Δ x _s, Δ y _s, therefore, following equation can be had:

$\frac{\∂\mathrm{\Φ}}{\∂\mathrm{\Δ}{x}_s}=0$ With $\frac{\∂\mathrm{\Φ}}{\∂\mathrm{\Δ}{y}_s}=0$ ... formula (5);

(3.4) solving equation (5), can obtain optimum solution Δ x _s, Δ y _s;

(3.5) motion vector setting three step search algorithm to obtain is as Δ x _i, Δ y _i, then the total motion vector Δ x finally obtained, Δ y is:

Δ x=Δ x _i+ Δ x _s, Δ y=Δ y _i+ Δ y _s... formula (6).

4. the lung 4D-CT image super-resolution rebuilding method based on fast sub-picture element estimation according to claim 1, it is characterized in that: described step (4) adopts iterative backprojection method to rebuild high resolving power lung 4D-CT image, iterative backprojection method be abbreviated as IBP, detailed process comprises:

(4.1) by the lung 3D-CT image interpolate enlarge to be reconstructed of low resolution, as initial high-resolution image H ⁽⁰⁾;

(4.2) set of a low-resolution image is obtained according to degradation model analog imaging process correspond respectively to original low-resolution image sequence k represents the quantity of image in original low-resolution image sequence, and n is iterations;

Concrete, in n-th iterative process, H ⁽ⁿ⁾simulation deteriorates to process represent as follows:

$L_k^{\left(n\right)}=\left(T_k\left(H^{\left(n\right)}\right)h\right)\↓s$ ... formula (7);

Wherein, T _krepresent from H to L _ktwo-dimensional geometry conversion, be the motion deformation field obtained in step (3); H is Gaussian Blur operator; ↓ _sit is down-sampling operator;

(4.3) error in judgement whether reach minimum value, if reach, then stop iteration, the H estimated in the past ⁽ⁿ⁾for final required super-resolution image; If do not reach, then enter step (4.4);

Error in judgement in above-mentioned steps (4.3) whether reaching minimum value can by error in judgement function e ⁽ⁿ⁾whether be less than setting threshold epsilon to carry out, the specific formula for calculation of error function is:

$e^{\left(n\right)}=\sqrt{\frac{1}{K}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\left|\right|L_k-L_k^{\left(n\right)}\left|\right|}_2^2}$ ... formula (8);

(4.4) according to error, Current high resolution image is upgraded, renewal process as shown in the formula:

$H^{(n+1)}=H^{\left(n\right)}+\frac{1}{K}{T}_k^{-1}\left(((L_k-L_k^{\left(n\right)})\↑s)p\right)$ ... formula (9);

In formula, ↑ s represents up-sampling operator; P represents back projection operator, depends on h and T _k;

(4.5) using upgrade after high-definition picture as initial high-resolution image, enter step (4.2).