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. The method comprises the steps that (1) an initial lung 4D-CT image is read, the lung 4D-CT image is composed of a plurality of lung 3D-CT images of different phases, and the lung 3D-CT image of any phase is selected to serve as the lung 3D-CT image to be reconstructed; (2) initial motion estimation is carried out on lung 3D-CT images, expect 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 an initial motion vector field among the lung 3D-CT images is obtained; (3) precision optimization is carried out on the obtained initial motion vector field to obtain a sub-pixel motion vector field; (4) on the basis of the sub-pixel motion field, the lung 3D-CT image to be reconstructed is reconstructed, and a reconstructed high-resolution lung 4D-CT image which corresponds to the lung 3D-CT image to be reconstructed and is of the same phase as the lung 3D-CT image to be reconstructed is obtained. By the adoption of the lung 4D-CT image super-resolution reconstruction method, the lung 4D-CT image resolution can be improved.

Description

Based on the lung 4D-CT image super-resolution rebuilding method of fast sub-picture element estimation

Technical field

The present invention relates to technical field of medical image processing, specifically refer to a kind of lung 4D-CT image super-resolution rebuilding method based on fast sub-picture element estimation.

Background technology

Lung 4D-CT image plays an important role in radiotherapy of lung cancer, and it can provide a comprehensive high accuracy radiation therapy respiratory movement and characterize, and contributes to tracking of knub motion, implements precise radiotherapy, and reduce the damage of normal tissue.But, due to the high-dose irradiation that CT is intrinsic, therefore the sampling that often can only reduce longitudinally (Z-direction) is to reduce lung 4D-CT sweep time in the hope of reducing radiant quantity, thus resolution between lung 4D-CT image layer is caused far below layer intrinsic resolution, to cause the significant anisotropy of data.When this makes to carry out many viewed in plan (hat sagittal plane etc.) to data, need carry out interpolation arithmetic to obtain correct display, this operation easily causes the fuzzy of image.

In super-resolution rebuilding process, between image, the estimation of sports ground is the principal element affecting reconstruction precision and speed.The super-resolution technique based on estimation proposed before applicant have employed the full searching moving estimation technique, estimates the sports ground between different frame image.The method can reconstruct more traditional interpolation method lung hat (arrow) shape face image more clearly, but its major defect is: speed is slow, and can only be fixing search step-length.So both affect reconstruction speed, also affect reconstruction precision.Therefore, provide a kind of lung 4D-CT image super-resolution rebuilding method based on fast sub-picture element estimation very necessary to overcome the deficiencies in the prior art.

Summary of the invention

The object of the present invention is to provide a kind of lung 4D-CT image super-resolution rebuilding method based on fast sub-picture element estimation, the method can improve lung 4D-CT image resolution ratio.

Object of the present invention realizes by following technique measures: 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.

In the present invention, adopt three step search algorithm to carry out initial motion estimation in described step (2), in search lung 4D-CT image, remove 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.This algorithm, by by slightly to the search pattern of essence, from search window centre point, is got the point group of 8 each search of some formation around, is then carried out matching primitives according to matching criterior, the match block central point finding error minimum by a fixed step size.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 that 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.S (v is made in region of search _x, v _y) the minimum motion vector of value is the optimal motion vector of current block.

In the present invention, utilize optical flow method to carry out precision optimizing to initial motion vectors in described step (3), 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);

In the present invention, 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).

Compared to the prior art, the present invention has following beneficial effect:

(1) due to the optimization of estimation of motion vectors speed and precision, computing time of the inventive method, comparatively full-search algorithm reduced nearly 20 times, was improvement highly significant.

(2) the high resolving power lung 4D-CT image rebuild of the present invention, the blood vessel in its pulmonary parenchyma and the brightness of perienchyma and sharpness all significantly strengthen.

Accompanying drawing explanation

Describe the present invention below in conjunction with the drawings and specific embodiments.

Fig. 1 is the process flow diagram of the lung 4D-CT image super-resolution rebuilding method that the present invention is based on fast sub-picture element estimation;

Fig. 2 is that the embodiment of the present invention one data 2 phase place 0 coronal-plane adopts distinct methods reconstructed results figure, is respectively from left to right to adopt bilinear interpolation method, based on the result figure that full searching moving algorithm for estimating and the inventive method are rebuild;

Fig. 3 is the enlarged diagram that the embodiment of the present invention one corresponds to Blocked portion in Fig. 2, is adopt bilinear interpolation method from left to right respectively, based on the partial enlarged drawing of the result figure that full searching moving algorithm for estimating and the inventive method are rebuild;

Fig. 4 is that the embodiment of the present invention one data 2 phase place 0 sagittal plane adopts distinct methods reconstructed results figure, is respectively from left to right to adopt bilinear interpolation method, based on the result figure that full searching moving algorithm for estimating and the inventive method are rebuild;

Fig. 5 is the enlarged diagram that the embodiment of the present invention one corresponds to Blocked portion in Fig. 4, is adopt bilinear interpolation method from left to right respectively, based on the partial enlarged drawing of the result figure that full searching moving algorithm for estimating and the inventive method are rebuild.

Embodiment

Embodiment one

Fig. 1 shows the idiographic flow of the inventive method.The processing procedure of the inventive method is described in detail below in conjunction with a set of publicly available lung 4D-CT data set.This data set is made up of 10 groups of lung 4D-CT data, often organizes packet containing 10 phase images.The present invention is based on the lung 4D-CT image super-resolution rebuilding method of fast sub-picture element estimation, concrete steps are as follows:

(1) initial lung 4D-CT image is read, these data are selected from the 2nd group of common data sets, image size is 256*256*112, image layer intrinsic resolution is 1.16mm, interlayer is distinguished as 2.5mm, this lung 4D-CT image is made up of lung's 3D-CT image of 10 outs of phase, selected phase 0 is preced with sagittal plane low resolution lung 3D-CT image as lung 3D-CT image to be reconstructed, the lung 3D-CT image chosen is what choose arbitrarily, choose the lung 3D-CT image of which phase place, then rebuilding the high resolving power lung 4D-CT image obtained in subsequent step (4) is exactly the identical lung 4D-CT image of the corresponding phase place of lung 3D-CT image to be reconstructed with this,

(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.

In the present invention, adopt three step search algorithm to carry out initial motion estimation in described step (2), in search lung 4D-CT image, remove 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.This algorithm, by by slightly to the search pattern of essence, from search window centre point, is got the point group of 8 each search of some formation around, is then carried out matching primitives according to matching criterior, the match block central point finding error minimum by a fixed step size.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 that 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.S (v is made in region of search _x, v _y) the minimum motion vector of value is the optimal motion vector of current block.

In the present invention, utilize optical flow method to carry out precision optimizing to initial motion vectors in described step (3), 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);

In this example, 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; ↓ s 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).

Phase place 0 coronal-plane low-resolution image reconstructed results is as shown in Fig. 2, Fig. 3.Adopt bilinear interpolation method from left to right respectively, based on the result figure that full searching moving algorithm for estimating and the inventive method are rebuild; Fig. 3 is the enlarged diagram corresponding to Blocked portion in Fig. 2.

Phase place 0 sagittal plane low-resolution image reconstructed results is as shown in Fig. 4, Fig. 5.Adopt bilinear interpolation method from left to right respectively, based on the result figure that full searching moving algorithm for estimating and the inventive method are rebuild; Fig. 5 is the enlarged diagram corresponding to Blocked portion in Fig. 4.

What deserves to be explained is, the inventive method reconstructed results is comparatively better based on full search reconstruction algorithm, not only makes image border and details obtain enhancing, especially reduces the error that estimation causes.As the 2nd row white arrow instruction in Fig. 2, Fig. 4, due to full searching moving algorithm for estimating, use fixed step size, the matching error of parts of images block can be caused comparatively large, and cause the superposition of process of reconstruction medial error, finally be shown as highlight regions in the result.And the motion estimation algorithm that the inventive method adopts can realize more accurate subpixel image as Block-matching, therefore there will not be the highlight regions caused by error iteration.

In addition, this example also index objective evaluation validity of the present invention by quantifying.This example adopts image averaging Gradient phase place 0 to be preced with sagittal reconstruction result.

Average gradient had both reflected minor detail contrast in image and texture variations feature, also reflects the sharpness of image.Average gradient is larger, show that image is more clear, contrast better, Hemifusus ternatanus better.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 formula, f (i, j), ▽ _if (i, j) and ▽ f _j(i, j) is respectively pixel gray scale and it is expert at, the gradient on column direction; M and N is respectively line number and the columns of image.

According to formula (10) to utilizing bilinear interpolation, the phase place 0 reconstructed based on full searching moving algorithm for estimating and context of methods is preced with sagittal plane high-definition picture, and be averaged gradient calculation respectively, and result is as shown in table 1.The hat sagittal reconstruction result images average gradient contrast table that table 1 is phase place 0 shown in Fig. 2 and Fig. 4.Visible, herein algorithm enlarges markedly compared with the average gradient value of bilinear interpolation method, and with based on entirely searching for compared with reconstructed results, also improve, image is more clear, Hemifusus ternatanus better.

Table 1: the hat sagittal reconstruction result images average gradient contrast table of phase place 0 shown in Fig. 2 and Fig. 4

Data	Bilinear interpolation	Full-search algorithm	The inventive method
				Phase place 0 coronal-plane	5.67	8.97	9.14

Phase place 0 sagittal plane	5.18	7.68	8.45
				Mean value	5.43	8.33	8.80

Finally, added up respectively in this example and adopted algorithms of different to obtain the consuming time of reconstructed results for same group of data, result is as shown in table 2.Table 2 is for adopting distinct methods to obtain the contrast table consuming time of reconstructed results for identical data.Visible, due to the optimization of estimation of motion vectors speed and precision, the comparatively full-search algorithm consuming time of the inventive method reduces nearly 20 times, is improvement highly significant.

Table 2: adopt distinct methods to obtain the contrast table consuming time of reconstructed results for identical data

Algorithm	Full-search algorithm	The inventive method
			Consuming time	40min	2min

All the elements in cumulated volume inventive embodiments one: Fig. 2, Fig. 3, brightness and the sharpness of blood vessel from visual effect shows its pulmonary parenchyma of high-definition picture that the inventive method is rebuild of Fig. 4, Fig. 5 and perienchyma all significantly strengthen, and there will not be the highlight regions caused by error iteration.Table 1 is again further from the validity of objective evaluation the inventive method quantizating index.Most importantly, table 2 shows that the inventive method is significantly reduced on consuming time, is the very significant improvement of tool.

Visible, compared with prior art, the present invention can obtain lung images more clearly, and picture structure obviously strengthens.Meanwhile, compared with full-search algorithm, operation time significantly reduces.

Finally it should be noted that embodiments of the present invention are not limited thereto, can modify according to actual needs, to adapt to different actual demands.Therefore, under stating basic fundamental thought prerequisite on the invention, according to the ordinary technical knowledge of this area and customary means to content of the present invention make the amendment of other various ways, replacement or change, all drop within rights 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).