CN104268914B - Reestablishing method of 4D-CT (Four Dimensional-Computed Tomography) different time phase sequence image - Google Patents

Abstract

The invention discloses a reestablishing method of a 4D-CT different time phase sequence image. The reestablishing method comprises the following steps: establishing an edge-protected multi-scale space, conducting multi-scale decomposition on a 3D-CT image on a therapeutic day under a free breathing state, determining an optimal time phase based on a maximal similarity measurement criterion, utilizing a free deformation model based on a spline B to realize deformation rectification of a 4D-CT different time phase image on a planned day, and reestablishing the 4D-CT image based on the 4D-CT image on the planned day and the 3D-CT image on the therapeutic day. The reestablishing method of the 4D-CT different time phase sequence image, disclosed by the invention, has the benefits that the 4D-CT image on the therapeutic day can be reestablished based on the 4D-CT image on the planned day and the 3D-CT image on the therapeutic day, so as to acquire the motion change law of tumor and endangered organ of a patient on the therapeutic day, and effectively avoid the damage to the body as the patient receives plenty of X-ray radiation again.

Description

A kind of method for reconstructing of 4d-ct difference phase sequence image

Technical field

The present invention relates to image processing field, more particularly, to a kind of method for reconstructing of 4d-ct difference phase sequence image.

Background technology

Clinical 4d-ct image brings larger radiation, x-ray to patient in implementation process, and tumour and soft group dynamic to swollen Knurl and the impact jeopardizing organ, but based on plan day, (before cancer patient's treatment, physics teacher and clinician utilize its 4d- for patient The radiotherapy planning that ct view data is done) 4d-ct radiotherapy planning distance treatment day (in cancer patient's Patients During Radiotherapy) view data Obtain the time interval in about 1 to 2 week, tumour or other tissue it may happen that large change.If it is heavy in each treatment day New scanning 4d-ct, patient will accept a large amount of radiation, x-ray again.

At present in treatment day, easily obtain 3d-ct or cbct (cone beam ct) image, these images are in patient certainly By obtain under breathing state.Although 4d-ct technology can effectively analyze respiratory movement to tumor target and the shadow jeopardizing organ Ring, but 4d-ct often obtains in plan day, and the tumour of patient or other physical trait may occur larger change in treatment day Change, existing method is that 3d-ct the or cbct image obtaining under the free breathing state of patient in treatment day carries out instructing pendulum position, Because patient respiratory motion or pendulum position etc. cause the motion of tumour irradiated site to cause guidance pendulum position to have error, and then affect to swollen The border of knurl and effective judgement of the characteristics of motion.

Content of the invention

The purpose of the present invention be exactly in order to solve the above problems it is proposed that a kind of 4d-ct difference phase sequence image weight The 3d-ct image of construction method, the 4d-ct of the method application plan day and treatment day to rebuild the 4d-ct image for the treatment of day, thus Obtain patient and treat day tumour and the motion change rule jeopardizing organ, realize to tumour and the motion change rule jeopardizing organ And the quantitative analysis on border.

To achieve these goals, the present invention adopts the following technical scheme that

A kind of method for reconstructing of 4d-ct difference phase sequence image, comprises the following steps:

The first step: remove plan day 4d-ct image therapeutic bed artifact, and remove noise in image or other shake into Point；

Second step: the multi-resolution decomposition based on the edge-protected multiscale space of tv-l1 is carried out to image, by using difference Smoothing parameter λ, obtain different respiratory phase comprises the how much different smoothed image of information content；

3rd step: joint spatial-temporal domain information deformable registration is carried out to the smoothed image of above-mentioned difference respiratory phase, by deformation Smoothed image after registration is defined as target image；

4th step: remove the noise of 3d-ct image under the free breathing state obtaining treatment day or other shake composition；

5th step: multi-resolution decomposition is carried out to 3d-ct image by the edge-protected multiscale space of tv-l1, and passes through Similarity measure determines optimum phase, is labeled as incomplete treatment day 4d-ct, and incomplete treatment day 4d-ct is defined For reference picture；

6th step: respectively the target image of different phases is registrated to reference using the free deformation model based on b batten In corresponding phase in image；Image after registration terminates is the treatment day 4d-ct image rebuild.

In described 3rd step, joint spatial-temporal domain information deformable registration is carried out to a series of smoothed image of difference respiratory phase Method particularly includes:

Described smoothed image to be represented with intensity function f (i, k), wherein, i and k is n (typically taking 10) respectively The spatially and temporally sample of individual respiratory phase；Described in sequence image time domain using the cycle Model in Time Domain based on b batten Correlation, the movement locus of tumour is modeled as the time domain smooth function based on b batten；Consider air-breathing end in respiratory movement simultaneously Begin to exhaling the snap back conversion in moment organ movement direction, set up segmentation constraint function it is ensured that whole respiratory movement track While smooth, keep breathing track in the accuracy of air-breathing end movement locus；

To express the deformation domain of the respiratory movement correlated series image of 4d-ct difference phase using dual-tree complex wavelet transform, Deformation domain is modeled as the function of wavelet coefficient；Set up deformation energy function by navier partial differential equation, by coarse and fine The optimal wavelet coefficient of representation space domain deformation characteristics is estimated on ground, and described wavelet coefficient is substituted into target image transformation model to mesh Logo image carries out deformation conversion.Final cycle Model in Time Domain (smooth track model and the piecewise smooth rail utilizing based on b batten Mark model) make with a kind of method of overall situation combining time domain and spatial information (si) of dual-tree complex wavelet transform deformation model coupling foundation For sample data.

Described the movement locus of tumour is modeled as time domain smooth function based on b batten particularly as follows:

$t_t(x,t)=x+\underset{l\&element;l}{\mathrm{\σ}}{b}_l{\mathrm{\ψ}}_l\left(t\right)$

Wherein, τ_t(x, t) represents t_yMoment when phase images in reference time t_yAnd the movement locus at the x of position, it is to close In the continuously smooth function of time t, t_yRepresent the y moment of t, ψ_l(t)=β^m(t/s-l), s ∈ r is the control point interval of time； b_l∈r³Coefficient for m rank b spline base function；β^mFor m rank b spline base function, β is b spline base function node, and l is b spline Basis Function counting variable；

The overall smooth track model of b batten cycle Model in Time Domain is:

$\frac{{\&partiald;}^z{t}_t(x,0)}{{\&partiald;t}^z}=\frac{{\&partiald;}^z{t}_t(x,t_e)}{{\&partiald;t}^z},z=[0,...,m-1]$

Wherein, t_tMovement locus at x for phase images when (x, 0) represented for 0 moment；t_t(x,t_e) represent t the e moment when Movement locus at x for the phase images.

The segmentation constraint function of described foundation particularly as follows:

$t_t^{*}(x,0)=t_t^{*}(x,t_e)$

Wherein,Represent sectionally smooth track at x for phase images during 0 moment,Represent the e moment of t When sectionally smooth track at x for the phase images；

For sectionally smooth track it is assumed that air-breathing foot couple answers t=0, single constraints can be set and apply at air-breathing end.

Described by deformation domain be modeled as wavelet coefficient function particularly as follows:

X'=x+u¹(x,y,z；c)

Y'=y+u²(x,y,z；c)

Z'=z+u³(x,y,z；c)

Wherein, (x, y, z) represents the space coordinates of reference picture, and (x ', y ', z ') represents the space coordinates of target image System, deformation domain u=(u¹,u²,u³), c is wavelet coefficient, and deformation domain u is the function of wavelet coefficient c.

Described set up deformation energy function by navier partial differential equation particularly as follows:

E (c)=inter (c)+w exter (c)

Wherein, w is weighting constant, adopts w to be constant 1, inter (c) represents internal force constraint function, exter (c) in experiment Represent external force constraint function.

In described 5th step, determine that by similarity measure the method for the optimum phase of 3d-ct image is: 3d-ct is schemed As and plan day 4d-ct in n when phase images carry out similarity measure assessment, accordingly select with to plan day 4d-ct closest N when phase images.

The invention has the beneficial effects as follows:

The present invention not only possess in traditional deformable registration method 4d-ct image difference respiratory phase sequence on spatial domain Correlation, and further contemplate its correlation in time domain, choose suitable wavelet basis function, by setting up rational shape simultaneously Become energy constraint function to meet time domain and the spatial domain constraint of 4d-ct image, effectively improve the performance of registration Algorithm.

4d-ct and the 3d-ct image for the treatment of day by using planning day to rebuild the 4d-ct image for the treatment of day, thus Obtain patient and treat day tumour and the motion change rule jeopardizing organ, effectively prevent patient and accept a large amount of x-ray spokes again Penetrate the injury that body is caused.

Brief description

Fig. 1 is the inventive method flow chart；

Fig. 2 is to be weighed by planning 3d-ct image under phase images and radiotherapy day free breathing state when 10 of day 4d-ct Build the schematic diagram of 4d-ct (the 10 phases) image for the treatment of day；

Fig. 3 (a)-(b) is to decompose based on the edge-protected multiscale space of tv-l1 to carry out decomposability with Gaussian scale-space Comparison diagram；Wherein,

Fig. 3 (a) the first row be tv-l1 multiscale space (first is original image, then respectively using parameter lambda= 0.26,0.17,0.1,0.08,0.06,0.04) multi-resolution decomposition is carried out to the 3d-ct image under free breathing state；

Fig. 3 (a) second row is corresponding residual image；

Fig. 3 (b) the first behavior Gaussian scale-space (first is original image, then respectively using parameter σ=2,4,8, 16,32,64,80) decompose；

Fig. 3 (b) second row is corresponding residual image；

Fig. 4 is the ct Image Multiscale being decomposed based on the edge-protected multiscale space of tv-l1 with Gaussian scale-space contrast Characteristics of decomposition comparison diagram；Wherein,

First behavior tv-l1 multi-resolution decomposition (from left to right be respectively original image, λ=0.7,0.45,0.3,0.2, Tv-l1 multi-resolution decomposition result when 0.15,0.12)；

Second behavior Gauss multi-resolution decomposition (from left to right respectively using parameter σ=2,4,8,16,32,64,80).

Specific embodiment:

The present invention will be further described with embodiment below in conjunction with the accompanying drawings:

By building edge-protected multiscale space, carry out many chis to treating the 3d-ct image under day free breathing state Degree decomposes, and the phase attribute after similarity measure function determines Image Multiscale decomposition.Meanwhile, in conjunction with based on b batten Free deformation grid model to design full-automatic registration approach, act on have different phase attributes radiotherapy day multiple dimensioned 3d-ct and the 4d-ct image of plan day, thus rebuild the treatment day 4d- that can describe radiotherapy same day tumour and jeopardize organ movement Ct image.Concrete steps are as shown in Figure 1:

The first step: therapeutic bed artifact is removed to plan day 4d-ct image, through initial multi-resolution decomposition, removes image In noise or other concussion composition.Then further multi-resolution decomposition is carried out to image, by using different smoothing parameters λ, you can obtain different degrees of smoothed image.

Second step: joint spatial-temporal domain information deformable registration is carried out to a series of smoothed image of difference respiratory phase, sets up Respiratory movement model, is defined as target image.

3rd step: the 3d-ct under free breathing state treatment day being obtained by tv-l1 edge-protected multiscale space Image carries out preliminary exposition, removes picture noise or other shake composition, in edge-protected multi-resolution decomposition further, and lead to Cross similarity measure and determine optimum phase, carry out similarity measure assessment, phase with phase images during 10 planning in day 4d-ct 10 images should be selected, be labeled as incomplete treatment day 4d-ct (10 phases) and be defined as reference picture.

4th step: input target image and reference picture.

5th step: using being intended to phase images (target figure when 10 of day 4d-ct based on the free deformation model of b batten Picture) it is registrated to the corresponding phase treated in day incomplete 4d-ct (reference picture) respectively.

6th step: export the treatment day 4d-ct image that the image after registration terminates is rebuild.

Fig. 2 show and is rebuild by planning 3d-ct under phase images and radiotherapy day free breathing state during day 4d-ct10 The schematic diagram of the 4d-ct (10 phases) for the treatment of day, wherein,

A () is 10 different respiratory phase images in plan day 4d-ct；

B () is the treatment day imperfect 4d-ct image (10 phases) being obtained by multi-resolution decomposition；

C () is the 3d-ct image under the free breathing state that treatment day obtains；

D () is the treatment day 4d-ct image (10 phases) rebuild；

1. joint spatial-temporal domain information carries out deformable registration, sets up breathing fortune to the image of 10 phases in plan day 4d-ct Movable model；

When 2. multi-resolution decomposition being carried out by edge-protected multiscale space, and determining optimum by similarity measure Phase, is labeled as incomplete treatment day 4d-ct, comprises 10 respiratory phase；

3. use phase images when being intended to 10 of day 4d-ct based on the free deformation model of b batten to be registrated to respectively to control Treat the corresponding phase in day incomplete 4d-ct；

4. phase images during difference in treatment day 4d-ct rebuild by algorithm.

For realizing from treating the information obtaining different respiratory phase in the 3d-ct image day free breathing state, this is specially Profit to set up edge-protected multiscale space using Nonlinear Diffusion model, carries out multi-resolution decomposition to treatment day 3d-ct image. In view of Nonlinear Diffusion model tv-l1, not only there is good edge-protected property, there are the many of suitable medical image simultaneously Scale Decomposition characteristic.Assume original image i₀By the larger contour images i of yardstick and less detail pictures v of yardstick (v (x)= i₀- i (x)) form, contour images, than detail pictures (comprising noise and concussion composition) more rule, solve the non-of contour images Linear diffusion model tv-l1 is represented by following energy function

$e(i,\mathrm{\λ})=\underset{\text{i}}{\min }\underset{\mathrm{\ω}}{\&integral;}|\▿i\left(x\right)|\mathrm{dx}+\mathrm{\λ}|i_0\left(x\right)-i\left(x\right)|\mathrm{dx}$

In conjunction with level set function correlation theory, above-mentioned formula is carried out analytical operation, can obtain and there is typical geometric meaning Solution.When carrying out multi-resolution decomposition using tv-l1, original image i is obtained by minimization energy functional e (i, λ)₀In yardstick λ On picture breakdown, i.e. i₀=i (λ)+v (λ), i (λ) image be mainly extract the edge obtaining having regular texture in image and Profile information, and v (λ) is then the noise being removed and little structural information.In this decomposable process, parameter lambda is to determine to put down Slippage degree and the control variables of degree of decomposition.

For example, if λ 1 > λ 2, in yardstick λ₁The image information obtaining in the case of decomposition is included in yardstick λ₂Decomposition situation Under the image information that obtains.Therefore, we can from the beginning of the yardstick of most original (λ=λ₁), obtain to original image i₀Many chis Degree decomposable process is that the profile and edge with different geometry sizes can be carried out selective decomposition by tv-l1, and then can build Multiscale space under Liru decomposes:

i₀=i (λ₁)+v(λ₁)；[i(λ₁),v(λ₁)]=e (i₀,λ₁)；λ₁＞ 0

i(λ₁)=i (λ₂)+v(λ₂)；[i(λ₂),v(λ₂)]=e (i (λ₁),λ₂)；λ₁＞ λ₂

i(λ₂)=i (λ₃)+v(λ₃)；[i(λ₃),v(λ₃)]=e (i (λ₂),λ₃)；λ₂＞ λ₃

i(λ_n-1)=i (λ_n)+v(λ_n)；[i(λ_n),v(λ_n)]=e (i (λ_n-1),λ_n)；λ_n-1＞ λ_n

Here multiscale space decomposition refers to carry out multiple dimensioned point according to the physical dimension size of contour structure in image Solution, and the gray scale with image and further feature are unrelated.On specific yardstick, by changing corresponding parameter lambda in energy function, The smooth removal corresponding less geometry image details of λ, other has profile and the marginal informations of relatively large geometry To be retained.Therefore, in image, geometry profile of different sizes will remain on different scale.Medical image Often it is made up of the tissue and organ of different geometry sizes, this property is suitable for multiple dimensioned point based on tv-l1 just Solution.

Tv-l1 multi-resolution decomposition is different from traditional multiscale transform decomposition, and tv-l1 is based on Nonlinear Diffusion, has preferable side Edge protective nature.Tv-l1 carries out multi-resolution decomposition according to the geometry size of profile in image, on the one hand can be selective Retain the edge in image and profile information, form the image of different scale；On the other hand on the image of particular dimensions, less The image detail of size is smoothed, and the structure of large-size is retained, and the edge of these structures is protected.And passing In the Gaussian scale-space of system, due to isotropic smooth features, in image, all of edge and profile information are all smoothed With fuzzy, on the image of some large scale, the image edge information that Gaussian scale-space obtains offset by the original of them Position.Concrete comparison diagram is as shown in Fig. 3 (a)-(b) and Fig. 4.

One group based on the decomposition of tv-l1 edge-protected multiscale space, by coarse and fine image, can be passed through with this just Control grid degree of roughness to be combined with each other come the deformation model to arrange differing complexity, to realize multiple dimensioned deformable registration frame Frame.On different scale space images, to control the performance recovering deformation using the free deformation grid of different fine dimensions The time being used with registration process.The deformation characteristics of the corresponding overall situation of coarse free deformation grid, and fine free deformation The deformation details of the corresponding local of grid, this property corresponds to the multi-scale image decomposing property of tv-l1 just.Many through tv-l1 After Scale Decomposition, coarse yardstick has few details, in this case, can be come using less free deformation grid Recover deformation, this recovers overall deformation for the higher robustness of holding and has advantage.

For fine dimension, due to comprising more details deformation characteristics it is therefore desirable to the free deformation grid of comparatively dense To recover details deformation, so could keep higher precision.When recovering deformation on a fine scale, it is to be obtained with coarse scale Deformation domain to recover deformation as initial value, therefore on a fine scale registration when still there is speed faster.Therefore adopt With based on tv-l1 edge-protected multiscale space resolution plan day 3d-ct image, in conjunction with based on b batten free deformation grid Lai Recover deformation and can preferably rebuild treatment day 4d-ct image.

For the 3d-ct image treated under day free breathing state, it obtains time more than a respiratory cycle.Therefore, It is the mixed image comprising whole respiratory cycle difference phase.Image-forming principle based on ct we be a reasonable assumption, 3d- In ct image, local detail region is less, and the possibility that it belongs to same respiratory phase is larger.If protected according to based on edge Shield multiscale space can carry out selective decomposition to different minutias in image, then can be the 3d-ct for the treatment of day Picture breakdown is a series of image representing different phases in the same respiratory cycle.

After the edge-protected metric space of tv-l1 obtains a series of multi-resolution decomposition images belonging to different phases, need Determine the phase attribute of these images.Clinically often 4d-ct image is divided into 10 phases according to individual breathing characteristic Image.Therefore, the 3d-ct picture breakdown for the treatment of day is the imperfect 4d-ct image representing 10 different phases by we.Pass through When calculating different from plan day 4d-ct, the similarity measure between phase images to determine the phase attribute of these imperfect images, Due to being a complete image and an imperfect image carries out similarity measurement herein, adopt mutual information here as similar Property measure function.Even this is because the similarity measure function incomplete situation of piece image wherein based on mutual information Under also enable preferable registration.Then, using incomplete difference phase ct image as reference picture, correspond phase Plan day ct image as floating image, carry out deformable registration operation using the free deformation model based on b batten, to rebuild Treatment day 4d-ct different when phase images.

Deformable registration algorithm is using the free deformation grid model based on b batten.In the deformable registration stage, before there is side The advantage of the multi-resolution decomposition property of edge protective nature is, on each specific yardstick, protected edge and profile can To effectively improve precision and the robustness of registration.Furthermore, it is contemplated that the automation of Clinical practice registration Algorithm requires, can pass through The automatic estimation model of smoothing parameter λ in Nonlinear Diffusion model is set up in existing clinical case storehouse.We will be to based on edge The multiscale space decomposition algorithm of protective nature, based on the free deformation model of b batten, the similarity measure letter based on mutual information Number is furtherd investigate, particularly on the methods and strategies that research registration accuracy, speed and robust performance can meet clinical practice. Meanwhile, all contain the information of each respiratory phase in treatment day 3d-ct or 3d-cbct, can in conjunction with plan day 4d-ct image Rebuild treatment day 4d-ct image, can achieve to tumour and the accurate analysis jeopardizing organ movement's rule.

Although the above-mentioned accompanying drawing that combines is described to the specific embodiment of the present invention, not model is protected to the present invention The restriction enclosed, one of ordinary skill in the art should be understood that on the basis of technical scheme, and those skilled in the art are not Need to pay the various modifications that creative work can make or deformation still within protection scope of the present invention.

Claims (6)

1. a kind of method for reconstructing of 4d-ct difference phase sequence image, is characterized in that, comprise the following steps:

The first step: remove plan day 4d-ct image therapeutic bed artifact, and remove the noise in image or other concussion compositions；

Second step: the multi-resolution decomposition based on the edge-protected multiscale space of tv-l1 is carried out to image, puts down by using different Sliding parameter λ, obtain different respiratory phase comprises the how much different smoothed image of information content；

3rd step: joint spatial-temporal domain information deformable registration is carried out to the smoothed image of above-mentioned difference respiratory phase, by deformable registration Smoothed image afterwards is defined as target image；Deformable registration method particularly includes: by described smoothed image with an intensity function f (i, k) representing, wherein, i and k is the spatially and temporally sample of n respiratory phase respectively；During using cycle based on b batten Domain model, to describe the correlation in sequence image time domain, the movement locus of tumour is modeled as smoothing based on the time domain of b batten Function；Consider the snap back conversion in air-breathing end to beginning moment organ movement direction of exhaling in respiratory movement simultaneously, set up segmentation Constraint function is it is ensured that while whole respiratory movement smooth trajectory, keep breathing track accurate in air-breathing end movement locus Property；To express the deformation domain of the respiratory movement correlated series image of 4d-ct difference phase using dual-tree complex wavelet transform, by deformation Domain is modeled as the function of wavelet coefficient；Set up deformation energy function by navier partial differential equation, estimated by coarse and fine The optimal wavelet coefficient of representation space domain deformation characteristics, described wavelet coefficient is substituted into target image transformation model to target image Carry out deformation conversion；

4th step: remove the noise of 3d-ct image under the free breathing state obtaining treatment day or other shake composition；

5th step: multi-resolution decomposition is carried out to 3d-ct image by the edge-protected multiscale space of tv-l1, and by similar Property estimate the optimum phase of determination, be labeled as incomplete treatment day 4d-ct, and incomplete treatment be defined as joining day 4d-ct Examine image；

6th step: respectively the target image of different phases is registrated to reference picture using the free deformation model based on b batten In corresponding phase in；Image after registration terminates is the treatment day 4d-ct image rebuild.

2. a kind of method for reconstructing of 4d-ct difference phase sequence image as claimed in claim 1, is characterized in that, described will swell The movement locus of knurl be modeled as time domain smooth function based on b batten particularly as follows:

$t_t(x,t)=x+\underset{l\&element;l}{\σ}{b}_l{\mathrm{\ψ}}_l\left(t\right)$

Wherein, t_t(x, t) represents t_yMoment when phase images in reference time t_yAnd the movement locus at the x of position, it is with regard to the time The continuously smooth function of t, t_yRepresent the y moment of t, ψ_l(t)=β^m(t/s-l), s ∈ r is the control point interval of time；b_l∈r³For The coefficient of m rank b spline base function；β^mFor m rank b spline base function, β is b spline base function node, and l counts for b spline base function Variable；

The overall smooth track model of b batten cycle Model in Time Domain is:

$\frac{{\∂}^z{t}_t(x,0)}{\∂t^z}=\frac{{\∂}^z{t}_t(x,t_e)}{\∂t^z},z=\[0,...,m-1\]$

Wherein, t_tMovement locus at x for phase images when (x, 0) represented for 0 moment；t_t(x,t_e) represent t the e moment when phasor As the movement locus at x.

3. a kind of method for reconstructing of 4d-ct difference phase sequence image as claimed in claim 1, is characterized in that, described foundation Segmentation constraint function particularly as follows:

t_t ^*(x, 0)=t_t ^*(x,t_e)

Wherein, t_t ^*Sectionally smooth track at x for phase images when (x, 0) represented for 0 moment, t_t ^*(x,t_e) represent t the e moment when Sectionally smooth track at x for the phase images；

t_t ^*For sectionally smooth track it is assumed that air-breathing foot couple answers t=0, single constraints can be set and apply at air-breathing end.

4. a kind of method for reconstructing of 4d-ct difference phase sequence image as claimed in claim 1, is characterized in that, described by shape

Variable domain be modeled as wavelet coefficient function particularly as follows:

X'=x+u¹(x,y,z；c)

Y'=y+u²(x,y,z；c)

Z'=z+u³(x,y,z；c)

Wherein, (x, y, z) represents the space coordinates of reference picture, and (x ', y ', z ') represents the space coordinates of target image, Deformation domain u=(u¹,u²,u³), c is wavelet coefficient, and deformation domain u is the function of wavelet coefficient c.

5. a kind of method for reconstructing of 4d-ct difference phase sequence image as claimed in claim 1, is characterized in that, described passes through Navier partial differential equation setting up deformation energy function particularly as follows:

E (c)=inter (c)+exter (c)

Wherein, inter (c) represents internal force constraint function, and exter (c) represents external force constraint function.

6. a kind of method for reconstructing of 4d-ct difference phase sequence image as claimed in claim 1, is characterized in that, the described 5th In step, determine that by similarity measure the method for the optimum phase of 3d-ct image is: by 3d-ct image and plan day 4d-ct N when phase images carry out similarity measure assessment, accordingly select and plan phase images during immediate n of day 4d-ct.