CN104867104A - Method for obtaining anatomical structural atlas for target mouse based on XCT image non-rigid registration - Google Patents

Abstract

The invention discloses a method for obtaining anatomical structural atlas for a target mouse based on XCT image non-rigid registration. The method includes basic steps of setting Digimouse model as referential mouse anatomical structural atlas, setting XCT images corresponding to Digimouse as referential images, obtaining and pretreating XCT images of the target mouse, building referential images to registration mapping matrix of the target images by non-rigid registration technology, acting registration mapping matrix onto the referential mouse anatomical structural atlas to build a target mouse anatomical structural atlas and mark tissue and organs of the target mouse. The registration method is high in precision, and tissue and organs of a mouse can be marked on a mouse XCT image in a simple and easy way through image registration method. The method is also applicable to other medical science application fields, such as brain structure research, i.e. the anatomical structure of a target brain can be obtained by standard brain anatomical images.

Description

Based on the target mouse anatomical structure collection of illustrative plates acquisition methods of the non-rigidity registration of XCT image

Technical field

The invention belongs to technical field of image processing, the finite element being specifically related to non-rigidity image registration and target rat tissue anatomical structure is demarcated.

Background technology

Current fluorescent molecular tomography (Fluorescence Molecular Tomography, FMT) the normally used uniform optical structural context of method will introduce the appreciable error in photon transport modeling, and the lifting of effective optical texture prior imformation to FMT reconstruction precision and sensitivity is significant.The foundation of optical texture and the acquisition of anatomical information closely related: on the one hand it is the physical geometry information imparting each region related optical parameter characteristic in anatomical structure, and it is the condition precedent that each area optical parameter obtains at body on the other hand ^[1,2].

Common anatomy imaging mode is used for optical texture acquisition and all has some limitations.High-Field toy Magnetic resonance imaging (Micro Magnetic Resonance Imaging, μM RI) has high gray resolution image to soft tissue, and image method can be utilized to obtain each soft tissue organs region.The people such as Dhenain utilize this technology to carry out imaging to multiple mice embryonic, obtain the anatomical structure collection of illustrative plates (Atlas) at different development stage mice embryonic ^[3].The people such as Segars utilize many group adult mice nuclear magnetic resonance datas to set up four-dimensional MOBY model, the dynamic dissection structural drawing spectrum of simulation mouse in the physiology courses such as heartbeat breathing ^[4].But μM RI cost is high, and imaging is consuming time longer, limit its mouse FMT test in application.X ray computer fault imaging (X-ray Computed Tomography, XCT) is as a kind of conventional anatomy imaging pattern, and its image taking speed is very fast and cost is moderate.But the resolution of X ray to soft tissue is lower, utilize XCT Iamge Segmentation soft tissue to have larger difficulty, if be aided with other image mode, then can obtain biological tissue's body anatomical structure more accurately.The people such as Dogdas utilize XCT, positron emission computerized tomography and cold cut chip technology to carry out multi-modality imaging to mouse, and set up Digimouse model on this basis, obtain the precise anatomical structure collection of illustrative plates of mouse ^[5].The method can obtain the anatomical information of mouse exactly, and for correlative study provides important reference significance, but the method imaging system used is complicated, and experimentation is loaded down with trivial details, and cost is higher, is unfavorable for equally adopting in the mouse FMT of routine tests.

[list of references]

[1]L.-H.Wu,W.-B.Wan and X.Wang et al,"Shape-parameterized diffuse opticaltomography holds promise for sensitivity enhancement of fluorescence molecular tomography,"Biomedical Optics Express 10,3640-3659(2014)。

[2]L.-H.Wu,H.-J.Zhao and X.Wang et al,"Enhancement of fluorescence moleculartomography with structural-prior-based diffuse optical tomography:combating opticalbackground uncertainty,"APPLIED OPTICS 53(30),6970-6982(2014)。

[3]M.Dhenain,S.W.Ruffins,and R.E.Jacobs,"Three-Dimensional Digital Mouse AtlasUsing High-Resolution MRI,"Developmental Biology 232,458-470(2001)。

[4]W.P.Segars,B.M.W.Tsui,and E.C.Frey et al,"Development of a 4-D Digital MousePhantom for Molecular Imaging Research,"Molecular Imaging and Biology 6(3),149–159(2004)。

[5]B.Dogdas,D.Stout and A.F.Chatziioannou et al,"Digimouse:a 3D whole body mouseatlas from CT and cryosection data,"PHYSICS IN MEDICINE AND BIOLOGY 52(3),577-587(2007)。

[6]H.Chui and A.Rangarajan,"A new point matching algorithm for non-rigid registration,"Computer Vision and Image Understanding 89,114-141(2003)。

[7]S.Lee,G.Wolberg,and S.Y.Shin,"Scattered Data Interpolation with MultilevelB-Splines,"IEEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS3(3),228-244(1997)。

Summary of the invention

For problems of the prior art, the present invention proposes a kind of target mouse anatomical structure collection of illustrative plates acquisition methods based on the non-rigidity registration of XCT image, and the inventive method is by Digimouse model ^[4]be chosen to be the standard anatomical structure of mouse, and by non-rigidity image registration algorithm, by XCT image registration corresponding for Digimouse model to target mouse XCT image, thus realize the demarcation to each histoorgan of target mouse.

In order to solve the problems of the technologies described above, a kind of target mouse anatomical structure collection of illustrative plates acquisition methods based on the non-rigidity registration of XCT image that the present invention proposes is: first, be with reference to mouse anatomical structure collection of illustrative plates by Digimouse model specification, be reference picture by XCT image setting corresponding for Digimouse simultaneously; Secondly, carry out XCT imaging to target mouse obtain target image and carry out pre-service; Then, non-rigidity image registration techniques is utilized to build the registration mapping matrix of reference picture to target image; Finally, registration mapping matrix is acted on reference to mouse anatomical structure collection of illustrative plates, construct target mouse anatomical structure collection of illustrative plates.

Further, the concrete steps that the present invention is based on the target mouse anatomical structure collection of illustrative plates acquisition methods of XCT image non-rigidity registration are as follows:

Step one, setting are with reference to mouse anatomical structure collection of illustrative plates and reference picture:

Be with reference to mouse anatomical structure collection of illustrative plates A by Digimouse configuration settings _r, be reference picture I by the XCT image setting of described Digimouse _r;

The acquisition of step 2, target mouse XCT image and pre-service:

Utilize XCT equipment to carry out whole body imaging to target mouse, obtain the XCT image of target mouse; Affined transformation is carried out to target mouse XCT image, by the head of target mouse XCT image small mouse and back towards being converted into and reference picture I _ridentical; Using the XCT image after conversion as target image I _t;

Step 3, structure preliminary registration mapping matrix M _cwith preliminary registration image I _c:

Image segmentation algorithm is utilized to be partitioned into reference picture I respectively _rin mouse bony areas and cuticle region, and target image I _tin mouse bony areas and cuticle region;

Edge detection algorithm is utilized to extract above-mentioned reference picture I respectively _rsmall mouse bony areas and cuticle region and target image I _tin mouse bony areas and cuticle region amount to the boundary profile in four regions, equiprobability sampling is carried out to the boundary profile in described four regions, thus calculates with reference to mouse skeleton character point set L _rb, with reference to mouse epidermis characteristic point set L _rs, target mouse skeleton character point set L _tband target mouse epidermis characteristic point set L _ts;

Utilize TPS-RPM (Thin-plate Spline Robust Point Matching, thin plate spline robust point registration) algorithm ^[6]program of increasing income, computing reference mouse skeleton character point set L _rbwith target mouse skeleton character point set L _tbbetween skeleton character point homography C _b; When calculating, by target mouse skeleton character point set L _tbbe set as target point set, with reference to mouse skeleton character point set L _rbbe set as point set subject to registration, and the initial temperature coefficient set in TPS-RPM simulated annealing used and stopping criterion for iteration, go out skeleton character point homography C by iterative computation _b; :

Utilize the program of increasing income of TPS-RPM algorithm, computing reference mouse epidermis characteristic point set L _rswith target mouse epidermis characteristic point set L _tsbetween epidermis characteristic point homography C _s; When calculating, by target mouse epidermis characteristic point set L _tsbe set as target point set, with reference to mouse epidermis characteristic point set L _rsbe set as point set subject to registration, and the initial temperature coefficient set in TPS-RPM simulated annealing used and stopping criterion for iteration, go out epidermis characteristic point homography C by iterative computation _s;

The skeleton character point homography C will tried to achieve respectively _band epidermis characteristic point homography C _sact on reference to mouse skeleton character point set L _rband with reference to mouse epidermis characteristic point set L _rs, obtain the skeleton character point set L after preliminary registration _cband the epidermis characteristic point set L after preliminary registration _cs;

Utilize the result of above-mentioned preliminary Characteristic points match of trying to achieve, build preliminary registration local displacement matrix P:

$P=\left[\begin{array}{cc}{L}_{\mathrm{rb}}& L_{\mathrm{rb}}-L_{\mathrm{cb}}\\ L_{\mathrm{rs}}& L_{\mathrm{rs}}-L_{\mathrm{cs}}\end{array}\right]=[x,y,z,\mathrm{\Δx},\mathrm{\Δy},\mathrm{\Δz}]---\left(3\right)$

Above-mentioned matrix P is converted into three groups of four-dimensional data point set P _x={ (x, y, z, △ x) }, P _y={ (x, y, z, △ y) }, P _z={ (x, y, z, △ z) }, is considered as the functional value of point (x, y, z), i.e. △ x=G respectively by Δ x, Δ y and Δ z ₁(x, y, z), △ y=G ₂(x, y, z), △ z=G ₃(x, y, z);

Adopt Multilevel B-splines fitting algorithm respectively to described three groups of four-dimensional data point set P _x={ (x, y, z, △ x) }, P _y={ (x, y, z, △ y) }, P _z={ (x, y, z, △ z) } carries out matching, for calculating described reference picture I _rin each pixel along the displacement on x, y, z three directions;

Utilize described Multilevel B-splines approximating method fitting data point set P _xthe process of={ (x, y, z, △ x) } is described in detail as follows:

Describedly to refer at many levels, utilize and be overlying on reference picture I _ron one group of cube control mesh Φ encrypted gradually ₀, Φ ₁..., Φ _k..., Φ _hsuccessively to the four-dimensional data point set that iteration upgrades carry out B-spline matching, and by required h layer fitting function sum as final Multilevel B-splines fitting function; Wherein, △ ⁰ξ=△ x, △ ^k+1ξ=△ ^kξ-g _k(x, y, z), g _k(x, y, z) is kth layer B-spline fitting result;

Described kth layer B-spline fit procedure is described below:

Suppose kth layer control mesh Φ _kbe of a size of K _x× K _y× K _z, then kth layer B-spline fitting function is shown below:

$g_k(x,y,z)=\underset{l=0}{\overset{3}{\mathrm{\Σ}}}\underset{m=0}{\overset{3}{\mathrm{\Σ}}}\underset{n=0}{\overset{3}{\mathrm{\Σ}}}{B}_l\left(d_x\right)B_m\left(d_y\right)B_n\left(d_z\right){\mathrm{\φ}}_{k,(l+i,m+j,n+k)}---\left(4\right)$

In formula (4), φ _{k, (l+i, m+j, n+k)}for being positioned at control mesh Φ _kmiddle coordinate is the Controlling vertex value of (l+i, m+j, n+k); L, m, n ∈ { 0,1,2,3}; B _l, B _mand B _nbe respectively l, m, n rank B-spline basis function, wherein the expression formula of 0 to 3 rank B-spline function is described below:

$\left\{\begin{array}{c}{B}_0\left(\mathrm{\δ}\right)={(1-\mathrm{\δ})}^3/6\\ B_1\left(\mathrm{\δ}\right)=(3{\mathrm{\δ}}^3-6{\mathrm{\δ}}^2+4)/6\\ B_2\left(\mathrm{\δ}\right)=(-3{\mathrm{\δ}}^3+3{\mathrm{\δ}}^2+3\mathrm{\δ}+1)/6\\ B_3\left(\mathrm{\δ}\right)={\mathrm{\δ}}^3/6\end{array}\right.---\left(5\right)$

In described formula (4), control mesh Φ _kin each Controlling vertex value by following two steps calculate:

A () calculates in each data point to control mesh Φ _kin the influence amount of each Controlling vertex value:

With in a data point p=(x _p, y _p, z _p, △ ^kξ _p) be described as follows for example:

Data point p=(x _p, y _p, z _p, △ ^kξ _p) to control mesh Φ _kin the influence matrix of each Controlling vertex show as one and be of a size of K _x× K _y× K _zmatrix Ψ _p; Easy for calculating, definition and matrix Ψ _ptwo measure-alike matrix Γ _pwith Ω _p; Described matrix Ψ _p, Γ _pwith Ω _pmiddle coordinate is that the element of (l+i, m+j, n+k) calculates respectively by formula (6):

$\left\{\begin{array}{c}{\mathrm{\ψ}}_{p,(l+i,m+j,n+k)}=\frac{{\mathrm{\γ}}_{p,(l+i,m+j,n+k)}{\mathrm{\Δ}}^k{\mathrm{\ξ}}_p}{{\mathrm{\ω}}_{p,(l+i,m+j,n+k)}}\\ {\mathrm{\γ}}_{p,(l+i,m+j,n+k)}=B_l\left(d_{\mathrm{xp}}\right)B_m\left(d_{\mathrm{yp}}\right)B_n\left(d_{\mathrm{zp}}\right)\\ {\mathrm{\ω}}_{p,(l+i,m+j,n+k)}={\mathrm{\Σ}}_{l=0}^3{\mathrm{\Σ}}_{m=0}^3{\mathrm{\Σ}}_{n=0}^3[B_l\left(d_{\mathrm{xp}}\right)B_m\left(d_{\mathrm{yp}}\right)B_n\left(d_{\mathrm{zp}}\right){]}^2\end{array}\right.---\left(6\right)$

In formula (6), l, m, n ∈ { 0,1,2,3};

At matrix Ψ _p, Γ _pwith Ω _pin, except described coordinate is that (l+i, m+j, n+k) amounts to all the other positions beyond 64 elements, ψ _p, γ _pwith ω _pbe 0;

B () asks for control mesh Φ _kin the value of each Controlling vertex

Comprehensively in each data point to control mesh Φ _kin the impact of each Controlling vertex value, ask for grid Φ processed _kin the value of each Controlling vertex; Control mesh Φ _kmiddle coordinate is the Controlling vertex φ of (a, b, c) _{k, (a, b, c)}value is:

In formula (7), γ _{p, (a, b, c)}, ω _{p, (a, b, c)}, ψ _{p, (a, b, c)}be respectively Ψ _p, Γ _pwith Ω _pmiddle coordinate is the element value of (a, b, c);

So far, kth layer B-spline fitting function g _k(x, y, z) establishes; Comprehensive each level fitting function, calculates Multilevel B-splines fitting function the Multilevel B-splines fitting function g (x, y, z) asked for described in utilization, calculates reference picture I _rin each pixel along the displacement of x-axis;

In like manner, the four-dimensional data point set P of Multilevel B-splines matching is utilized _y={ (x, y, z, △ y) } and P _z={ (x, y, z, △ z) }, thus calculate reference picture I _rin each pixel along y, z-axis in displacement, build reference picture I thus _rpreliminary registration mapping matrix M _c;

Utilize described preliminary registration mapping matrix M _c, oppositely solve preliminary registration image I _c; When building preliminary registration image, the assignment of gray scale adopts tri-linear interpolation methods;

Step 4, build meticulous registration mapping matrix M _f:

Utilize image segmentation algorithm, extract preliminary registration image I _cin mouse skin region and bony areas, and utilize edge detection algorithm to extract the profile of described mouse skin region and bony areas respectively; The profile of mouse skin and bone is superposed mutually, and uses rectangular node to sample, obtain one group of preliminary registration characteristics of image point set L' thus _c;

Block matching method is utilized to ask for preliminary registration characteristics of image point set L' _cin the correspondence position of each unique point on target image, with L' _cin any point p _l=(x _p, y _p, z _p) be example, by as follows for described process prescription:

At preliminary registration image I _cin with coordinate (x _p, y _p, z _p) centered by choose and be of a size of N ₁× N ₁× N ₁cube neighborhood T, at target image I _tin with coordinate (x _p, y _p, z _p) centered by choose N ₂× N ₂× N ₂cube neighborhood S, wherein N ₂>N ₁; Using T as template, S is as region of search, and in S region, search and T have the subregion s of maximum similarity ₁, and by s ₁central point p _l' as putting p _lcorrespondence position;

By that analogy, preliminary registration characteristics of image point set L' is found out successively _cmiddle each point is at target image I _ton correspondence position, construct meticulous registration features point set L thus _f;

Utilize preliminary registration characteristics of image point set L' _cand meticulous registration features point set L _fbuild meticulous registration local displacement matrix Q:

Q＝[L' _c,L' _c-L _f]＝[x,y,z,Δx,Δy,Δz] (8)

Above-mentioned matrix Q is converted into three groups of four-dimensional data point set Q _x={ (x, y, z, △ x) }, Q _y={ (x, y, z, △ y) }, Q _z={ (x, y, z, △ z) }; Consistent with Multilevel B-splines fit procedure in step 3, respectively to described three groups of four-dimensional data point set Q _x={ (x, y, z, △ x) }, Q _y={ (x, y, z, △ y) } carries out Multilevel B-splines matching, thus calculates preliminary registration image I respectively _cin each pixel along the displacement of x, y, z three axis, build preliminary registration image I thus _cmeticulous registration mapping matrix M _f;

Step 5, establishing target mouse anatomical structure collection of illustrative plates:

By described with reference to mouse anatomical structure collection of illustrative plates A _rbe projected to pixel coordinate system, obtain the reference mouse anatomical structure collection of illustrative plates under pixel coordinate system in order successively by preliminary registration mapping matrix M _cand meticulous registration mapping matrix M _fact on the reference mouse anatomical structure collection of illustrative plates under described pixel coordinate system make the reference mouse anatomical structure collection of illustrative plates under described pixel coordinate system produce the distortion identical with step 3 preliminary registration and the meticulous registration process of step 4, obtain the registration mouse anatomical structure collection of illustrative plates under pixel coordinate system by the registration mouse anatomical structure collection of illustrative plates under described pixel coordinate system under being projected to physical coordinates system, obtain the registration mouse anatomical structure collection of illustrative plates A under physical coordinates system _f, the registration mouse anatomical structure collection of illustrative plates A under described physical coordinates system _fbe target mouse anatomical structure collection of illustrative plates.

Compared with prior art, the invention has the beneficial effects as follows:

1. the present invention only uses XCT single mode imaging mode to mouse imaging, and experiment is simple, and cost is lower;

2. the soft tissue identification difficulty that method proposed by the invention can effectively avoid XCT image single mode formation method to cause, simple realizes the identified problems of mouse tissue organ on XCT image;

3. the two step registrations that the present invention uses can realize good registration accuracy, thus provide condition for the correct mark of mouse anatomical structure;

4. the present invention also can use other mouse anatomical structure collection of illustrative plates as a reference in specific implementation process, obtains to the anatomical structure collection of illustrative plates realizing being in target mouse under different figure or developmental stage;

5. main thought of the present invention is equally applicable to other field of medical applications as the research of human brain structure, namely utilizes standard human brain anatomical images and the method for the invention, the acquisition of the anatomical structure of realize target human brain.

Accompanying drawing explanation

Fig. 1 is the target mouse anatomical structure collection of illustrative plates acquisition methods block diagram based on the non-rigidity registration of XCT image;

Fig. 2 is the concrete implementing procedure figure of the target mouse anatomical structure collection of illustrative plates acquisition methods based on XCT image non-rigidity registration that the present invention proposes;

Fig. 3 is the data flow diagram of the present invention in specific implementation process.

Embodiment

Be described in further detail technical solution of the present invention below in conjunction with the drawings and specific embodiments, described specific embodiment only explains the present invention, not in order to limit the present invention.

The basic step that Fig. 1 shows the target mouse anatomical structure collection of illustrative plates acquisition methods that the present invention is based on the non-rigidity registration of XCT image is:

First, be with reference to mouse anatomical structure collection of illustrative plates by Digimouse configuration settings, be reference picture by XCT image setting corresponding for Digimouse simultaneously;

Secondly, carry out XCT imaging to target mouse obtain target image and carry out pre-service;

Then, non-rigidity image registration techniques is utilized to build the registration mapping matrix of reference picture to target image;

Finally, registration mapping matrix is acted on reference to mouse anatomical structure collection of illustrative plates, construct target mouse anatomical structure collection of illustrative plates; Registration mapping matrix can be selected to act on reference picture, construct registering images, in order to evaluation image registration accuracy.

Because different mouse exists larger difference when imaging, first described method for registering carries out preliminary registration adjustment mouse general configuration, secondly uses meticulous registration adjustment image detail part.Therefore, the target mouse anatomical structure collection of illustrative plates acquisition methods that the present invention is based on the non-rigidity registration of XCT image can be refined as 5 steps in implementation process.As shown in Figure 2, the data stream in implementation process as shown in Figure 3 for the process flow diagram of described implementation process.5 steps that the present invention is based in the target mouse anatomical structure collection of illustrative plates acquisition methods specific implementation process of the non-rigidity registration of XCT image are described in detail as follows:

Step one, setting are with reference to mouse anatomical structure collection of illustrative plates and reference picture:

By the Digimouse model proposed in people's documents such as Dogdas ^[5]be set as with reference to mouse anatomical structure collection of illustrative plates A _r, be reference picture I by the XCT image setting of described Digimouse _r;

So far, digitized mouse reference model specification is complete;

The acquisition of step 2, target mouse XCT image and pre-service:

Utilize XCT equipment to carry out whole body imaging to target mouse, obtain the XCT image of target mouse; Affined transformation is carried out to target mouse XCT image, by the head of target mouse XCT image small mouse and back towards being converted into and reference picture I _ridentical; Using the XCT image after conversion as target image I _t;

So far, the target image I of target mouse is come from _tobtain complete;

Step 3, structure preliminary registration mapping matrix M _cwith preliminary registration image I _c:

Because target mouse also exists larger difference with reference to the posture of mouse, size and body internal cavity position, therefore first the present invention carries out preliminary images registration, to adjust the difference of two width image small mouse sizes and figure.The main process of described step 3 is reference picture I _rand target image I _tcarry out preliminary registration, construct and control reference picture I _rproduce the preliminary registration mapping matrix M of non-rigidity deformation _c, by M _cact on reference picture I _r, construct preliminary registration image I _c.Specific implementation process is described in detail as follows:

(3-1) reference picture I is extracted respectively _rand target image I _tfeature point set

Image segmentation algorithm is utilized to be partitioned into reference picture I respectively _rin mouse bony areas and cuticle region, and target image I _tin mouse bony areas and cuticle region;

Edge detection algorithm is utilized to extract above-mentioned reference picture I respectively _rsmall mouse bony areas and cuticle region and target image I _tin mouse bony areas and cuticle region amount to the boundary profile in four regions, equiprobability sampling is carried out to the boundary profile in described four regions, thus calculates with reference to mouse skeleton character point set L _rb, with reference to mouse epidermis characteristic point set L _rs, target mouse skeleton character point set L _tband target mouse epidermis characteristic point set L _ts;

(3-2) based on the Characteristic points match of TPS-RPM algorithm

Utilize the TPS-RPM algorithm that Haili Chui proposes ^[6]program of increasing income, computing reference mouse skeleton character point set L _rbwith target mouse skeleton character point set L _tbbetween skeleton character point homography C _b, with reference to mouse epidermis characteristic point set L _rswith target mouse epidermis characteristic point set L _tsbetween epidermis characteristic point homography C _s;

Described with reference to mouse skeleton character point set L _rbwith target mouse skeleton character point set L _tbbetween skeleton character point homography C _band with reference to mouse epidermis characteristic point set L _rswith target mouse epidermis characteristic point set L _tsbetween epidermis characteristic point homography C _sbe fuzzy homography, in fuzzy homography, each element is the floating number between [0,1], in order to describe the power of degree of correspondence between 2;

Wherein, at computing reference mouse skeleton character point set L _rbwith target mouse skeleton character point set L _tbbetween skeleton character point homography C _btime, by target mouse skeleton character point set L _tbbe set as target point set, with reference to mouse skeleton character point set L _rbbe set as point set subject to registration, and the initial temperature coefficient set in TPS-RPM simulated annealing used and stopping criterion for iteration, the extreme point by the following energy equation of iterative computation:

$E(C_b,f)=\underset{i=1}{\overset{N_{\mathrm{rb}}}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_{\mathrm{tb}}}{\mathrm{\Σ}}}{c}_{b,\mathrm{ij}}{\left|\right|l_{\mathrm{rb},i}-f\left(l_{\mathrm{tb},j}\right)\left|\right|}^2+\mathrm{\λ}{\left|\right|\mathrm{Lf}\left|\right|}^2+T\underset{i=1}{\overset{N_{\mathrm{rb}}}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_{\mathrm{tb}}}{\mathrm{\Σ}}}{c}_{b,\mathrm{ij}}\log c_{b,\mathrm{ij}}-\mathrm{\γ}\underset{i=1}{\overset{N_{\mathrm{rb}}}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_{\mathrm{tb}}}{\mathrm{\Σ}}}{c}_{b,\mathrm{ij}}---\left(1\right)$

In formula (1), c _{b, ij}for homography C _bin element; N _rbfor reference mouse skeleton character point set L _rbcomprise the number of unique point, N _tbfor target mouse skeleton character point set L _tbcomprise the number of point patterns; l _{rb, i}for reference mouse skeleton character point set L _rbin the coordinate of i-th unique point, l _{tb, j}for target mouse skeleton character point set L _tbin the coordinate of a jth unique point; F is thin plate spline function; || Lf|| ²for the smoothness constraint to f, wherein L is differentiating operator; T is for controlling the temperature coefficient of fuzzy degree of correspondence in simulated annealing process; λ is priori smoothness weights; γ is Lu Bang Control Sampled-Data weights; Along with temperature coefficient T reduces gradually, homography C _bwith thin plate spline function f alternating iteration; Homography C _bfinal iteration result be required skeleton character point homography C _b;

At computing reference mouse epidermis characteristic point set L _rswith target mouse epidermis characteristic point set L _tsbetween epidermis characteristic point homography C _stime, by target mouse epidermis characteristic point set L _tsbe set as target point set, with reference to mouse epidermis characteristic point set L _rsbe set as point set subject to registration, and the initial temperature coefficient set in TPS-RPM simulated annealing used and stopping criterion for iteration, the extreme point by the following energy equation of iterative computation:

$E(C_s,f)=\underset{i=1}{\overset{N_{\mathrm{rs}}}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_{\mathrm{ts}}}{\mathrm{\Σ}}}{c}_{s,\mathrm{ij}}{\left|\right|l_{\mathrm{rs},i}-f\left(l_{\mathrm{ts},j}\right)\left|\right|}^2+\mathrm{\λ}{\left|\right|\mathrm{Lf}\left|\right|}^2+T\underset{i=1}{\overset{N_{\mathrm{rs}}}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_{\mathrm{ts}}}{\mathrm{\Σ}}}{c}_{s,\mathrm{ij}}\log c_{s,\mathrm{ij}}-\mathrm{\γ}\underset{i=1}{\overset{N_{\mathrm{rs}}}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_{\mathrm{ts}}}{\mathrm{\Σ}}}{c}_{s,\mathrm{ij}}---\left(2\right)$

In formula (2), c _{s, ij}for homography C _sin element; N _rsfor reference mouse epidermis characteristic point set L _rscomprise the number of unique point, N _tsfor target mouse epidermis characteristic point set L _tscomprise the number of point patterns; l _{rs, i}for reference mouse epidermis characteristic point set L _rsin the coordinate of i-th point, l _{ts, j}for target mouse epidermis characteristic point set L _tsin the coordinate of jth point; Remaining variables implication is identical with formula (1); Along with temperature coefficient T reduces gradually, homography C _swith thin plate spline function f alternating iteration; Homography C _sfinal iteration result be required epidermis characteristic point homography C _s;

The skeleton character point homography C will tried to achieve respectively _band epidermis characteristic point homography C _sact on reference to mouse skeleton character point set L _rband with reference to mouse epidermis characteristic point set L _rs, obtain the skeleton character point set L after preliminary registration _cband the epidermis characteristic point set L after preliminary registration _cs;

(3-3) structure of preliminary registration local displacement matrix P

Utilize the result of above-mentioned preliminary Characteristic points match of trying to achieve, build preliminary registration local displacement matrix P:

$P=\left[\begin{array}{cc}{L}_{\mathrm{rb}}& L_{\mathrm{rb}}-L_{\mathrm{cb}}\\ L_{\mathrm{rs}}& L_{\mathrm{rs}}-L_{\mathrm{cs}}\end{array}\right]=[x,y,z,\mathrm{\Δx},\mathrm{\Δy},\mathrm{\Δz}]---\left(3\right)$

(3-4) preliminary registration mapping matrix M _cstructure

Above-mentioned matrix P is converted into three groups of four-dimensional data point set P _x={ (x, y, z, △ x) }, P _y={ (x, y, z, △ y) }, P _z={ (x, y, z, △ z) }, is considered as the functional value of point (x, y, z), i.e. △ x=G respectively by Δ x, Δ y and Δ z ₁(x, y, z), △ y=G ₂(x, y, z), △ z=G ₃(x, y, z);

Adopt Multilevel B-splines fitting algorithm respectively to described three groups of four-dimensional data point set P _x={ (x, y, z, △ x) }, P _y={ (x, y, z, △ y) }, P _z={ (x, y, z, △ z) } carries out matching, for calculating described reference picture I _rin each pixel along the displacement on x, y, z three directions;

In the present invention, Multilevel B-splines fitting algorithm is the three-dimensional data fitting algorithm that Seungyong Lee proposes ^[7]expansion on four-dimentional space.Utilize described Multilevel B-splines approximating method fitting data point set P _xthe process of={ (x, y, z, △ x) } is described in detail as follows:

Describedly to refer at many levels, utilize and be overlying on reference picture I _ron one group of cube control mesh Φ encrypted gradually ₀, Φ ₁..., Φ _k..., Φ _hsuccessively to the four-dimensional data point set that iteration upgrades carry out B-spline matching, and by required h layer fitting function sum as final Multilevel B-splines fitting function; Wherein, △ ⁰ξ=△ x, △ ^k+1ξ=△ ^kξ-g _k(x, y, z), g _k(x, y, z) is kth layer B-spline fitting result;

Described kth layer B-spline fit procedure is described below:

Suppose kth layer control mesh Φ _kbe of a size of K _x× K _y× K _z, then kth layer B-spline fitting function is shown below:

$g_k(x,y,z)=\underset{l=0}{\overset{3}{\mathrm{\Σ}}}\underset{m=0}{\overset{3}{\mathrm{\Σ}}}\underset{n=0}{\overset{3}{\mathrm{\Σ}}}{B}_l\left(d_x\right)B_m\left(d_y\right)B_n\left(d_z\right){\mathrm{\φ}}_{k,(l+i,m+j,n+k)}---\left(4\right)$

In formula (4), φ _{k, (l+i, m+j, n+k)}for being positioned at control mesh Φ _kmiddle coordinate is the Controlling vertex value of (l+i, m+j, n+k); L, m, n ∈ { 0,1,2,3}; B _l, B _mand B _nbe respectively l, m, n rank B-spline basis function, wherein the expression formula of 0 to 3 rank B-spline function is described below:

$\left\{\begin{array}{c}{B}_0\left(\mathrm{\δ}\right)={(1-\mathrm{\δ})}^3/6\\ B_1\left(\mathrm{\δ}\right)=(3{\mathrm{\δ}}^3-6{\mathrm{\δ}}^2+4)/6\\ B_2\left(\mathrm{\δ}\right)=(-3{\mathrm{\δ}}^3+3{\mathrm{\δ}}^2+3\mathrm{\δ}+1)/6\\ B_3\left(\mathrm{\δ}\right)={\mathrm{\δ}}^3/6\end{array}\right.---\left(5\right)$

In described formula (4), control mesh Φ _kin each Controlling vertex value by following two steps calculate:

A () calculates in each data point to control mesh Φ _kin the influence amount of each Controlling vertex value:

With in a data point p=(x _p, y _p, z _p, △ ^kξ _p) be described as follows for example:

Data point p=(x _p, y _p, z _p, △ ^kξ _p) to control mesh Φ _kin the influence matrix of each Controlling vertex show as one and be of a size of K _x× K _y× K _zmatrix Ψ _p; Easy for calculating, definition and matrix Ψ _ptwo measure-alike matrix Γ _pwith Ω _p; Described matrix Ψ _p, Γ _pwith Ω _pmiddle coordinate is that the element of (l+i, m+j, n+k) calculates respectively by formula (6):

$\left\{\begin{array}{c}{\mathrm{\ψ}}_{p,(l+i,m+j,n+k)}=\frac{{\mathrm{\γ}}_{p,(l+i,m+j,n+k)}{\mathrm{\Δ}}^k{\mathrm{\ξ}}_p}{{\mathrm{\ω}}_{p,(l+i,m+j,n+k)}}\\ {\mathrm{\γ}}_{p,(l+i,m+j,n+k)}=B_l\left(d_{\mathrm{xp}}\right)B_m\left(d_{\mathrm{yp}}\right)B_n\left(d_{\mathrm{zp}}\right)\\ {\mathrm{\ω}}_{p,(l+i,m+j,n+k)}={\mathrm{\Σ}}_{l=0}^3{\mathrm{\Σ}}_{m=0}^3{\mathrm{\Σ}}_{n=0}^3[B_l\left(d_{\mathrm{xp}}\right)B_m\left(d_{\mathrm{yp}}\right)B_n\left(d_{\mathrm{zp}}\right){]}^2\end{array}\right.---\left(6\right)$

In formula (6), l, m, n ∈ { 0,1,2,3};

At matrix Ψ _p, Γ _pwith Ω _pin, except described coordinate is that (l+i, m+j, n+k) amounts to all the other positions beyond 64 elements, ψ _p, γ _pwith ω _pbe 0;

B () asks for control mesh Φ _kin the value of each Controlling vertex

Comprehensively in each data point to control mesh Φ _kin the impact of each Controlling vertex value, ask for grid Φ processed _kin the value of each Controlling vertex; Control mesh Φ _kmiddle coordinate is the Controlling vertex φ of (a, b, c) _{k, (a, b, c)}value is:

In formula (7), γ _{p, (a, b, c)}, ω _{p, (a, b, c)}, ψ _{p, (a, b, c)}be respectively Ψ _p, Γ _pwith Ω _pmiddle coordinate is the element value of (a, b, c);

So far, kth layer B-spline fitting function g _k(x, y, z) establishes; Comprehensive each level fitting function, calculates Multilevel B-splines fitting function the Multilevel B-splines fitting function g (x, y, z) asked for described in utilization, calculates reference picture I _rin each pixel along the displacement of x-axis;

In like manner, the four-dimensional data point set P of Multilevel B-splines matching is utilized _y={ (x, y, z, △ y) } and P _z={ (x, y, z, △ z) }, thus calculate reference picture I _rin each pixel along y, z-axis in displacement, build reference picture I thus _rpreliminary registration mapping matrix M _c;

(3-5) preliminary registration image I _cstructure

Utilize described preliminary registration mapping matrix M _c, oppositely solve preliminary registration image I _c; When building preliminary registration image, the assignment of gray scale adopts tri-linear interpolation methods;

So far, preliminary registration mapping matrix M _c, with preliminary registration image I _cobtain complete.

Step 4, build meticulous registration mapping matrix M _fand meticulous registering images I _f:

After preliminary registration, preliminary registration image I _cwith target image I _tstill there is bigger difference, in order to improve registration accuracy, needing to carry out more meticulous registration.The main process of described step 4 is preliminary registration image I _cand target image I _tcarry out meticulous registration, construct and control preliminary registration image I _cproduce the meticulous registration mapping matrix M of non-rigidity deformation _f, by M _fact on preliminary registration image I _c, construct meticulous registering images I _cf.Specific implementation process is described in detail as follows:

(4-1) preliminary registration characteristics of image point set is again extracted

Utilize image segmentation algorithm, extract preliminary registration image I _cin mouse skin region and bony areas, and utilize edge detection algorithm to extract the profile of described mouse skin region and bony areas respectively; The profile of mouse skin and bone is superposed mutually, and uses rectangular node to sample, obtain one group of preliminary registration characteristics of image point set L' thus _c;

(4-2) based on the Characteristic points match of block matching method

Block matching method is utilized to ask for preliminary registration characteristics of image point set L' _cin the correspondence position of each unique point on target image, with L' _cin any point p _l=(x _p, y _p, z _p) be example, by as follows for described process prescription:

At preliminary registration image I _cin with coordinate (x _p, y _p, z _p) centered by choose and be of a size of N ₁× N ₁× N ₁cube neighborhood T, at target image I _tin with coordinate (x _p, y _p, z _p) centered by choose N ₂× N ₂× N ₂cube neighborhood S, wherein N ₂>N ₁; Using T as template, S is as region of search, and in S region, search and T have the subregion s of maximum similarity ₁, and by s ₁central point p _l' as putting p _lcorrespondence position;

By that analogy, preliminary registration characteristics of image point set L' is found out successively _cmiddle each point is at target image I _ton correspondence position, construct meticulous registration features point set L thus _f;

(4-3) structure of meticulous registration local displacement matrix Q

Utilize preliminary registration characteristics of image point set L' _cand meticulous registration features point set L _fbuild meticulous registration local displacement matrix Q:

Q＝[L' _c,L' _c-L _f]＝[x,y,z,Δx,Δy,Δz] (8)

(4-4) meticulous registration mapping matrix M _fstructure

Above-mentioned matrix Q is converted into three groups of four-dimensional data point set Q _x={ (x, y, z, △ x) }, Q _y={ (x, y, z, △ y) }, Q _z={ (x, y, z, △ z) }; Consistent with Multilevel B-splines fit procedure in step 3, respectively to described three groups of four-dimensional data point set Q _x={ (x, y, z, △ x) }, Q _y={ (x, y, z, △ y) } carries out Multilevel B-splines matching, thus calculates preliminary registration image I respectively _cin each pixel along the displacement of x, y, z three axis, build preliminary registration image I thus _cmeticulous registration mapping matrix M _f;

So far, meticulous registration mapping matrix M _fobtain complete.

(4-5) meticulous registering images I _fstructure (optional)

After this process, can select meticulous registration mapping matrix M _fact on preliminary registration image I _c, construct meticulous registering images I _f; When building preliminary registration image, the assignment of gray scale adopts tri-linear interpolation methods equally; By judging final registering images (i.e. meticulous registering images I _f) and target image I _tsimilarity to evaluate the degree of functioning of registration;

Step 5, establishing target mouse anatomical structure collection of illustrative plates:

In non-rigidity process of image registration, by preliminary registration mapping matrix M _cand meticulous registration mapping matrix M _fcontrol action, reference picture I _rbe deformed into meticulous registering images I _f, meticulous registering images I _fwith target image I _tthere is higher similarity.Therefore, by preliminary registration mapping matrix M _cand meticulous registration mapping matrix M _fact on the anatomical structure collection of illustrative plates A with reference to mouse successively _r, then the anatomical structure collection of illustrative plates of establishing target mouse can be similar to.But, under building on physical coordinates system with reference to anatomical structure collection of illustrative plates (in units of mm), and above-mentioned two registration mapping matrixes build on pixel coordinate system under (in units of pixel), therefore need in the process A _rcarry out coordinate transform.Specific implementation process is described in detail as follows:

(5-1) by described reference mouse anatomical structure collection of illustrative plates A _rbe projected to pixel coordinate system, obtain the reference mouse anatomical structure collection of illustrative plates under pixel coordinate system

(5-2) in order successively by preliminary registration mapping matrix M _cand meticulous registration mapping matrix M _fact on the reference mouse anatomical structure collection of illustrative plates under described pixel coordinate system make the reference mouse anatomical structure collection of illustrative plates under described pixel coordinate system produce the distortion identical with step 3 preliminary registration and the meticulous registration process of step 4, obtain the registration mouse anatomical structure collection of illustrative plates under pixel coordinate system

(5-3) by the registration mouse anatomical structure collection of illustrative plates under described pixel coordinate system under being projected to physical coordinates system, obtain the registration mouse anatomical structure collection of illustrative plates A under physical coordinates system _f, the registration mouse anatomical structure collection of illustrative plates A under described physical coordinates system _fbe target mouse anatomical structure collection of illustrative plates.

Although invention has been described by reference to the accompanying drawings above; but the present invention is not limited to above-mentioned embodiment; above-mentioned embodiment is only schematic; instead of it is restrictive; those of ordinary skill in the art is under enlightenment of the present invention; when not departing from present inventive concept, can also make a lot of distortion, these all belong within protection of the present invention.

Claims (3)

1. the target mouse anatomical structure collection of illustrative plates acquisition methods based on the non-rigidity registration of XCT image, it is characterized in that, its basic step is: first, is with reference to mouse anatomical structure collection of illustrative plates by Digimouse model specification, is reference picture by XCT image setting corresponding for Digimouse simultaneously; Secondly, carry out XCT imaging to target mouse obtain target image and carry out pre-service; Then, non-rigidity image registration techniques is utilized to build the registration mapping matrix of reference picture to target image; Finally, registration mapping matrix is acted on reference to mouse anatomical structure collection of illustrative plates, construct target mouse anatomical structure collection of illustrative plates.

2., according to claim 1 based on the target mouse anatomical structure collection of illustrative plates acquisition methods of the non-rigidity registration of XCT image, it is characterized in that, concrete steps are as follows:

Step one, setting are with reference to mouse anatomical structure collection of illustrative plates and reference picture:

Be with reference to mouse anatomical structure collection of illustrative plates A by Digimouse configuration settings _r, be reference picture I by the XCT image setting of described Digimouse _r;

The acquisition of step 2, target mouse XCT image and pre-service:

Utilize XCT equipment to carry out whole body imaging to target mouse, obtain the XCT image of target mouse; Affined transformation is carried out to target mouse XCT image, by the head of target mouse XCT image small mouse and back towards being converted into and reference picture I _ridentical; Using the XCT image after conversion as target image I _t;

Step 3, structure preliminary registration mapping matrix M _cwith preliminary registration image I _c:

Image segmentation algorithm is utilized to be partitioned into reference picture I respectively _rin mouse bony areas and cuticle region, and target image I _tin mouse bony areas and cuticle region;

Edge detection algorithm is utilized to extract above-mentioned reference picture I respectively _rsmall mouse bony areas and cuticle region and target image I _tin mouse bony areas and cuticle region amount to the boundary profile in four regions, equiprobability sampling is carried out to the boundary profile in described four regions, thus calculates with reference to mouse skeleton character point set L _rb, with reference to mouse epidermis characteristic point set L _rs, target mouse skeleton character point set L _tband target mouse epidermis characteristic point set L _ts;

Utilize the program of increasing income of TPS-RPM algorithm, computing reference mouse skeleton character point set L _rbwith target mouse skeleton character point set L _tbbetween skeleton character point homography C _b; When calculating, by target mouse skeleton character point set L _tbbe set as target point set, with reference to mouse skeleton character point set L _rbbe set as point set subject to registration, and the initial temperature coefficient set in TPS-RPM simulated annealing used and stopping criterion for iteration, go out skeleton character point homography C by iterative computation _b;

Utilize the program of increasing income of TPS-RPM algorithm, computing reference mouse epidermis characteristic point set L _rswith target mouse epidermis characteristic point set L _tsbetween epidermis characteristic point homography C _s; When calculating, by target mouse epidermis characteristic point set L _tsbe set as target point set, with reference to mouse epidermis characteristic point set L _rsbe set as point set subject to registration, and the initial temperature coefficient set in TPS-RPM simulated annealing used and stopping criterion for iteration, go out epidermis characteristic point homography C by iterative computation _s;

The skeleton character point homography C will tried to achieve respectively _band epidermis characteristic point homography C _sact on reference to mouse skeleton character point set L _rband with reference to mouse epidermis characteristic point set L _rs, obtain the skeleton character point set L after preliminary registration _cband the epidermis characteristic point set L after preliminary registration _cs;

Utilize the result of above-mentioned preliminary Characteristic points match of trying to achieve, build preliminary registration local displacement matrix P:

$P=\left[\begin{array}{cc}{L}_{\mathrm{rb}}& L_{\mathrm{rb}}-L_{\mathrm{cb}}\\ L_{\mathrm{rs}}& L_{\mathrm{rs}}-L_{\mathrm{cs}}\end{array}\right]=[x,y,z,\mathrm{\Δx},\mathrm{\Δy},\mathrm{\Δz}]---\left(3\right)$

Above-mentioned matrix P is converted into three groups of four-dimensional data point set P _x={ (x, y, z, △ x) }, P _y={ (x, y, z, △ y) }, P _z={ (x, y, z, △ z) }, is considered as the functional value of point (x, y, z), i.e. △ x=G respectively by Δ x, Δ y and Δ z ₁(x, y, z), △ y=G ₂(x, y, z), △ z=G ₃(x, y, z);

Adopt Multilevel B-splines fitting algorithm respectively to described three groups of four-dimensional data point set P _x={ (x, y, z, △ x) }, P _y={ (x, y, z, △ y) }, P _z={ (x, y, z, △ z) } carries out matching, for calculating described reference picture I _rin each pixel along the displacement on x, y, z three directions;

Utilize described Multilevel B-splines approximating method fitting data point set P _xthe process of={ (x, y, z, △ x) } is described in detail as follows:

Describedly to refer at many levels, utilize and be overlying on reference picture I _ron one group of cube control mesh Φ encrypted gradually ₀, Φ ₁..., Φ _k..., Φ _hsuccessively to the four-dimensional data point set that iteration upgrades carry out B-spline matching, and by required h layer fitting function sum as final Multilevel B-splines fitting function; Wherein, △ ⁰ξ=△ x, △ ^k+1ξ=△ ^kξ-g _k(x, y, z), g _k(x, y, z) is kth layer B-spline fitting result;

Described kth layer B-spline fit procedure is described below:

Suppose kth layer control mesh Φ _kbe of a size of K _x× K _y× K _z, then kth layer B-spline fitting function is shown below:

$g_k(x,y,z)=\underset{l=0}{\overset{3}{\mathrm{\Σ}}}\underset{m=0}{\overset{3}{\mathrm{\Σ}}}\underset{n=0}{\overset{3}{\mathrm{\Σ}}}{B}_l\left(d_x\right)B_m\left(d_y\right)B_n\left(d_z\right){\mathrm{\φ}}_{k,(l+i,m+j,n+k)}---\left(4\right)$

In formula (4), φ _{k, (l+i, m+j, n+k)}for being positioned at control mesh Φ _kmiddle coordinate is the Controlling vertex value of (l+i, m+j, n+k); L, m, n ∈ { 0,1,2,3}; B _l, B _mand B _nbe respectively l, m, n rank B-spline basis function, wherein the expression formula of 0 to 3 rank B-spline function is described below:

$\left\{\begin{array}{c}{B}_0\left(\mathrm{\δ}\right)={(1-\mathrm{\δ})}^3/6\\ B_1\left(\mathrm{\δ}\right)=(3{\mathrm{\δ}}^3-6{\mathrm{\δ}}^2+4)/6\\ B_2\left(\mathrm{\δ}\right)=(-3{\mathrm{\δ}}^3+3{\mathrm{\δ}}^2+3\mathrm{\δ}+1)/6\\ B_3\left(\mathrm{\δ}\right)={\mathrm{\δ}}^3/6\end{array}\right.---\left(5\right)$

In described formula (4), control mesh Φ _kin each Controlling vertex value by following two steps calculate:

A () calculates in each data point to control mesh Φ _kin the influence amount of each Controlling vertex value:

With in a data point p=(x _p, y _p, z _p, △ ^kξ _p) be described as follows for example:

Data point p=(x _p, y _p, z _p, △ ^kξ _p) to control mesh Φ _kin the influence matrix of each Controlling vertex show as one and be of a size of K _x× K _y× K _zmatrix Ψ _p; Easy for calculating, definition and matrix Ψ _ptwo measure-alike matrix Γ _pwith Ω _p; Described matrix Ψ _p, Γ _pwith Ω _pmiddle coordinate is that the element of (l+i, m+j, n+k) calculates respectively by formula (6):

$\left\{\begin{array}{c}{\mathrm{\ψ}}_{p,(l+i,m+j,n+k)}=\frac{{\mathrm{\γ}}_{p,(l+i,m+j,n+k)}{\mathrm{\Δ}}^k{\mathrm{\ξ}}_p}{{\mathrm{\ω}}_{p,(l+i,m+j,n+k)}}\\ {\mathrm{\γ}}_{p,(l+i,m+j,n+k)}=B_l\left(d_{\mathrm{xp}}\right)B_m\left(d_{\mathrm{yp}}\right)B_n\left(d_{\mathrm{zp}}\right)\\ {\mathrm{\ω}}_{p,(l+i,m+j,n+k)}={\mathrm{\Σ}}_{l=0}^3{\mathrm{\Σ}}_{m=0}^3{\mathrm{\Σ}}_{n=0}^3{\left[B_l\left(d_{\mathrm{xp}}\right)B_m\left(d_{\mathrm{yp}}\right)B_n\left(d_{\mathrm{zp}}\right)\right]}^2\end{array}\right.---\left(6\right)$

In formula (6), l, m, n ∈ { 0,1,2,3};

At matrix Ψ _p, Γ _pwith Ω _pin, except described coordinate is that (l+i, m+j, n+k) amounts to all the other positions beyond 64 elements, ψ _p, γ _pwith ω _pbe 0;

B () asks for control mesh Φ _kin the value of each Controlling vertex

Comprehensively in each data point to control mesh Φ _kin the impact of each Controlling vertex value, ask for grid Φ processed _kin the value of each Controlling vertex; Control mesh Φ _kmiddle coordinate is the Controlling vertex φ of (a, b, c) _{k, (a, b, c)}value is:

In formula (7), γ _{p, (a, b, c)}, ω _{p, (a, b, c)}, ψ _{p, (a, b, c)}be respectively Ψ _p, Γ _pwith Ω _pmiddle coordinate is the element value of (a, b, c);

So far, kth layer B-spline fitting function g _k(x, y, z) establishes; Comprehensive each level fitting function, calculates Multilevel B-splines fitting function the Multilevel B-splines fitting function g (x, y, z) asked for described in utilization, calculates reference picture I _rin each pixel along the displacement of x-axis;

In like manner, the four-dimensional data point set P of Multilevel B-splines matching is utilized _y={ (x, y, z, △ y) } and P _z={ (x, y, z, △ z) }, thus calculate reference picture I _rin each pixel along y, z-axis in displacement, build reference picture I thus _rpreliminary registration mapping matrix M _c;

Utilize described preliminary registration mapping matrix M _c, oppositely solve preliminary registration image I _c; When building preliminary registration image, the assignment of gray scale adopts tri-linear interpolation methods;

Step 4, build meticulous registration mapping matrix M _f:

Utilize image segmentation algorithm, extract preliminary registration image I _cin mouse skin region and bony areas, and utilize edge detection algorithm to extract the profile of described mouse skin region and bony areas respectively; The profile of mouse skin and bone is superposed mutually, and uses rectangular node to sample, obtain one group of preliminary registration characteristics of image point set L' thus _c;

Block matching method is utilized to ask for preliminary registration characteristics of image point set L' _cin the correspondence position of each unique point on target image, with L' _cin any point p _l=(x _p, y _p, z _p) be example, by as follows for described process prescription:

At preliminary registration image I _cin with coordinate (x _p, y _p, z _p) centered by choose and be of a size of N ₁× N ₁× N ₁cube neighborhood T, at target image I _tin with coordinate (x _p, y _p, z _p) centered by choose N ₂× N ₂× N ₂cube neighborhood S, wherein N ₂>N ₁; Using T as template, S is as region of search, and in S region, search and T have the subregion s of maximum similarity ₁, and by s ₁central point p ' _las a p _lcorrespondence position;

By that analogy, preliminary registration characteristics of image point set L' is found out successively _cmiddle each point is at target image I _ton correspondence position, construct meticulous registration features point set L thus _f;

Utilize preliminary registration characteristics of image point set L' _cand meticulous registration features point set L _fbuild meticulous registration local displacement matrix Q:

Q＝[L' _c,L' _c-L _f]＝[x,y,z,Δx,Δy,Δz] (8)

Above-mentioned matrix Q is converted into three groups of four-dimensional data point set Q _x={ (x, y, z, △ x) }, Q _y={ (x, y, z, △ y) }, Q _z={ (x, y, z, △ z) }; Consistent with Multilevel B-splines fit procedure in step 3, respectively to described three groups of four-dimensional data point set Q _x={ (x, y, z, △ x) }, Q _y={ (x, y, z, △ y) } carries out Multilevel B-splines matching, thus calculates preliminary registration image I respectively _cin each pixel along the displacement of x, y, z three axis, build preliminary registration image I thus _cmeticulous registration mapping matrix M _f;

Step 5, establishing target mouse anatomical structure collection of illustrative plates:

By described with reference to mouse anatomical structure collection of illustrative plates A _rbe projected to pixel coordinate system, obtain the reference mouse anatomical structure collection of illustrative plates under pixel coordinate system in order successively by preliminary registration mapping matrix M _cand meticulous registration mapping matrix M _fact on the reference mouse anatomical structure collection of illustrative plates under described pixel coordinate system make the reference mouse anatomical structure collection of illustrative plates under described pixel coordinate system produce the distortion identical with step 3 preliminary registration and the meticulous registration process of step 4, obtain the registration mouse anatomical structure collection of illustrative plates under pixel coordinate system by the registration mouse anatomical structure collection of illustrative plates under described pixel coordinate system under being projected to physical coordinates system, obtain the registration mouse anatomical structure collection of illustrative plates A under physical coordinates system _f, the registration mouse anatomical structure collection of illustrative plates A under described physical coordinates system _fbe target mouse anatomical structure collection of illustrative plates.

3. according to claim 2 based on the target mouse anatomical structure collection of illustrative plates acquisition methods of the non-rigidity registration of XCT image, it is characterized in that, in step 3, utilize the program of increasing income of TPS-RPM algorithm to ask for skeleton character point homography C _bwith epidermis characteristic point homography C _sdetailed process as follows:

Described with reference to mouse skeleton character point set L _rbwith target mouse skeleton character point set L _tbbetween skeleton character point homography C _band with reference to mouse epidermis characteristic point set L _rswith target mouse epidermis characteristic point set L _tsbetween epidermis characteristic point homography C _sbe fuzzy homography, in fuzzy homography, each element is the floating number between [0,1], in order to describe the power of degree of correspondence between 2;

Described TPS-RPM algorithm is at computing reference mouse skeleton character point set L _rbwith target mouse skeleton character point set L _tbbetween skeleton character point homography C _btime, utilize simulated annealing to be calculated as follows the minimum point of energy equation:

$E(C_b,f)=\underset{i=1}{\overset{N_{\mathrm{rb}}}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_{\mathrm{tb}}}{\mathrm{\Σ}}}{c}_{b,\mathrm{ij}}{\left|\right|l_{\mathrm{rb},i}-f\left(l_{\mathrm{tb},j}\right)\left|\right|}^2+\mathrm{\λ}{\left|\right|\mathrm{Lf}\left|\right|}^2+T\underset{i=1}{\overset{N_{\mathrm{rb}}}{\text{\Σ}}}\underset{j=1}{\overset{N_{\mathrm{tb}}}{\mathrm{\Σ}}}{c}_{b,\mathrm{ij}}\log c_{b,\mathrm{ij}}-\mathrm{\γ}\underset{i=1}{\overset{N_{\mathrm{rb}}}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_{\mathrm{tb}}}{\mathrm{\Σ}}}{c}_{b,\mathrm{ij}}---\left(1\right)$

In formula (1), c _{b, ij}for homography C _bin element; N _rbfor reference mouse skeleton character point set L _rbcomprise the number of unique point, N _tbfor target mouse skeleton character point set L _tbcomprise the number of point patterns; l _{rb, i}for reference mouse skeleton character point set L _rbin the coordinate of i-th unique point, l _{tb, j}for target mouse skeleton character point set L _tbin the coordinate of a jth unique point; F is thin plate spline function; || Lf|| ²for the smoothness constraint to f, wherein L is differentiating operator; T is for controlling the temperature coefficient of fuzzy degree of correspondence in simulated annealing process; λ is priori smoothness weights; γ is Lu Bang Control Sampled-Data weights; Along with temperature coefficient T reduces gradually, homography C _bwith thin plate spline function f alternating iteration; Homography C _bfinal iteration result be required skeleton character point homography C _b;

Described TPS-RPM algorithm is at computing reference mouse epidermis characteristic point set L _rswith target mouse epidermis characteristic point set L _tsbetween epidermis characteristic point homography C _stime, utilize simulated annealing to be calculated as follows the minimum point of energy equation:

$E(C_s,f)=\underset{i=1}{\overset{N_{\mathrm{rs}}}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_{\mathrm{ts}}}{\mathrm{\Σ}}}{c}_{s,\mathrm{ij}}{\left|\right|l_{\mathrm{rs},i}-f\left(l_{\mathrm{ts},j}\right)\left|\right|}^2+\mathrm{\λ}{\left|\right|\mathrm{Lf}\left|\right|}^2+T\underset{i=1}{\overset{N_{\mathrm{rs}}}{\text{\Σ}}}\underset{j=1}{\overset{N_{\mathrm{ts}}}{\mathrm{\Σ}}}{c}_{s,\mathrm{ij}}\log c_{s,\mathrm{ij}}-\mathrm{\γ}\underset{i=1}{\overset{N_{\mathrm{rs}}}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_{\mathrm{ts}}}{\mathrm{\Σ}}}{c}_{s,\mathrm{ij}}---\left(2\right)$

In formula (2), c _{s, ij}for homography C _sin element; N _rsfor reference mouse epidermis characteristic point set L _rscomprise the number of unique point, N _tsfor target mouse epidermis characteristic point set L _tscomprise the number of point patterns; l _{rs, i}for reference mouse epidermis characteristic point set L _rsin the coordinate of i-th point, l _{ts, j}for target mouse epidermis characteristic point set L _tsin the coordinate of jth point; Remaining variables implication is identical with formula (1); Along with temperature coefficient T reduces gradually, homography C _swith thin plate spline function f alternating iteration; Homography C _sfinal iteration result be required epidermis characteristic point homography C _s.