CN100446034C - Image elastic registrating method based on limited sampling global optimisation - Google Patents

Abstract

The invention opens an elastic image registration method based on limited sampling and global optimization. The steps of the registration are: First read the deformation image; compute the maximum frequency of defined target function, perform the average sampling calculation on the deformation coefficients with two times of the maximum frequency of target function; get the sequence of deformation images; then read into the reference images; deform the deformation images to finish the image registration with the deformation function that corresponds to minimum target function value. Furthermore, start from the coefficient value of deformation function in the rough registration; search the global minimum point in target function with the gradient descent algorithm.

Description

Image elastic registrating method based on limited sampling global optimisation

Technical field

The present invention relates to the graph image conversion in the plane of delineation, be specifically related to a kind of image elastic registrating method of global optimization, be applicable to fields such as image comparison, data fusion, mutation analysis and Target Recognition.

Background technology

Image registration is meant for piece image seeks a kind of spatial alternation, makes corresponding point on it and another width of cloth image reach coupling on the space.This coupling be meant content identical o'clock on two matching images, identical locus is arranged.

Method for registering images commonly used has a variety of, but their step of registration is identical substantially, mainly may further comprise the steps:

(1) constructs the warping function that image is out of shape.Existing method for registering can be divided into two kinds of Rigid Registration method and elastic registrating methods.The warping function that uses in the Rigid Registration method is a kind of rigid deformation function, and it can only carry out simple relatively conversion such as translation, rotation to figure.If pattern distortion is serious, just need to introduce the elastic deformation function during registration, corresponding method for registering is elastic registrating method just.The elastic deformation function has enough versatilities, can approach nonlinear transformation arbitrarily.The elastic deformation function can be divided into two classes substantially: a class is a non-parametric model.In this model, image is regarded as the resilient film of a slice, is out of shape under the effect of external force and internal force, and finally reaches balance.External force is determined that by the difference of reference picture and deformation pattern internal force is determined by the intensity peace slippage degree of film.The most famous in this model is fluid model, also has diffusion model, light stream model etc. in addition.Another kind of is parameterized model, and model need use some parameters and basis function to represent, the computation process of warping function is the CALCULATION OF PARAMETERS process just.The kind of basis function is a lot, and common basis function comprises that polynomial expression, harmonic function, classification basis function, small echo are with B batten etc.

(2) structure needs the similarity criterion of the reference picture and the deformation pattern of registration.The similarity criterion is also referred to as objective function sometimes, and it is to weigh similarity between reference picture and the deformation pattern.Objective function commonly used has mutual information, mean square deviation etc.

(3) objective function is optimized calculating, reaches its extreme value up to objective function.This process constantly is out of shape with warping function deformation pattern exactly, and when objective function reached extreme value, when promptly the similarity between the deformation pattern after reference picture and the distortion was maximum, registration process had also just been finished.Mostly optimization method commonly used in the above-mentioned registration process is some local optimization methods, as Newton method, half Newton method, conjugate vector method etc.The shortcoming of this method is the global extremum that can not guarantee to find objective function, thereby can not guarantee the correctness of registration results.Also can use global optimization approach for elastic registrating method, corresponding global optimization method has simulated annealing, genetic method, tunneling method, method of exhaustion, method such as a lot of somes methods at random at random.Wherein, simulated annealing and genetic method are oversize computing time, and can not guarantee to obtain at last the overall situation and separate; The shortcoming of tunneling method is exactly that its optimization objects has certain limitation, and versatility is bad; And various methods based on limit at random do not have definite calculation methods yet on the quantity Calculation of random point, obtaining the accurate overall situation separates, have only the quantity that strengthens random point, this just causes computing time long, can not guarantee the correctness of last registration results simultaneously.

Summary of the invention

Technical matters to be solved by this invention is the calculated amount that reduces the global extremum point of calculating target function, shortens the time of the elastic registrating of image, improves the robustness of image registration.

It is as described below that the present invention solves the technical solution of above-mentioned technical matters.

A kind of image elastic registrating method based on limited sampling global optimisation, this method is made up of following steps:

(1) reads in deformation pattern, calculate its frequency spectrum, calculate its highest frequency ω by this frequency spectrum with fast fourier transform _m

(2), calculate the highest frequency of objective function again by this frequency spectrum with the frequency spectrum of fourier transform method calculating target function E (C) with respect to deformation coefficient C.This highest frequency is ω _Max=2 ω _mMax (g ' (x, a)), wherein ω _mBe the highest frequency of deformation pattern in the step (1), (x a) is the derivative of warping function to g '.When batten frequency n=1, (x is 1 C) to g '.

(3) the highest frequency ω of 2 objective functions that obtain set by step _Max2 times of spaces to deformation coefficient C correspondence carry out uniform sampling and calculate, obtain the value c of a series of different deformation coefficient C _iPairing warping function g (x, c _i), respectively deformation pattern is out of shape with resulting these warping functions again, obtain deformation pattern series;

(4) read in reference picture, calculate the mean square deviation of the deformation pattern series that reference picture and step (3) obtain respectively, obtain target function value data set E (c _i), relatively the size of these target function values finds wherein minimum target function value E (c ^*), with this minimum target functional value pairing warping function g (x, c ^*) deformation pattern is out of shape, obtain the deformation pattern of thick registration.

(5) the target function value E (c of the minimum that is found with step (4) ^*) pairing deformation coefficient value c ^*Be starting point, with the minimum value of gradient descent method ferret out function E (C)

Use the pairing warping function of this minimum value again Deformation pattern is out of shape, just obtains the deformation pattern of smart registration.

According to prior art as can be known, in fact just two width of cloth images subject to registration, wherein a width of cloth is as the reference image in image registration, and another width of cloth is out of shape deformation pattern as benchmark with reference picture as deformation pattern, makes the two coupling.In view of the above, the deformation pattern behind the registration can be regarded as by the reference picture distortion and obtains.Therefore, the present invention is the same with prior art, will define the warping function and the objective function of contact reference picture and deformation pattern before registration earlier.The define method of warping function and objective function has a variety of, and the present invention's B batten of sampling defines warping function and objective function, is defined as follows:

(1) warping function with B batten structure is:

$g(x,C)=x+\underset{i\∈I_c}{\mathrm{\Σ}}{C}_i{\mathrm{\β}}_n(x/2^w-i)$

Wherein, C is a deformation coefficient, β _nBe the vector product of n B batten, promptly ${\mathrm{\β}}_n\left(x\right)={\mathrm{\Π}}_{k=1}^m{\mathrm{\β}}_n\left(x_k\right),$ X is the m dimensional vector, 2 ^wBe the nodal pitch of B batten, the w representational level, the flexible yardstick between promptly adjacent level is 2, I _cBe the node space set of B batten, i is that node is at node space I _cIn displacement.X is the coordinate of original image, and g (x) is a corresponding coordinate figure after conversion.

(2) establish reference picture and deformation pattern is respectively f _r(x) and f _t(x), (it can be expressed as f with reference picture with the deformation pattern after the warping function distortion _tThe mean square deviation of (g (x, C))) is as objective function, and then objective function can be defined as:

$E\left(C\right)=\frac{1}{\left|\right|I\left|\right|}\underset{x\&Element;\mathrm{<figure-callout\; id="2"\; label="I\; e\; i"\; filenames="C20061012379600131.png"\; state="\{\{state\}\}">I</figure-callout>}}{\mathrm{\&Sigma;}}{e}_i^{}{}_2^{}=\frac{1}{\left|\right|I\left|\right|}\underset{x\&Element;I}{\mathrm{\&Sigma;}}{(f_t\left(g(x,C)\right)-f_r\left(x\right))}^2$

Wherein, I is the pairing space of image coordinate, and (x C) is warping function to g.

According to sampling thheorem, when sample frequency during more than or equal to the twice of a function highest frequency, this function can come out with the complete recovery of these sampled points.The method of the invention, read in deformation pattern earlier, calculate the highest frequency of the objective function of definition, with 2 times of the highest frequency of objective function deformation coefficient being carried out uniform sampling again calculates, obtain warping function series, respectively deformation pattern is out of shape with the warping function series that obtains again, obtain deformation pattern series, read in reference picture then, with the pairing warping function of minimum target functional value deformation pattern is out of shape, finishes thick registration, and then be starting point with the pairing deformation coefficient value of the target function value of minimum, with the point of the global minimum of gradient descent method ferret out function, finish the registration of image.This shows that the present invention adopts limited sampled point just can obtain the point of the global minimum of objective function according to sampling thheorem, can realize the elastic registrating of image, so it is little to have amount of calculation, the remarkable technique effect that the registration time is short.In addition, the present invention has also overcome stochastic sampling needs the operator rule of thumb to carry out manual intervention, is difficult to guarantee the deficiency of registration same effect.

Description of drawings

Fig. 1 is the process flow diagram of the method for the invention;

Fig. 2 is that the sequence number picked out from the MRI figure of one group of heart perfusion is 5 and 3 adjacent two width of cloth image graph, and respectively named as Fig. 2 a and Fig. 2 b, wherein Fig. 2 a is as the reference image, and Fig. 2 b is as deformation pattern; Fig. 2 c is the image after image shown in Fig. 2 b adopts method for registering distortion of the present invention; Fig. 2 d is the error image of image shown in image shown in Fig. 2 b and Fig. 2 c.

Fig. 3 is the influence curve figure of noise to the registration accuracy rate.

Fig. 4 is for when noise variance is 0.1, method for registering of the present invention and with reference to the registration results comparison diagram of method for registering to same deformation pattern, wherein, Fig. 4 a is a reference picture; Fig. 4 b is the deformation pattern of Fig. 4 a at the noise that adds variance 0.1; Fig. 4 c is for being reference picture with Fig. 4 a, when Fig. 4 b is deformation pattern, and the image that adopts the inventive method registration to obtain; Fig. 4 d is for being reference picture with Fig. 4 a, when Fig. 4 b is deformation pattern, and the image that adopts the reference method registration to obtain.

Embodiment

Below by four non-limiting examples invention is elaborated.

Example 1 (registration of two dimensional image):

Needs for medical research, estimate the kinematic parameter of human heart, the inventor has gathered the MRI image of one group of heart perfusion from General Hospital of Tianjin Medical Univ., totally 20 width of cloth, size is 256 * 256 pixels, and therefrom the image of Xuan Zeing is that sequence number is that 5 and 3 adjacent two width of cloth images are reference picture and deformation pattern (referring to Fig. 2 a and Fig. 2 b).The test platform of registration is a PC, is operating system with WindowXP, and Matlab6.5 is a programming language.Concrete registration process is as described below:

Define warping function and objective function before the registration earlier, concrete steps are as follows:

(1) establishes reference picture and deformation pattern and be defined as f respectively _r(x) and f _t(x), wherein variable x is a bivector, then with the contact reference picture of B batten structure and the warping function of deformation pattern is:

$(x,C)=x+\underset{i\&Element;I_c}{\mathrm{\&Sigma;}}{C}_i{\mathrm{\&beta;}}_n(x/2^w-i)$

Wherein C is a deformation coefficient, β _nBe the vector product of n B batten, promptly ${\mathrm{\&beta;}}_n\left(x\right)={\mathrm{\&Pi;}}_{k=1}^m{\mathrm{\&beta;}}_n\left(x_k\right).$ Here the B batten is taken as batten one time.

(2) establish reference picture and deformation pattern equally and be defined as f respectively _r(x) and f _t(x), (it can be expressed as f with reference picture with the deformation pattern after the warping function distortion _tThe mean square deviation of (g (x, C))) is as objective function, and then objective function is defined as:

$E\left(C\right)=\frac{1}{\left|\right|I\left|\right|}\underset{x\&Element;\mathrm{<figure-callout\; id="2"\; label="I\; e\; i"\; filenames="C20061012379600131.png"\; state="\{\{state\}\}">I</figure-callout>}}{\mathrm{\&Sigma;}}{e}_i^{}{}_2^{}=\frac{1}{\left|\right|I\left|\right|}\underset{x\&Element;I}{\mathrm{\&Sigma;}}{(f_t\left(g(x,C)\right)-f_r\left(x\right))}^2$

Wherein x is a bivector.Make objective function E (C) obtain the just corresponding best distortion of warping function g of extreme value, promptly ask g=argmin _{G ∈ G}E (g (x, C)), wherein G is the pairing space of warping function g;

After performing the preceding preparation of registration, can carry out registration, concrete steps are as follows:

(1) read in deformation pattern shown in Fig. 2 b, this image can be expressed as f _t(x), wherein x is a bivector, calculates f with the fast Flourier deformation gauge _t(x) frequency spectrum is:

$F_t\left(\mathrm{\&omega;}\right)={\&Integral;}_{-\&infin;}^{+\&infin;}{f}_t\left(x\right)\&CenterDot;e^{-\mathrm{j\&omega;C}}\mathrm{dC}$

Calculate f by this frequency spectrum _t(x) highest frequency ω _m

${\mathrm{\&omega;}}_m=\underset{F\left(\mathrm{\&omega;}\right)\&NotEqual;0}{\max }\left(\mathrm{\&omega;}\right)$

For the medical image of reality, ω _mCan be very not big, through counting statistics ω _mBe generally between 3～10.For Fig. 2 b, ω as calculated _m Be 5.

(2) with the frequency spectrum of fourier transform method calculating target function E (C) with respect to deformation coefficient C, the highest frequency that is calculated objective function by this frequency spectrum is:

ω _max＝2ω _m·max(g′(x，C))

ω in this formula _mBe the highest frequency of deformation pattern in the step (1), (x C) is the derivative of warping function to g '.The frequency n of B batten is taken as at 1 o'clock among the present invention, at this moment the maximum frequency ω of objective function E (C) _MaxBe 2 ω _m

The highest frequency ω of objective function _Max=2 ω _mMax (the concrete derivation of g ' (x, C)) is as described below:

Objective function E (C) makes Fourier transform for deformation coefficient C:

$E\left(\mathrm{\&omega;}\right)={\&Integral;}_{-\&infin;}^{+\&infin;}\underset{x\&Element;I}{\mathrm{\&Sigma;}}{(f_r\left(x\right)-f_t\left(g(x,C)\right))}^2\&CenterDot;e^{-\mathrm{j\&omega;C}}\mathrm{dC}$

$={\&Integral;}_{-\&infin;}^{+\&infin;}\underset{x\&Element;I}{\mathrm{\&Sigma;}}(f_r^2\left(x\right)-2f_r\left(x\right)f_t\left(g(x,C)\right)+f_t^2\left(g(x,C)\right))\&CenterDot;e^{-\mathrm{j\&omega;C}}\mathrm{dC}$

$=\underset{x\&Element;I}{\mathrm{\&Sigma;}}{\&Integral;}_{-\&infin;}^{+\&infin;}(f_r^2\left(x\right)-2f_r\left(x\right)f_t\left(g(x,C)\right)+f_t^2\left(g(x,C)\right))\&CenterDot;e^{-\mathrm{j\&omega;C}}\mathrm{dC}$

$=\underset{x\&Element;I}{\mathrm{\&Sigma;}}({\&Integral;}_{-\&infin;}^{+\&infin;}{f}_r^2\left(x\right)e^{-\mathrm{j\&omega;C}}\mathrm{dC}-2{\&Integral;}_{-\&infin;}^{+\&infin;}{f}_r\left(x\right)f_t\left(g(x,C)\right)e^{-\mathrm{j\&omega;C}}\mathrm{dC}+{\&Integral;}_{-\&infin;}^{+\&infin;}{f}_t^2\left(g(x,C)\right)e^{-\mathrm{j\&omega;C}}\mathrm{dC})$

Following formula has three, wherein first

$\underset{x\&Element;I}{\mathrm{\&Sigma;}}{\&Integral;}_{-\&infin;}^{+\&infin;}{f}_r^2\left(x\right)e^{-\mathrm{j\&omega;C}}\mathrm{dC}=2\mathrm{\&pi;}\&CenterDot;\mathrm{\&delta;}\left(\mathrm{\&omega;}\right)\&CenterDot;\underset{x\&Element;I}{\mathrm{\&Sigma;}}{f}_r^2\left(x\right)$

Second is

$\underset{x\&Element;I}{\mathrm{\&Sigma;}}{\&Integral;}_{-\&infin;}^{+\&infin;}{f}_r\left(x\right)f_t\left(g(x,C)\right)e^{-\mathrm{j\&omega;C}}\mathrm{dC}=\underset{x\&Element;I}{\mathrm{\&Sigma;}}{f}_r\left(x\right)\&CenterDot;{\&Integral;}_{-\&infin;}^{+\&infin;}{f}_t\left(g(x,C)\right)e^{-\mathrm{j\&omega;C}}\mathrm{dC}$

The highest frequency of following formula is ∫ just _-∞ ^{+ ∞}f _t(g (x, C)) e ^{-j ω C}The highest frequency of dC.Calculate ∫ now _-∞ ^{+ ∞}f _t(g (x, C)) e ^{-j ω C}The highest frequency of dC.

(x C) in arbitrfary point C=a place Taylor expansion, gets g

$g(x,C)=g(x,a)+g^{\&prime;}(x,a)\&CenterDot;(C-a)+\frac{g^{\&prime;\&prime;}(x,a)}{2!}{(C-a)}^2+\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;$

Use first approximation,

g(x，C)＝g(x，a)+g′(x，a)·(C-a)＝(g(x，C)-g′(x，a)·a)+g′(x，a)·C

In step (1), f _t(x) frequency spectrum is F (ω), and its highest frequency is ω _m, f then _t(absolute value of the frequency spectrum of g ' (x, a) C+b) is

Its highest frequency is ω _mG ' (x, a).Because a is the arbitrfary point, so ∫ _-∞ ^{+ ∞}f _t(g (x, C)) e ^{-j ω C}The highest frequency of dC is ω _mMax (g ' (x, C)).

The 3rd is ∫ _-∞ ^{+ ∞}f _t ²(g (x, C)) e ^{-j ω C}DC.Because do not need to ask the frequency spectrum of following formula, just need ask its highest frequency, and f _t(highest frequency of g (x, C)) is ω _m(g ' (x, C)) is so f for max _t ²(highest frequency of g (x, C)) is 2 ω _mMax (g ' (x, C)).

Above three maximum frequency, we as can be seen the highest frequency of objective function E (C) be 2 ω _mMax (g ' (x, C)).

By g (x, definition C) can know, g (x, C) deformation coefficient in is C, so g ' (x, C)=β _n(x/2 ^w-i) β _n(y/2 ^w-j), β _n(x/2 ^w-maximal value i) is relevant with the size of n.When n=1, β _n(x/2 ^w-i) maximal value is 1, so the maximum frequency ω of objective function E (C) _MaxBe 2 ω _m

(3) highest frequency of objective function E (C) is ω _Max=2 ω _m, ω wherein _mHighest frequency for image 2b.So with 2 ω _MaxSample is exactly with 4 ω _mDeformation coefficient C is carried out uniform sampling, and at this moment sampling interval is $T=\frac{\mathrm{\&pi;}}{2{\mathrm{\&omega;}}_m},$ For image 2, T is 0.314..In the field of definition of C (size of field of definition is with relevant to the estimation of anamorphose size, and generally in 10, this example gets 10), adopt a point every 0.314, obtain a series of concrete value c after the sampling _i(i ∈ M, M is the point set of sampling), corresponding also just have a series of concrete warping function g (x, a c _i), with these warping functions image shown in Fig. 2 b is out of shape, obtain deformation pattern series.

(4) read in reference picture shown in Fig. 2 a, calculate the mean square deviation of the deformation pattern series that it and step (3) obtain respectively, just obtain a series of E (c _i) value, relatively these E (c _i) size, therefrom find out their minimum value E (c ^*), Dui Ying warping function g (x, c at this moment ^*) just be a reasonable distortion to image shown in Fig. 2 b (behind the thick registration).

(target function value E (the c of 5 minimums that found with step (4) ^*) corresponding deformation coefficient value c ^*As starting point, with the global minimum point of gradient descent method ferret out function E (C) Detailed process is described below:

With warping function g (x, the C that obtains in step (4) lining ^*) in deformation coefficient C ^*Be initial value.In the i step of iterative process, we are according to existing deformation coefficient c _iCalculate the step-length of upgrading $\mathrm{\&Delta;}{c}_i=-{\mathrm{\&mu;}}_i{\&dtri;}_{c}E\left(c_i\right),$ If at deformation coefficient is c=c _i+ Δ c _iSituation under, the value of E is littler, this step upgrades successfully so, the parameter value before replacing with new parameter value, i.e. c _I+1=c _iTen Δ c _i, μ _I+1=μ _iIf the value of E has increased on the contrary, just get μ _I+1=μ _i/ 2, compute gradient value once more $\mathrm{\&Delta;}{c}_i=-{\mathrm{\&mu;}}_{i+1}{\&dtri;}_{c}E\left(c_i\right);$ This iterative process continues always, until a certain threshold value of Δ c less than prior setting.At c _iAfter determining, warping function has also just been determined.The deformation pattern of the thick registration that step (4) is obtained with this warping function is out of shape, and just obtains the deformation pattern of smart registration.

In above-mentioned optimizing process, need to calculate the partial derivative vector of E (C)

The single order partial derivative of E (C) is

$\frac{\&PartialD;E}{\&PartialD;c_{j,m}}=\frac{1}{\left|\right|I\left|\right|}\underset{i\&Element;I}{\mathrm{\&Sigma;}}\frac{\&PartialD;e_i}{\&PartialD;f_{w\left(i\right)}}\frac{\&PartialD;f_t^c\left(x\right)}{\&PartialD;x_m}|x=g\left(i\right)\frac{\&PartialD;g_m\left(i\right)}{\&PartialD;c_{j.m}}$

Wherein

$\frac{\&PartialD;e_i}{\&PartialD;f_w\left(i\right)}=2{(f_w\left(i\right)-f_r\left(i\right))}^2$

$\frac{\&PartialD;g_m}{\&PartialD;c_{j,m}}={\mathrm{\&beta;}}_{n_m}(x/h-j)$

$\frac{\&PartialD;f_t^c}{\&PartialD;x_m}\left(x\right)=\underset{k\&Element;I}{\mathrm{\&Sigma;}}{b}_k{\mathrm{\&beta;}}_n^{\&prime;}(x_m-k_m)\underset{l=1,l\&NotEqual;m}{\overset{N}{\mathrm{\&Pi;}}}{\mathrm{\&beta;}}_n(x_l-k_l)$

Obtain the global minimum point with objective function E (C) The time pairing warping function

Deformation pattern is out of shape, obtains the deformation pattern of smart registration, image is as seeing shown in Fig. 2 c shown in Fig. 2 b behind the registration.The error image of image shown in image shown in Fig. 2 b and Fig. 2 c shown in Fig. 2 d, the sports ground of the estimated heart of its expression, the dirty sports ground of picture centre is very clear as can be seen.

Example 2 (3-D view registration):

The used 3-D view of present embodiment is the MRI image of choosing from Institute of Medical Information of biomedical engineering institute of Nanfang Medical Univ image data base, and the image size is 128 * 128 * 128 pixels.The test platform of registration is a PC, is operating system with windowXP, and Matlab6.5 is a programming language.Concrete registration process is as described below:

In the present embodiment, the define method of warping function and objective function is identical with example 1, is specially:

(1) establishes reference picture and deformation pattern and divide other f _rBe defined as f _r(x) and f _t(x), wherein x is a tri-vector.Construct the warping function of getting in touch reference picture and deformation pattern with the B batten, obtain the parametrization elastic model:

$g(x,C)=x+\underset{i\&Element;I_c}{\mathrm{\&Sigma;}}{C}_i{\mathrm{\&beta;}}_n(x/2^w-i)$

Wherein C is a deformation coefficient, β _nBe the vector product of n B batten, promptly ${\mathrm{\&beta;}}_n\left(x\right)={\mathrm{\&Pi;}}_{k=1}^m{\mathrm{\&beta;}}_n\left(x_k\right).$ Here the B batten is taken as batten one time.

(2) establish reference picture and deformation pattern equally and be defined as f respectively _r(x) and f _t(x), (it can be expressed as f with reference picture with the deformation pattern after the warping function distortion _tThe mean square deviation of (g (x, C))) is as objective function, and then objective function is defined as:

$E\left(C\right)=\frac{1}{\left|\right|I\left|\right|}\underset{x\&Element;I}{\mathrm{\&Sigma;}}{e}_i^2=\frac{1}{\left|\right|I\left|\right|}\underset{x\&Element;I}{\mathrm{\&Sigma;}}{\left(f_t\left(g(x,C)\right)f_r\left(x\right)\right)}^2$

After performing the preceding preparation of registration, can carry out registration, concrete steps are as follows:

(1) read in the 3 D deformation image, this image can be expressed as f _t(x), wherein x is a tri-vector, calculates its frequency spectrum with fast fourier transform $F_t\left(\mathrm{\&omega;}\right)={\&Integral;}_{-\&infin;}^{+\&infin;}{f}_t\left(x\right)\&CenterDot;e^{-\mathrm{j\&omega;C}}\mathrm{dC},$ Calculate its highest frequency by this frequency spectrum ${\mathrm{\&omega;}}_m=\underset{F\left(\mathrm{\&omega;}\right)\&NotEqual;0}{\max }\left(\mathrm{\&omega;}\right).$

(2), calculate the highest frequency ω of objective function by this frequency spectrum with the frequency spectrum of fourier transform method calculating target function E (C) with respect to deformation coefficient C _Max=2 ω _mMax (g ' (x, a)).ω in this formula _mBe the highest frequency of deformation pattern in the step (1), (x a) is the derivative of warping function to g '.Concrete derivation is identical with example 1.

By g (x, definition C) can know, g (x, C) deformation coefficient in is C, so g ' in the present embodiment (x, C)=β _n(x/2 ^w-i) β _n(y/2 ^w-j) β _n(z/2 ^w-k).When n=1, β _n(x/2 ^w-i) maximal value is 1, so the maximum frequency ω of objective function E (C) _MaxBe 2 ω _m

(3) highest frequency of objective function E (C) is ω _Max=2 ω _m, ω wherein _mHighest frequency for deformation pattern.So with 2 ω _MaxSample is exactly with 4 ω _mDeformation coefficient C is carried out uniform sampling (with respect to two dimensional image, the difference of 3-D view is that sampling is to carry out) in three-dimensional, obtain a series of concrete value c after the sampling _i(i ∈ M, M is the point set of sampling), corresponding also just have a series of concrete warping function g (x, a c _i), be out of shape with these warping function deformation patterns, obtain deformation pattern series.

(4) read in reference picture, calculate the mean square deviation of the deformation pattern series that it and step (3) obtain respectively, just obtain a series of E (c _i) value, relatively these E (c _i) size of value, therefrom find out their minimum value E (c ^*), finish thick registration.

(5) the target function value E (c of the minimum that is found with step (4) ^*) corresponding deformation coefficient value c ^*As starting point, with the global minimum point of gradient descent method ferret out function E (C)

Finish smart registration.Detailed process is with example 1.

Example 3 (registration accuracy rate contrast experiment):

Owing to,, from this laboratory picture library, chosen 16 width of cloth images at random as the reference image so the inventor adopts artificial deformation's method to produce deformation pattern being out of shape the test effect that better to check the method for the invention under the known condition.Selected reference picture wherein, CT image 8 width of cloth, MR image 8 width of cloth, and respectively they are carried out identical deformation process.The size of original image is 256 * 256 pixels.The content of distortion comprise along X and Y-axis respectively translation 10 each pixel, 10 degree and with the distortion elastic deformation of software photoshop realization turn clockwise.With the reference method of the present invention contrast be that Kybic J and Unser M were in the elastic registrating method based on the B batten in the article " Fast Parametric Elastic Image Registration " that " IEEE Trans.Image Processing " the 12nd volume 11 phase 1427-1442 pages or leaves are delivered in 2003.Comparative result is as shown in table 1.By table 1 as seen, algorithm of the present invention has been realized overall 93.75% accuracy rate, and relatively poor with reference to algorithm registration effect, accuracy rate has only 31.25%.The registration effect of algorithm of the present invention is significantly better than the reference algorithm.

Result's contrast of table 116 group test

	The CT image	The MR image	Summation
				Test sample quantity	8	8	16
Registration correct number (the inventive method)	8	7	15
				Accuracy rate (the inventive method)	100％	87.5％	93.75％
Registration correct number (reference method)	3	2	5
				Accuracy rate (reference method)	37.7％	25％	31.25％

Example 4 (antinoise test):

1, robustness comparison test

This people has selected 5 groups from 16 groups of images of example 3, adopt the inventive method and example 3 described reference methods to carry out the registration experiment respectively.In 5 groups of images of experiment, the noise that adds varying strength is checked the robustness of two kinds of methods to noise.The noise that adds is a Gaussian noise, and noise variance is respectively 0.1 and 0.5.The relation of registration accuracy rate and noise variance as shown in Figure 3.As can be seen from Figure 3, the registration accuracy rate of the inventive method on different noise levels all will be higher than with reference to algorithm, and the inventive method is stronger to the robustness of noise.

2, accuracy comparative experiments

For the actual effect of clearer demonstration registration, the inventor has chosen a CT image (shown in Fig. 4 a) again, and to have added variance be 0.1 noise (seeing Fig. 4 b), carries out the registration experiment with the inventive method and example 3 described reference methods respectively then.The registration results of the inventive method is shown in Fig. 4 c, and the registration results of reference method is shown in Fig. 4 d.Can be clear that from test findings algorithm of the present invention has been realized correct registration, and be wrong to the distortion of image in the lower left corner with reference to algorithm.

Claims (1)

1, a kind of image elastic registrating method based on limited sampling global optimisation, this method is made up of following steps:

(1) reads in deformation pattern, calculate its frequency spectrum, calculate its highest frequency again by this frequency spectrum with fast fourier transform;

(2), calculate the highest frequency of objective function again by this frequency spectrum with the frequency spectrum of fourier transform method calculating target function with respect to the deformation coefficient in the warping function;

(3) elder generation carries out uniform sampling calculating in 2 times of spaces to the deformation coefficient correspondence of the highest frequency of 2 objective functions that obtain set by step, obtains warping function series, respectively deformation pattern is out of shape with the warping function series that obtains again, obtains deformation pattern series;

(4) read in reference picture, calculate the mean square deviation of the deformation pattern series that reference picture and step (3) obtain respectively, obtain the target function value data set, the size that compares these target function values, find wherein minimum target functional value, with the pairing warping function of this minimum target functional value deformation pattern is out of shape, obtains the deformation pattern of thick registration;

(5) the pairing deformation coefficient value of target function value of the minimum that is found with step (4) is a starting point, global minimum with gradient descent method ferret out function, with the pairing warping function of this minimum value deformation pattern is out of shape again, just obtains the deformation pattern of smart registration;

Wherein, described warping function and objective function are defined as follows described:

A. use B batten structural deformation function:

$g(x,C)=x+\underset{i\&Element;I_c}{\mathrm{\&Sigma;}}{C}_i{\mathrm{\&beta;}}_n(x/2^w-i)$

Wherein, C is a deformation coefficient, β _nBe the vector product of n B batten, promptly ${\mathrm{\&beta;}}_n\left(x\right)={\mathrm{\&Pi;}}_{k=1}^m{\mathrm{\&beta;}}_n\left(x_k\right),$ X is the m dimensional vector, and w representational level, Ic are the node space set of B batten, and i is the displacement of node in node space Ic;

B. establish reference picture and deformation pattern and be defined as f respectively _r(x) and f _t(x), with reference picture and with the mean square deviation of the deformation pattern after the warping function distortion as objective function:

$E\left(C\right)=\frac{1}{\left|\right|I\left|\right|}\underset{x\&Element;I}{\mathrm{\&Sigma;}}{e}_i^2=\frac{1}{\left|\right|I\left|\right|}\underset{x\&Element;I}{\mathrm{\&Sigma;}}{(f_t\left(g(x,C)\right)-f_r\left(x\right))}^2$

Wherein, I is the pairing space of image coordinate, and (x C) is warping function to g.

Images