CN104867104B - Target mouse anatomical structure collection of illustrative plates acquisition methods based on the non-rigidity registration of XCT images - Google Patents

Target mouse anatomical structure collection of illustrative plates acquisition methods based on the non-rigidity registration of XCT images Download PDF

Info

Publication number
CN104867104B
CN104867104B CN201510259870.5A CN201510259870A CN104867104B CN 104867104 B CN104867104 B CN 104867104B CN 201510259870 A CN201510259870 A CN 201510259870A CN 104867104 B CN104867104 B CN 104867104B
Authority
CN
China
Prior art keywords
mrow
msub
mouse
point set
target
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201510259870.5A
Other languages
Chinese (zh)
Other versions
CN104867104A (en
Inventor
高峰
万文博
赵会娟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tianjin University
Original Assignee
Tianjin University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Tianjin University filed Critical Tianjin University
Priority to CN201510259870.5A priority Critical patent/CN104867104B/en
Publication of CN104867104A publication Critical patent/CN104867104A/en
Application granted granted Critical
Publication of CN104867104B publication Critical patent/CN104867104B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00Geometric image transformations in the plane of the image
    • G06T3/14Transformations for image registration, e.g. adjusting or mapping for alignment of images
    • G06T3/147Transformations for image registration, e.g. adjusting or mapping for alignment of images using affine transformations
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/12Edge-based segmentation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10081Computed x-ray tomography [CT]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Image Processing (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)

Abstract

The invention discloses a kind of target mouse anatomical structure collection of illustrative plates acquisition methods based on the non-rigidity registration of XCT images, its basic step is:First, it with reference to mouse anatomical structure collection of illustrative plates, while by XCT image settings corresponding to Digimouse is reference picture to be by Digimouse model specifications;Secondly, XCT imagings are carried out to target mouse to obtain target image and pre-process;Then, reference picture is built to the registering mapping matrix of target image using non-rigidity image registration techniques;Finally, registering mapping matrix is acted on reference to mouse anatomical structure collection of illustrative plates, target mouse anatomical structure collection of illustrative plates is constructed, so as to realize the demarcation to each histoorgan of target mouse.Method for registering precision used herein is high, can by the method for image registration it is simple and easy realize mark of the mouse tissue organ on mouse XCT images;It is proposed method of the present invention is equally applicable to other field of medical applications such as human brain structure and studied, you can to realize the acquisition of the anatomical structure of target human brain using standard human brain anatomical images.

Description

Target mouse anatomical structure collection of illustrative plates acquisition methods based on the non-rigidity registration of XCT images
Technical field
The invention belongs to technical field of image processing, and in particular to non-rigidity image registration and target rat tissue anatomical structure Finite element demarcation.
Background technology
Fluorescent molecular tomography (Fluorescence Molecular Tomography, FMT) method generally makes at present The appreciable error that uniform optical structural context will be introduced in photon transport modeling, effective optical texture prior information pair The lifting of FMT reconstruction precisions and sensitivity is significant.The foundation of optical texture and the acquisition of anatomical information are close It is related:On the one hand it has been to confer to the physical geometry information of each region related optical parameter characteristic in anatomical structure, on the other hand It is the prerequisite that each area optical parameter obtains in body[1,2]
Common anatomy imaging mode is used for optical texture acquisition and had some limitations.High-Field toy nuclear-magnetism is total to Imaging (Micro Magnetic Resonance Imaging, μM RI) of shaking has high gray resolution image, Neng Gouli to soft tissue Each soft tissue organs region is obtained with image method.Dhenain et al. is carried out using the technology to multiple mice embryonics Imaging, obtains the anatomical structure collection of illustrative plates (Atlas) in different development stage mice embryonic[3].Segars et al. utilizes form more Year, mouse nuclear magnetic resonance data established four-dimensional MOBY models, and dynamic of the simulation mouse in the physiology courses such as heartbeat breathing dissects knot Composition is composed[4].However, μM RI costs are high, imaging is time-consuming longer, limits its application in mouse FMT experiments.Roentgenometer Calculation machine fault imaging (X-ray Computed Tomography, XCT) as a kind of conventional anatomy imaging pattern, its into As speed and cost is moderate.But X ray is relatively low to the resolution ratio of soft tissue, have using XCT images segmentation soft tissue Larger difficulty, if being aided with other image modes, it more can accurately obtain biological tissue's body anatomical structure.Dogdas et al. profits Multi-modality imaging is carried out to mouse with XCT, positron emission computerized tomography and cold cut chip technology, and built on this basis Vertical Digimouse models, obtain the precise anatomical structure collection of illustrative plates of mouse[5].The method can obtain the dissection of mouse exactly Structural information, important reference significance is provided for correlative study, but imaging system used in this method is complicated, and experimentation is numerous Trivial, cost is higher, is equally unfavorable for using in the mouse FMT experiments of routine.
[bibliography]
[1]L.-H.Wu,W.-B.Wan and X.Wang et al,"Shape-parameterized diffuse optical tomography 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 molecular tomography with structural-prior-based diffuse optical tomography: combating optical background uncertainty,"APPLIED OPTICS 53(30),6970-6982 (2014)。
[3]M.Dhenain,S.W.Ruffins,and R.E.Jacobs,"Three-Dimensional Digital Mouse Atlas Using 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 Mouse Phantom 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 mouse atlas 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 Multilevel B-Splines,"IEEE TRANSACTIONS ON VISUALIZATION AND COMPUTER GRAPHICS3(3),228-244(1997)。
The content of the invention
For problems of the prior art, the present invention proposes a kind of target mouse based on the non-rigidity registration of XCT images Anatomical structure collection of illustrative plates acquisition methods, the inventive method is by Digimouse models[4]The standard anatomical structure of mouse is chosen to be, and is led to Non- rigidity image registration algorithm is crossed, by XCT image registrations corresponding to Digimouse models to target mouse XCT images, so as to realize Demarcation to each histoorgan of target mouse.
In order to solve the above-mentioned technical problem, a kind of target mouse solution based on the non-rigidity registration of XCT images proposed by the present invention Cuing open structure collection of illustrative plates acquisition methods is:First, it is with reference to mouse anatomical structure collection of illustrative plates by Digimouse model specifications, simultaneously will XCT image settings corresponding to Digimouse are reference picture;Secondly, XCT imaging acquisition target images are carried out to target mouse to go forward side by side Row pretreatment;Then, reference picture is built to the registering mapping matrix of target image using non-rigidity image registration techniques;Most Eventually, registering mapping matrix is acted on reference to mouse anatomical structure collection of illustrative plates, constructs target mouse anatomical structure collection of illustrative plates.
Further, the tool of target mouse anatomical structure collection of illustrative plates acquisition methods of the present invention based on the non-rigidity registration of XCT images Body step is as follows:
Step 1: setting refers to mouse anatomical structure collection of illustrative plates and reference picture:
It is with reference to mouse anatomical structure collection of illustrative plates A by Digimouse configuration settingsr, the XCT images of the Digimouse are set It is set to reference picture Ir
Step 2: the acquisition and pretreatment of target mouse XCT images:
Whole body imaging is carried out to target mouse using XCT equipment, obtains the XCT images of target mouse;Target mouse XCT images are entered Row affine transformation, the head of mouse and the direction at back in target mouse XCT images are converted into and reference picture IrIt is identical;It will become XCT images after changing are as target image It
Step 3: structure preliminary registration mapping matrix McWith preliminary registration image Ic
It is partitioned into reference picture I respectively using image segmentation algorithmrIn mouse bony areas and cuticle region, and target Image ItIn mouse bony areas and cuticle region;
Extract above-mentioned reference picture I respectively using edge detection algorithmrMiddle mouse bony areas and cuticle region and target Image ItIn mouse bony areas and cuticle region amount to the boundary profile in four regions, the border in four regions is taken turns Exterior feature carries out equiprobability sampling, and mouse skeleton character point set L is referred to so as to calculaterb, with reference to mouse epidermis characteristic point set Lrs, target mouse Skeleton character point set LtbAnd target mouse epidermis characteristic point set Lts
Using TPS-RPM, (Thin-plate Spline Robust Point Matching, thin plate spline robust point are matched somebody with somebody It is accurate) algorithm[6]Program of increasing income, calculating refer to mouse skeleton character point set LrbWith target mouse skeleton character point set LtbBetween bone Characteristic point homography Cb;When calculating, by target mouse skeleton character point set LtbIt is set as target point set, it is special mouse bone will be referred to Levy point set LrbIt is set as point set subject to registration, and sets initial temperature coefficient and iteration in simulated annealing used in TPS-RPM End condition, by iterating to calculate out skeleton character point homography Cb;:
Using the program of increasing income of TPS-RPM algorithms, calculating refers to mouse epidermis characteristic point set LrsWith target mouse epidermis characteristic point Collect LtsBetween epidermis characteristic point homography Cs;When calculating, by target mouse epidermis characteristic point set LtsIt is set as target point set, Mouse epidermis characteristic point set L will be referred torsIt is set as point set subject to registration, and sets initial in simulated annealing used in TPS-RPM Temperature coefficient and stopping criterion for iteration, by iterating to calculate out epidermis characteristic point homography Cs
Respectively by the skeleton character tried to achieve point homography CbAnd epidermis characteristic point homography CsAct on and refer to mouse bone Feature point set LrbAnd with reference to mouse epidermis characteristic point set Lrs, obtain the skeleton character point set L after preliminary registrationcbAnd after preliminary registration Epidermis characteristic point set Lcs
Utilize the result of the above-mentioned preliminary Characteristic points match tried to achieve, structure preliminary registration local displacement matrix P:
Above-mentioned matrix P is converted into three groups of 4 D data point set Px={ (x, y, z, △ x) }, Py={ (x, y, z, △ y) }, Pz={ (x, y, z, △ z) }, Δ x, Δ y and Δ z are respectively seen as to point (x, y, z) functional value, i.e. △ x=G1(x, y, z), △ y =G2(x, y, z), △ z=G3(x,y,z);
Using Multilevel B-splines fitting algorithm respectively to three groups of 4 D data point sets Px={ (x, y, z, △ x) }, Py ={ (x, y, z, △ y) }, Pz={ (x, y, z, △ z) } is fitted, for calculating the reference picture IrIn each pixel along Displacement on three directions of x, y, z;
Utilize the Multilevel B-splines approximating method fitting data point set Px={ (x, y, z, △ x) } process is retouched in detail State as follows:
It is described to refer at many levels, using being overlying on reference picture IrOn one group of cube control grid Φ gradually encrypted0, Φ1,...,Φk,...,ΦhSuccessively to the 4 D data point set of iteration renewalB-spline fitting is carried out, and By required h layers fitting function sumAs final Multilevel B-splines fitting function;Wherein, △0ξ=△ x, △k+1ξ=△kξ-gk(x, y, z), gk(x, y, z) is kth layer B-spline fitting result;
The kth layer B-spline fit procedure is described as follows:
Assuming that kth layer control grid ΦkSize is Kx×Ky×Kz, then kth layer B-spline fitting function be shown below:
In formula (4), φk,(l+i,m+j,n+k)For positioned at control grid ΦkMiddle coordinate is The control node value of (l+i, m+j, n+k);l,m,n∈{0,1,2,3};Bl、BmAnd BnRespectively l, m, n rank B-spline base letter Number, wherein the expression formula of 0 to 3 rank B-spline function is described as follows:
In the formula (4), grid Φ is controlledkIn the value of each control node calculated by following two step:
(a) calculateIn each data point to controlling grid ΦkIn each control node value shadow Ring amount:
WithIn a data point p=(xp,yp,zp,△kξp) exemplified by be described as follows:
Data point p=(xp,yp,zp,△kξp) to controlling grid ΦkIn the influence matrix of each control node show as one Individual size is Kx×Ky×KzMatrix Ψp;It is easy to calculate, definition and matrix ΨpTwo matrix Γ of size identicalpWith Ωp; The matrix Ψp、ΓpWith ΩpMiddle coordinate is calculated by formula (6) respectively for the element of (l+i, m+j, n+k):
In formula (6), l, m, n ∈ { 0,1,2,3 };
In matrix Ψp、ΓpWith ΩpIn, except the coordinate is that (l+i, m+j, n+k) amounts to remaining position beyond 64 elements Put, ψp、γpWith ωpIt is 0;
(b) ask for controlling grid ΦkIn each control node value
It is comprehensiveIn each data point to controlling grid ΦkIn each control node value influence, Ask for grid Φ processedkIn each control node value;Control grid ΦkMiddle coordinate is the control node φ of (a, b, c)k,(a,b,c)Take It is worth and is:
In formula (7), γp,(a,b,c)、ωp,(a,b,c)、ψp,(a,b,c)Respectively Ψp、ΓpWith ΩpMiddle coordinate is the member of (a, b, c) Element value;
So far, kth layer B-spline fitting function gk(x, y, z) is established;Comprehensive each level fitting function, calculates Multilevel B Spline-fit functionUsing the Multilevel B-splines fitting function g (x, y, z) asked for, Calculate reference picture IrIn displacement of each pixel along x-axis;
Similarly, it is fitted 4 D data point set P using Multilevel B-splinesy={ (x, y, z, △ y) } and Pz={ (x, y, z, △ Z) }, so as to calculating reference picture IrIn displacement of each pixel along y, z-axis in, thus build reference picture IrIt is preliminary Registering mapping matrix Mc
Utilize the preliminary registration mapping matrix Mc, reversely solve preliminary registration image Ic;In structure preliminary registration image When, the assignment of gray scale uses tri-linear interpolation methods;
Step 4: the fine registering mapping matrix M of structuref
Using image segmentation algorithm, preliminary registration image I is extractedcIn mouse skin region and bony areas, and profit Extract the profile of the mouse skin region and bony areas respectively with edge detection algorithm;By mouse skin and the wheel of bone Exterior feature is overlapped mutually, and is sampled using rectangular mesh, thus obtains one group of preliminary registration characteristics of image point set L'c
Preliminary registration characteristics of image point set L' is asked for using block matching methodcIn each pair of the characteristic point on target image Position is answered, with L'cIn any point pl=(xp,yp,zp) exemplified by, the process description is as follows:
In preliminary registration image IcIn with coordinate (xp,yp,zp) centered on choose size be N1×N1×N1Cube neighborhood T, in target image ItIn with coordinate (xp,yp,zp) centered on choose N2×N2×N2Cube neighborhood S, wherein N2>N1;T is made For template, S is as region of search, and search and T have the subregion s of maximum similarity in S regions1, and by s1Central point pl′ As point plCorrespondence position;
By that analogy, preliminary registration characteristics of image point set L' is searched out successivelycMiddle each point is in target image ItOn corresponding position Put, thus construct fine registration features point set Lf
Utilize preliminary registration characteristics of image point set L'cAnd fine registration features point set LfThe fine registering local displacement square of structure Battle array Q:
Q=[L'c,L'c-Lf]=[x, y, z, Δ x, Δ y, Δ z] (8)
Above-mentioned matrix Q is converted into three groups of 4 D data point set Qx={ (x, y, z, △ x) }, Qy={ (x, y, z, △ y) }, Qz={ (x, y, z, △ z) };It is consistent with Multilevel B-splines fit procedure in step 3, respectively to three groups of 4 D datas point Collect Qx={ (x, y, z, △ x) }, Qy={ (x, y, z, △ y) } carries out Multilevel B-splines fitting, tentatively matches somebody with somebody so as to calculate respectively Quasi- image IcIn displacement of each pixel along three axial directions of x, y, z, thus build preliminary registration image IcIt is fine registration mapping Matrix Mf
Step 5: structure target mouse anatomical structure collection of illustrative plates:
By described with reference to mouse anatomical structure collection of illustrative plates ArProjection obtains the reference mouse solution under pixel coordinate system to pixel coordinate system Cut open structure collection of illustrative platesIn order successively by preliminary registration mapping matrix McAnd fine registering mapping matrix MfAct on the pixel Reference mouse anatomical structure collection of illustrative plates under coordinate systemMake the reference mouse anatomical structure collection of illustrative plates under the pixel coordinate systemProduce Deformed with step 3 preliminary registration and the fine registration process identical of step 4, obtain the registering mouse dissection knot under pixel coordinate system Composition is composedBy the registering mouse anatomical structure collection of illustrative plates under the pixel coordinate systemUnder projection to physical coordinates system, thing is obtained Manage the registering mouse anatomical structure collection of illustrative plates A under coordinate systemf, the registering mouse anatomical structure collection of illustrative plates A under the physical coordinates systemfAs mesh Mark mouse anatomical structure collection of illustrative plates.
Compared with prior art, the beneficial effects of the invention are as follows:
1. the present invention is used only XCT single modes imaging mode and mouse is imaged, experiment is simple, and cost is relatively low;
Soft tissue caused by 2. method proposed by the invention can effectively avoid XCT image single mode imaging methods recognizes Difficulty, simple and easy realizes mark problem of the mouse tissue organ on XCT images;
3. the two step registrations that the present invention uses can realize preferable registration accuracy, so as to be the correct of mouse anatomical structure Offer condition is provided;
4. other mouse anatomical structure collection of illustrative plates also can be used as reference in the present invention in specific implementation process, to realize The anatomical structure collection of illustrative plates of target mouse obtains under different figures or developmental stage;
Studied 5. the main thought of the present invention is equally applicable to other field of medical applications such as human brain structure, that is, utilize mark Quasi- human brain anatomical images and the method for the invention, realize the acquisition of the anatomical structure of target human brain.
Brief description of the drawings
Fig. 1 is the target mouse anatomical structure collection of illustrative plates acquisition methods block diagram based on the non-rigidity registration of XCT images;
Fig. 2 is the target mouse anatomical structure collection of illustrative plates acquisition methods proposed by the present invention based on the non-rigidity registration of XCT images Flow chart is embodied;
Fig. 3 is DFD of the present invention in specific implementation process.
Embodiment
Technical solution of the present invention is described in further detail with specific embodiment below in conjunction with the accompanying drawings, described is specific Only the present invention is explained for embodiment, is not intended to limit the invention.
Fig. 1 shows the base of target mouse anatomical structure collection of illustrative plates acquisition methods of the present invention based on the non-rigidity registration of XCT images This step is:
First, it is with reference to mouse anatomical structure collection of illustrative plates, while by XCT corresponding to Digimouse by Digimouse configuration settings Image setting is reference picture;
Secondly, XCT imagings are carried out to target mouse to obtain target image and pre-process;
Then, reference picture is built to the registering mapping matrix of target image using non-rigidity image registration techniques;
Finally, registering mapping matrix is acted on reference to mouse anatomical structure collection of illustrative plates, constructs target mouse anatomical structure collection of illustrative plates; It may be selected registering mapping matrix acting on reference picture, registering image constructed, to evaluation image registration accuracy.
Because different mouse have larger difference in imaging, the method for registering carries out preliminary registration adjustment mouse first General configuration, secondly adjust image detail part using fine registration.Therefore, the present invention is based on the non-rigidity registration of XCT images Target mouse anatomical structure collection of illustrative plates acquisition methods can be refined as 5 steps in implementation process.The flow chart of the implementation process is such as Shown in Fig. 2, the data flow in implementation process is as shown in Figure 3.Target mouse of the present invention based on the non-rigidity registration of XCT images is dissected 5 steps in structure collection of illustrative plates acquisition methods specific implementation process describe in detail as follows:
Step 1: setting refers to mouse anatomical structure collection of illustrative plates and reference picture:
By the Digimouse models proposed in Dogdas et al. documents[5]It is set as referring to mouse anatomical structure collection of illustrative plates Ar, It is reference picture I by the XCT image settings of the Digimouser
So far, digitized mouse reference model specification finishes;
Step 2: the acquisition and pretreatment of target mouse XCT images:
Whole body imaging is carried out to target mouse using XCT equipment, obtains the XCT images of target mouse;Target mouse XCT images are entered Row affine transformation, the head of mouse and the direction at back in target mouse XCT images are converted into and reference picture IrIt is identical;It will become XCT images after changing are as target image It
So far, the target image I of target mouse is come fromtAcquisition finishes;
Step 3: structure preliminary registration mapping matrix McWith preliminary registration image Ic
Due to target mouse with posture, size and the internal empty cavity position with reference to mouse there is larger difference, therefore this hair It is bright to carry out preliminary images registration first, to adjust mouse size and the difference of figure in two images.The step 3 it is main Process is to reference picture IrAnd target image ItPreliminary registration is carried out, constructs control reference picture IrProduce non-rigidity shape The preliminary registration mapping matrix M of changec, by McAct on reference picture Ir, construct preliminary registration image Ic.Specific implementation process It is described in detail as follows:
(3-1) extracts reference picture I respectivelyrAnd target image ItFeature point set
It is partitioned into reference picture I respectively using image segmentation algorithmrIn mouse bony areas and cuticle region, and target Image ItIn mouse bony areas and cuticle region;
Extract above-mentioned reference picture I respectively using edge detection algorithmrMiddle mouse bony areas and cuticle region and target Image ItIn mouse bony areas and cuticle region amount to the boundary profile in four regions, the border in four regions is taken turns Exterior feature carries out equiprobability sampling, and mouse skeleton character point set L is referred to so as to calculaterb, with reference to mouse epidermis characteristic point set Lrs, target mouse Skeleton character point set LtbAnd target mouse epidermis characteristic point set Lts
The Characteristic points match of (3-2) based on TPS-RPM algorithms
The TPS-RPM algorithms proposed using Haili Chui[6]Program of increasing income, calculating refer to mouse skeleton character point set Lrb With target mouse skeleton character point set LtbBetween skeleton character point homography Cb, with reference to mouse epidermis characteristic point set LrsWith target mouse Epidermis characteristic point set LtsBetween epidermis characteristic point homography Cs
The reference mouse skeleton character point set LrbWith target mouse skeleton character point set LtbBetween skeleton character point correspond to square Battle array CbAnd with reference to mouse epidermis characteristic point set LrsWith target mouse epidermis characteristic point set LtsBetween epidermis characteristic point homography Cs To obscure homography, each element is the floating number between [0,1] in fuzzy homography, to describe between 2 points The power of degree of correspondence;
Wherein, mouse skeleton character point set L is referred in calculatingrbWith target mouse skeleton character point set LtbBetween skeleton character Point homography CbWhen, by target mouse skeleton character point set LtbIt is set as target point set, mouse skeleton character point set L will be referred torbIf It is set to point set subject to registration, and sets initial temperature coefficient and stopping criterion for iteration in simulated annealing used in TPS-RPM, leads to Cross the extreme point for iterating to calculate following energy equation:
In formula (1), cb,ijFor homography CbIn element;NrbTo refer to mouse skeleton character point set LrbIncluded characteristic point Number, NtbFor target mouse skeleton character point set LtbThe number of included point feature;lrb,iTo refer to mouse skeleton character point set Lrb In ith feature point coordinate, ltb,jFor target mouse skeleton character point set LtbIn j-th of characteristic point coordinate;F is thin plate Spline function;||Lf||2For the smoothness constraint to f, wherein L is differential operator;T be simulated annealing during be used for control mould Paste the temperature coefficient of degree of correspondence;λ is priori smoothness weights;γ is Lu Bang Control Sampled-Data weights;As temperature coefficient T gradually subtracts It is small, homography CbWith thin plate spline function f alternating iterations;Homography CbFinal iteration result be required skeleton character Point homography Cb
Mouse epidermis characteristic point set L is referred to calculatingrsWith target mouse epidermis characteristic point set LtsBetween epidermis characteristic point it is corresponding Matrix CsWhen, by target mouse epidermis characteristic point set LtsIt is set as target point set, mouse epidermis characteristic point set L will be referred torsIt is set as treating Registering point set, and initial temperature coefficient and stopping criterion for iteration in simulated annealing used in TPS-RPM are set, pass through iteration Calculate the extreme point of following energy equation:
In formula (2), cs,ijFor homography CsIn element;NrsTo refer to mouse epidermis characteristic point set LrsIncluded characteristic point Number, NtsFor target mouse epidermis characteristic point set LtsThe number of included point feature;lrs,iTo refer to mouse epidermis characteristic point set Lrs In i-th point of coordinate, lts,jFor target mouse epidermis characteristic point set LtsIn j-th point of coordinate;Remaining variables implication with Formula (1) is identical;As temperature coefficient T is gradually reduced, homography CsWith thin plate spline function f alternating iterations;Homography Cs's Final iteration result is required epidermis characteristic point homography Cs
Respectively by the skeleton character tried to achieve point homography CbAnd epidermis characteristic point homography CsAct on and refer to mouse bone Feature point set LrbAnd with reference to mouse epidermis characteristic point set Lrs, obtain the skeleton character point set L after preliminary registrationcbAnd after preliminary registration Epidermis characteristic point set Lcs
(3-3) preliminary registration local displacement matrix P structure
Utilize the result of the above-mentioned preliminary Characteristic points match tried to achieve, structure preliminary registration local displacement matrix P:
(3-4) preliminary registration mapping matrix McStructure
Above-mentioned matrix P is converted into three groups of 4 D data point set Px={ (x, y, z, △ x) }, Py={ (x, y, z, △ y) }, Pz={ (x, y, z, △ z) }, Δ x, Δ y and Δ z are respectively seen as to point (x, y, z) functional value, i.e. △ x=G1(x, y, z), △ y =G2(x, y, z), △ z=G3(x,y,z);
Using Multilevel B-splines fitting algorithm respectively to three groups of 4 D data point sets Px={ (x, y, z, △ x) }, Py ={ (x, y, z, △ y) }, Pz={ (x, y, z, △ z) } is fitted, for calculating the reference picture IrIn each pixel along Displacement on three directions of x, y, z;
Multilevel B-splines fitting algorithm is the three-dimensional data fitting algorithm that Seungyong Lee are proposed in the present invention[7] Extension on space-time.Utilize the Multilevel B-splines approximating method fitting data point set Px={ (x, y, z, △ x) } mistake Journey is described in detail as follows:
It is described to refer at many levels, using being overlying on reference picture IrOn one group of cube control grid Φ gradually encrypted0, Φ1,...,Φk,...,ΦhSuccessively to the 4 D data point set of iteration renewalB-spline fitting is carried out, and By required h layers fitting function sumAs final Multilevel B-splines fitting function;Wherein, △0ξ=△ x, △k+1ξ=△kξ-gk(x, y, z), gk(x, y, z) is kth layer B-spline fitting result;
The kth layer B-spline fit procedure is described as follows:
Assuming that kth layer control grid ΦkSize is Kx×Ky×Kz, then kth layer B-spline fitting function be shown below:
In formula (4), φk,(l+i,m+j,n+k)For positioned at control grid ΦkMiddle coordinate is The control node value of (l+i, m+j, n+k);l,m,n∈{0,1,2,3};Bl、BmAnd BnRespectively l, m, n rank B-spline base letter Number, wherein the expression formula of 0 to 3 rank B-spline function is described as follows:
In the formula (4), grid Φ is controlledkIn the value of each control node calculated by following two step:
(a) calculateIn each data point to controlling grid ΦkIn each control node value shadow Ring amount:
WithIn a data point p=(xp,yp,zp,△kξp) exemplified by be described as follows:
Data point p=(xp,yp,zp,△kξp) to controlling grid ΦkIn the influence matrix of each control node show as one Individual size is Kx×Ky×KzMatrix Ψp;It is easy to calculate, definition and matrix ΨpTwo matrix Γ of size identicalpWith Ωp; The matrix Ψp、ΓpWith ΩpMiddle coordinate is calculated by formula (6) respectively for the element of (l+i, m+j, n+k):
In formula (6), l, m, n ∈ { 0,1,2,3 };
In matrix Ψp、ΓpWith ΩpIn, except the coordinate is that (l+i, m+j, n+k) amounts to remaining position beyond 64 elements Put, ψp、γpWith ωpIt is 0;
(b) ask for controlling grid ΦkIn each control node value
It is comprehensiveIn each data point to controlling grid ΦkIn each control node value influence, Ask for grid Φ processedkIn each control node value;Control grid ΦkMiddle coordinate is the control node φ of (a, b, c)k,(a,b,c)Take It is worth and is:
In formula (7), γp,(a,b,c)、ωp,(a,b,c)、ψp,(a,b,c)Respectively Ψp、ΓpWith ΩpMiddle coordinate is the member of (a, b, c) Element value;
So far, kth layer B-spline fitting function gk(x, y, z) is established;Comprehensive each level fitting function, calculates Multilevel B Spline-fit functionUsing the Multilevel B-splines fitting function g (x, y, z) asked for, Calculate reference picture IrIn displacement of each pixel along x-axis;
Similarly, it is fitted 4 D data point set P using Multilevel B-splinesy={ (x, y, z, △ y) } and Pz={ (x, y, z, △ Z) }, so as to calculating reference picture IrIn displacement of each pixel along y, z-axis in, thus build reference picture IrIt is preliminary Registering mapping matrix Mc
(3-5) preliminary registration image IcStructure
Utilize the preliminary registration mapping matrix Mc, reversely solve preliminary registration image Ic;In structure preliminary registration image When, the assignment of gray scale uses tri-linear interpolation methods;
So far, preliminary registration mapping matrix Mc, with preliminary registration image IcAcquisition finishes.
Step 4: the fine registering mapping matrix M of structurefAnd fine registering image If
After preliminary registration, preliminary registration image IcWith target image ItStill there is bigger difference, in order to improve registration Precision is, it is necessary to carry out more fine registration.The main process of the step 4 is to preliminary registration image IcAnd target figure As ItFine registration is carried out, constructs control preliminary registration image IcProduce the fine registering mapping matrix M of non-rigidity deformationf, will MfAct on preliminary registration image Ic, construct fine registering image Icf.Specific implementation process is described in detail as follows:
(4-1) extracts preliminary registration characteristics of image point set again
Using image segmentation algorithm, preliminary registration image I is extractedcIn mouse skin region and bony areas, and profit Extract the profile of the mouse skin region and bony areas respectively with edge detection algorithm;By mouse skin and the wheel of bone Exterior feature is overlapped mutually, and is sampled using rectangular mesh, thus obtains one group of preliminary registration characteristics of image point set L'c
The Characteristic points match of (4-2) based on block matching method
Preliminary registration characteristics of image point set L' is asked for using block matching methodcIn each pair of the characteristic point on target image Position is answered, with L'cIn any point pl=(xp,yp,zp) exemplified by, the process description is as follows:
In preliminary registration image IcIn with coordinate (xp,yp,zp) centered on choose size be N1×N1×N1Cube neighborhood T, in target image ItIn with coordinate (xp,yp,zp) centered on choose N2×N2×N2Cube neighborhood S, wherein N2>N1;T is made For template, S is as region of search, and search and T have the subregion s of maximum similarity in S regions1, and by s1Central point pl′ As point plCorrespondence position;
By that analogy, preliminary registration characteristics of image point set L' is searched out successivelycMiddle each point is in target image ItOn corresponding position Put, thus construct fine registration features point set Lf
(4-3) fine registering local displacement matrix Q structure
Utilize preliminary registration characteristics of image point set L'cAnd fine registration features point set LfThe fine registering local displacement square of structure Battle array Q:
Q=[L'c,L'c-Lf]=[x, y, z, Δ x, Δ y, Δ z] (8)
(4-4) fine registering mapping matrix MfStructure
Above-mentioned matrix Q is converted into three groups of 4 D data point set Qx={ (x, y, z, △ x) }, Qy={ (x, y, z, △ y) }, Qz={ (x, y, z, △ z) };It is consistent with Multilevel B-splines fit procedure in step 3, respectively to three groups of 4 D datas point Collect Qx={ (x, y, z, △ x) }, Qy={ (x, y, z, △ y) } carries out Multilevel B-splines fitting, tentatively matches somebody with somebody so as to calculate respectively Quasi- image IcIn displacement of each pixel along three axial directions of x, y, z, thus build preliminary registration image IcIt is fine registration mapping Matrix Mf
So far, fine registering mapping matrix MfAcquisition finishes.
(4-5) fine registering image IfStructure (optional)
After this process, it can select fine registering mapping matrix MfAct on preliminary registration image Ic, construct fine Registering image If;When building preliminary registration image, the assignment of gray scale equally uses tri-linear interpolation methods;It is final by judging Registering image (i.e. fine registering image If) and target image ItSimilarity with evaluate registration significant degree;
Step 5: structure target mouse anatomical structure collection of illustrative plates:
In non-rigidity process of image registration, pass through preliminary registration mapping matrix McAnd fine registering mapping matrix MfControl Make and use, reference picture IrIt is deformed into fine registering image If, fine registering image IfWith target image ItWith higher similar Degree.Therefore, by preliminary registration mapping matrix McAnd fine registering mapping matrix MfThe anatomical structure collection of illustrative plates with reference to mouse is acted on successively Ar, then can the approximate anatomical structure collection of illustrative plates for building target mouse.However, built on reference to anatomical structure collection of illustrative plates under physical coordinates system (in units of mm), and above-mentioned two registration mapping matrix builds under pixel coordinate system (in units of pixel), therefore herein During need to ArCarry out coordinate transform.Specific implementation process is described in detail as follows:
(5-1) is by described with reference to mouse anatomical structure collection of illustrative plates ArProjection obtains the ginseng under pixel coordinate system to pixel coordinate system Examine mouse anatomical structure collection of illustrative plates
(5-2) is in order successively by preliminary registration mapping matrix McAnd fine registering mapping matrix MfAct on the pixel Reference mouse anatomical structure collection of illustrative plates under coordinate systemMake the reference mouse anatomical structure collection of illustrative plates under the pixel coordinate systemProduce Deformed with step 3 preliminary registration and the fine registration process identical of step 4, obtain the registering mouse dissection knot under pixel coordinate system Composition is composed
(5-3) is by the registering mouse anatomical structure collection of illustrative plates under the pixel coordinate systemUnder projection to physical coordinates system, obtain Registering mouse anatomical structure collection of illustrative plates A under physical coordinates systemf, the registering mouse anatomical structure collection of illustrative plates A under the physical coordinates systemfAs Target mouse anatomical structure collection of illustrative plates.
Although above in conjunction with accompanying drawing, invention has been described, and the invention is not limited in above-mentioned specific implementation Mode, above-mentioned embodiment is only schematical, rather than restricted, and one of ordinary skill in the art is at this Under the enlightenment of invention, without deviating from the spirit of the invention, many variations can also be made, these belong to the present invention's Within protection.

Claims (1)

1. a kind of target mouse anatomical structure collection of illustrative plates acquisition methods based on the non-rigidity registration of XCT images, it is characterised in that its is basic Step is:First, it is with reference to mouse anatomical structure collection of illustrative plates, while by XCT corresponding to Digimouse by Digimouse model specifications Image setting is reference picture;Secondly, XCT imagings are carried out to target mouse to obtain target image and pre-process;Then, utilize Non- rigidity image registration techniques build reference picture to the registering mapping matrix of target image;Finally, registering mapping matrix is made For referring to mouse anatomical structure collection of illustrative plates, target mouse anatomical structure collection of illustrative plates is constructed;Comprise the following steps that:
Step 1: setting refers to mouse anatomical structure collection of illustrative plates and reference picture:
It is with reference to mouse anatomical structure collection of illustrative plates A by Digimouse configuration settingsr, it is ginseng by the XCT image settings of the Digimouse Examine image Ir
Step 2: the acquisition and pretreatment of target mouse XCT images:
Whole body imaging is carried out to target mouse using XCT equipment, obtains the XCT images of target mouse;Target mouse XCT images are imitated Conversion is penetrated, the head of mouse and the direction at back in target mouse XCT images are converted into and reference picture IrIt is identical;After converting XCT images as target image It
Step 3: structure preliminary registration mapping matrix McWith preliminary registration image Ic
It is partitioned into reference picture I respectively using image segmentation algorithmrIn mouse bony areas and cuticle region, and target image ItIn mouse bony areas and cuticle region;
Extract above-mentioned reference picture I respectively using edge detection algorithmrMiddle mouse bony areas and cuticle region and target image It In mouse bony areas and cuticle region amount to the boundary profile in four regions, the boundary profile in four regions is carried out Equiprobability is sampled, and mouse skeleton character point set L is referred to so as to calculaterb, with reference to mouse epidermis characteristic point set Lrs, target mouse bone it is special Levy point set LtbAnd target mouse epidermis characteristic point set Lts
Using the program of increasing income of TPS-RPM algorithms, calculating refers to mouse skeleton character point set LrbWith target mouse skeleton character point set Ltb Between skeleton character point homography Cb;When calculating, by target mouse skeleton character point set LtbIt is set as target point set, will joins Examine mouse skeleton character point set LrbIt is set as point set subject to registration, and sets the initial temperature in simulated annealing used in TPS-RPM Coefficient and stopping criterion for iteration, by iterating to calculate out skeleton character point homography Cb
Using the program of increasing income of TPS-RPM algorithms, calculating refers to mouse epidermis characteristic point set LrsWith target mouse epidermis characteristic point set Lts Between epidermis characteristic point homography Cs;When calculating, by target mouse epidermis characteristic point set LtsIt is set as target point set, will joins Examine mouse epidermis characteristic point set LrsIt is set as point set subject to registration, and sets the initial temperature in simulated annealing used in TPS-RPM Coefficient and stopping criterion for iteration, by iterating to calculate out epidermis characteristic point homography Cs
Skeleton character point homography C is asked for using the program of increasing income of TPS-RPM algorithmsbWith epidermis characteristic point homography Cs's Detailed process is as follows:
The reference mouse skeleton character point set LrbWith target mouse skeleton character point set LtbBetween skeleton character point homography Cb And with reference to mouse epidermis characteristic point set LrsWith target mouse epidermis characteristic point set LtsBetween epidermis characteristic point homography CsIt is mould Homography is pasted, each element is the floating number between [0,1] in fuzzy homography, to corresponding between describing at 2 points The power of degree;
The TPS-RPM algorithms refer to mouse skeleton character point set L in calculatingrbWith target mouse skeleton character point set LtbBetween bone Characteristic point homography CbWhen, the minimum point of energy equation is calculated as follows using simulated annealing:
<mrow> <mi>E</mi> <mrow> <mo>(</mo> <msub> <mi>C</mi> <mi>b</mi> </msub> <mo>,</mo> <mi>f</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>r</mi> <mi>b</mi> </mrow> </msub> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>t</mi> <mi>b</mi> </mrow> </msub> </munderover> <msub> <mi>c</mi> <mrow> <mi>b</mi> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>|</mo> <mo>|</mo> <msub> <mi>l</mi> <mrow> <mi>r</mi> <mi>b</mi> <mo>,</mo> <mi>i</mi> </mrow> </msub> <mo>-</mo> <mi>f</mi> <mrow> <mo>(</mo> <msub> <mi>l</mi> <mrow> <mi>t</mi> <mi>b</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>+</mo> <mi>&amp;lambda;</mi> <mo>|</mo> <mo>|</mo> <mi>L</mi> <mi>f</mi> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>+</mo> <mi>T</mi> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>r</mi> <mi>b</mi> </mrow> </msub> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>t</mi> <mi>b</mi> </mrow> </msub> </munderover> <msub> <mi>c</mi> <mrow> <mi>b</mi> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mi>log</mi> <mi> </mi> <msub> <mi>c</mi> <mrow> <mi>b</mi> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>-</mo> <mi>&amp;gamma;</mi> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>r</mi> <mi>b</mi> </mrow> </msub> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>t</mi> <mi>b</mi> </mrow> </msub> </munderover> <msub> <mi>c</mi> <mrow> <mi>b</mi> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
In formula (1), cb,ijFor homography CbIn element;NrbTo refer to mouse skeleton character point set LrbOf included characteristic point Number, NtbFor target mouse skeleton character point set LtbThe number of included point feature;lrb,iTo refer to mouse skeleton character point set LrbIn The coordinate of ith feature point, ltb,jFor target mouse skeleton character point set LtbIn j-th of characteristic point coordinate;F is thin plate spline Function;||Lf||2For the smoothness constraint to f, wherein L is differential operator;T is to be used to control fuzzy pair during simulated annealing Answer the temperature coefficient of degree;λ is priori smoothness weights;γ is Lu Bang Control Sampled-Data weights;As temperature coefficient T is gradually reduced, Homography CbWith thin plate spline function f alternating iterations;Homography CbFinal iteration result be required skeleton character point pair Answer Matrix Cb
The TPS-RPM algorithms refer to mouse epidermis characteristic point set L in calculatingrsWith target mouse epidermis characteristic point set LtsBetween epidermis Characteristic point homography CsWhen, the minimum point of energy equation is calculated as follows using simulated annealing:
<mrow> <mi>E</mi> <mrow> <mo>(</mo> <msub> <mi>C</mi> <mi>s</mi> </msub> <mo>,</mo> <mi>f</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>r</mi> <mi>s</mi> </mrow> </msub> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>t</mi> <mi>s</mi> </mrow> </msub> </munderover> <msub> <mi>c</mi> <mrow> <mi>s</mi> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>|</mo> <mo>|</mo> <msub> <mi>l</mi> <mrow> <mi>r</mi> <mi>s</mi> <mo>,</mo> <mi>i</mi> </mrow> </msub> <mo>-</mo> <mi>f</mi> <mrow> <mo>(</mo> <msub> <mi>l</mi> <mrow> <mi>t</mi> <mi>s</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>+</mo> <mi>&amp;lambda;</mi> <mo>|</mo> <mo>|</mo> <mi>L</mi> <mi>f</mi> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>+</mo> <mi>T</mi> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>r</mi> <mi>s</mi> </mrow> </msub> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>t</mi> <mi>s</mi> </mrow> </msub> </munderover> <msub> <mi>c</mi> <mrow> <mi>s</mi> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mi>log</mi> <mi> </mi> <msub> <mi>c</mi> <mrow> <mi>s</mi> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>-</mo> <mi>&amp;gamma;</mi> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>r</mi> <mi>s</mi> </mrow> </msub> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mrow> <mi>t</mi> <mi>s</mi> </mrow> </msub> </munderover> <msub> <mi>c</mi> <mrow> <mi>s</mi> <mo>,</mo> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>
In formula (2), cs,ijFor homography CsIn element;NrsTo refer to mouse epidermis characteristic point set LrsOf included characteristic point Number, NtsFor target mouse epidermis characteristic point set LtsThe number of included point feature;lrs,iTo refer to mouse epidermis characteristic point set LrsIn I-th point of coordinate, lts,jFor target mouse epidermis characteristic point set LtsIn j-th point of coordinate;Remaining variables implication and formula (1) In it is identical;As temperature coefficient T is gradually reduced, homography CsWith thin plate spline function f alternating iterations;Homography CsMost Whole iteration result is required epidermis characteristic point homography Cs
Respectively by the skeleton character tried to achieve point homography CbAnd epidermis characteristic point homography CsAct on and refer to mouse skeleton character Point set LrbAnd with reference to mouse epidermis characteristic point set Lrs, obtain the skeleton character point set L after preliminary registrationcbAnd the table after preliminary registration Skin feature point set Lcs
Utilize the result of the above-mentioned preliminary Characteristic points match tried to achieve, structure preliminary registration local displacement matrix P:
<mrow> <mi>P</mi> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>L</mi> <mrow> <mi>r</mi> <mi>b</mi> </mrow> </msub> </mtd> <mtd> <mrow> <msub> <mi>L</mi> <mrow> <mi>r</mi> <mi>b</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>L</mi> <mrow> <mi>c</mi> <mi>b</mi> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <msub> <mi>L</mi> <mrow> <mi>r</mi> <mi>s</mi> </mrow> </msub> </mtd> <mtd> <mrow> <msub> <mi>L</mi> <mrow> <mi>r</mi> <mi>s</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>L</mi> <mrow> <mi>c</mi> <mi>s</mi> </mrow> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <mo>&amp;lsqb;</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>,</mo> <mi>&amp;Delta;</mi> <mi>x</mi> <mo>,</mo> <mi>&amp;Delta;</mi> <mi>y</mi> <mo>,</mo> <mi>&amp;Delta;</mi> <mi>z</mi> <mo>&amp;rsqb;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>
Above-mentioned matrix P is converted into three groups of 4 D data point set Px=(x, y, z, Δ x) }, Py=(x, y, z, Δ y) }, Pz= { (x, y, z, Δ z) }, Δ x, Δ y and Δ z are respectively seen as to point (x, y, z) functional value, i.e. Δ x=G1(x, y, z), Δ y=G2 (x, y, z), Δ z=G3(x,y,z);
Using Multilevel B-splines fitting algorithm respectively to three groups of 4 D data point sets Px=(x, y, z, Δ x) }, Py= (x, y, z, Δ y) }, Pz=(x, y, z, Δ z) } it is fitted, for calculating the reference picture IrIn each pixel along x, Y, the displacement on tri- directions of z;
Utilize the Multilevel B-splines approximating method fitting data point set PxThe process of={ (x, y, z, Δ x) } is described in detail such as Under:
It is described to refer at many levels, using being overlying on reference picture IrOn one group of cube control grid Φ gradually encrypted0, Φ1,...,Φk,...,ΦhSuccessively to the 4 D data point set of iteration renewalB-spline fitting is carried out, and By required h layers fitting function sumAs final Multilevel B-splines fitting function;Wherein, Δ0ξ=Δ x, Δk+1ξ=Δkξ-gk(x, y, z), gk(x, y, z) is kth layer B-spline fitting result;
The kth layer B-spline fit procedure is described as follows:
Assuming that kth layer control grid ΦkSize is Kx×Ky×Kz, then kth layer B-spline fitting function be shown below:
<mrow> <msub> <mi>g</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>0</mn> </mrow> <mn>3</mn> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>0</mn> </mrow> <mn>3</mn> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>0</mn> </mrow> <mn>3</mn> </munderover> <msub> <mi>B</mi> <mi>l</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>)</mo> </mrow> <msub> <mi>B</mi> <mi>m</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> <mo>)</mo> </mrow> <msub> <mi>B</mi> <mi>n</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mi>z</mi> </msub> <mo>)</mo> </mrow> <msub> <mi>&amp;phi;</mi> <mrow> <mi>k</mi> <mo>,</mo> <mrow> <mo>(</mo> <mi>l</mi> <mo>+</mo> <mi>i</mi> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mi>j</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>
In formula (4), φk,(l+i,m+j,n+k)For positioned at control grid ΦkMiddle coordinate is The control node value of (l+i, m+j, n+k);l,m,n∈{0,1,2,3};Bl、BmAnd BnRespectively l, m, n rank B-spline base letter Number, wherein the expression formula of 0 to 3 rank B-spline function is described as follows:
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>B</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>&amp;delta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <mi>&amp;delta;</mi> <mo>)</mo> </mrow> <mn>3</mn> </msup> <mo>/</mo> <mn>6</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>B</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>&amp;delta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mrow> <mo>(</mo> <mn>3</mn> <msup> <mi>&amp;delta;</mi> <mn>3</mn> </msup> <mo>-</mo> <mn>6</mn> <msup> <mi>&amp;delta;</mi> <mn>2</mn> </msup> <mo>+</mo> <mn>4</mn> <mo>)</mo> </mrow> <mo>/</mo> <mn>6</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>B</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>&amp;delta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mrow> <mo>(</mo> <mo>-</mo> <mn>3</mn> <msup> <mi>&amp;delta;</mi> <mn>3</mn> </msup> <mo>+</mo> <mn>3</mn> <msup> <mi>&amp;delta;</mi> <mn>2</mn> </msup> <mo>+</mo> <mn>3</mn> <mi>&amp;delta;</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>/</mo> <mn>6</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>B</mi> <mn>3</mn> </msub> <mrow> <mo>(</mo> <mi>&amp;delta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>&amp;delta;</mi> <mn>3</mn> </msup> <mo>/</mo> <mn>6</mn> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>
In the formula (4), grid Φ is controlledkIn the value of each control node calculated by following two step:
(a) calculateIn each data point to controlling grid ΦkIn each control node value influence amount:
WithIn a data point p=(xp,yp,zpkξp) exemplified by be described as follows:
Data point p=(xp,yp,zpkξp) to controlling grid ΦkIn the influence matrix of each control node show as a chi Very little is Kx×Ky×KzMatrix Ψp;It is easy to calculate, definition and matrix ΨpTwo matrix Γ of size identicalpWith Ωp;It is described Matrix Ψp、ΓpWith ΩpMiddle coordinate is calculated by formula (6) respectively for the element of (l+i, m+j, n+k):
<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>&amp;psi;</mi> <mrow> <mi>p</mi> <mo>,</mo> <mrow> <mo>(</mo> <mi>l</mi> <mo>+</mo> <mi>i</mi> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mi>j</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <msub> <mi>&amp;gamma;</mi> <mrow> <mi>p</mi> <mo>,</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>+</mo> <mi>i</mi> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mi>j</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </msub> <msup> <mi>&amp;Delta;</mi> <mi>k</mi> </msup> <msub> <mi>&amp;xi;</mi> <mi>p</mi> </msub> </mrow> <msub> <mi>&amp;omega;</mi> <mrow> <mi>p</mi> <mo>,</mo> <mrow> <mo>(</mo> <mi>l</mi> <mo>+</mo> <mi>i</mi> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mi>j</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </msub> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>&amp;gamma;</mi> <mrow> <mi>p</mi> <mo>,</mo> <mrow> <mo>(</mo> <mi>l</mi> <mo>+</mo> <mi>i</mi> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mi>j</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </msub> <mo>=</mo> <msub> <mi>B</mi> <mi>l</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mrow> <mi>x</mi> <mi>p</mi> </mrow> </msub> <mo>)</mo> </mrow> <msub> <mi>B</mi> <mi>m</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mrow> <mi>y</mi> <mi>p</mi> </mrow> </msub> <mo>)</mo> </mrow> <msub> <mi>B</mi> <mi>n</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mrow> <mi>z</mi> <mi>p</mi> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>&amp;omega;</mi> <mrow> <mi>p</mi> <mo>,</mo> <mrow> <mo>(</mo> <mi>l</mi> <mo>+</mo> <mi>i</mi> <mo>,</mo> <mi>m</mi> <mo>+</mo> <mi>j</mi> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </msub> <mo>=</mo> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>0</mn> </mrow> <mn>3</mn> </msubsup> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>0</mn> </mrow> <mn>3</mn> </msubsup> <msubsup> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>0</mn> </mrow> <mn>3</mn> </msubsup> <msup> <mrow> <mo>&amp;lsqb;</mo> <msub> <mi>B</mi> <mi>l</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mrow> <mi>x</mi> <mi>p</mi> </mrow> </msub> <mo>)</mo> </mrow> <msub> <mi>B</mi> <mi>m</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mrow> <mi>y</mi> <mi>p</mi> </mrow> </msub> <mo>)</mo> </mrow> <msub> <mi>B</mi> <mi>n</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mrow> <mi>z</mi> <mi>p</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> <mn>2</mn> </msup> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>
In formula (6), l, m, n ∈ { 0,1,2,3 };
In matrix Ψp、ΓpWith ΩpIn, except the coordinate is that (l+i, m+j, n+k) amounts to remaining position beyond 64 elements, ψp、 γpWith ωpIt is 0;
(b) ask for controlling grid ΦkIn each control node value
It is comprehensiveIn each data point to controlling grid ΦkIn each control node value influence, ask for Grid Φ processedkIn each control node value;Control grid ΦkMiddle coordinate is the control node φ of (a, b, c)k,(a,b,c)Value For:
In formula (7), γp,(a,b,c)、ωp,(a,b,c)、ψp,(a,b,c)Respectively Ψp、ΓpWith ΩpMiddle coordinate is the element of (a, b, c) Value;
So far, kth layer B-spline fitting function gk(x, y, z) is established;Comprehensive each level fitting function, calculates Multilevel B-splines Fitting functionUsing the Multilevel B-splines fitting function g (x, y, z) asked for, calculate Go out reference picture IrIn displacement of each pixel along x-axis;
Similarly, it is fitted 4 D data point set P using Multilevel B-splinesy=(x, y, z, Δ y) } and Pz=(x, y, z, Δ z) }, So as to calculate reference picture IrIn displacement of each pixel along y, z-axis in, thus build reference picture IrPreliminary registration Mapping matrix Mc
Utilize the preliminary registration mapping matrix Mc, reversely solve preliminary registration image Ic;When building preliminary registration image, The assignment of gray scale uses tri-linear interpolation methods;
Step 4: the fine registering mapping matrix M of structuref
Using image segmentation algorithm, preliminary registration image I is extractedcIn mouse skin region and bony areas, and utilize edge Detection algorithm extracts the profile of the mouse skin region and bony areas respectively;The profile of mouse skin and bone is mutual Superposition, and sampled using rectangular mesh, thus obtain one group of preliminary registration characteristics of image point set L'c
Preliminary registration characteristics of image point set L' is asked for using block matching methodcIn each corresponding position of the characteristic point on target image Put, with L'cIn any point pl=(xp,yp,zp) exemplified by, the process description is as follows:
In preliminary registration image IcIn with coordinate (xp,yp,zp) centered on choose size be N1×N1×N1Cube neighborhood T, in mesh Logo image ItIn with coordinate (xp,yp,zp) centered on choose N2×N2×N2Cube neighborhood S, wherein N2>N1;Using T as template, S is as region of search, and search and T have the subregion s of maximum similarity in S regions1, and by s1Central point pl' as point pl Correspondence position;
By that analogy, preliminary registration characteristics of image point set L' is searched out successivelycMiddle each point is in target image ItOn correspondence position, Thus fine registration features point set L is constructedf
Utilize preliminary registration characteristics of image point set L'cAnd fine registration features point set LfThe fine registering local displacement matrix Q of structure:
Q=[L'c,L'c-Lf]=[x, y, z, Δ x, Δ y, Δ z] (8)
Above-mentioned matrix Q is converted into three groups of 4 D data point set Qx=(x, y, z, Δ x) }, Qy=(x, y, z, Δ y) }, Qz= {(x,y,z,Δz)};It is consistent with Multilevel B-splines fit procedure in step 3, respectively to three groups of 4 D data point sets Qx =(x, y, z, Δ x) }, Qy=(x, y, z, Δ y) } Multilevel B-splines fitting is carried out, so as to calculate preliminary registration figure respectively As IcIn displacement of each pixel along three axial directions of x, y, z, thus build preliminary registration image IcFine registering mapping matrix Mf
Step 5: structure target mouse anatomical structure collection of illustrative plates:
By described with reference to mouse anatomical structure collection of illustrative plates ArProjection obtains the reference mouse dissection knot under pixel coordinate system to pixel coordinate system Composition is composedIn order successively by preliminary registration mapping matrix McAnd fine registering mapping matrix MfAct on the pixel coordinate Reference mouse anatomical structure collection of illustrative plates under systemMake the reference mouse anatomical structure collection of illustrative plates under the pixel coordinate systemProduce and walk Rapid three preliminary registration and the fine registration process identical deformation of step 4, obtain the registering mouse anatomical structure figure under pixel coordinate system SpectrumBy the registering mouse anatomical structure collection of illustrative plates under the pixel coordinate systemUnder projection to physical coordinates system, physical coordinates are obtained Registering mouse anatomical structure collection of illustrative plates A under systemf, the registering mouse anatomical structure collection of illustrative plates A under the physical coordinates systemfAs target mouse solution Cut open structure collection of illustrative plates.
CN201510259870.5A 2015-05-20 2015-05-20 Target mouse anatomical structure collection of illustrative plates acquisition methods based on the non-rigidity registration of XCT images Expired - Fee Related CN104867104B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201510259870.5A CN104867104B (en) 2015-05-20 2015-05-20 Target mouse anatomical structure collection of illustrative plates acquisition methods based on the non-rigidity registration of XCT images

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201510259870.5A CN104867104B (en) 2015-05-20 2015-05-20 Target mouse anatomical structure collection of illustrative plates acquisition methods based on the non-rigidity registration of XCT images

Publications (2)

Publication Number Publication Date
CN104867104A CN104867104A (en) 2015-08-26
CN104867104B true CN104867104B (en) 2017-12-15

Family

ID=53912921

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201510259870.5A Expired - Fee Related CN104867104B (en) 2015-05-20 2015-05-20 Target mouse anatomical structure collection of illustrative plates acquisition methods based on the non-rigidity registration of XCT images

Country Status (1)

Country Link
CN (1) CN104867104B (en)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105212936A (en) * 2015-09-09 2016-01-06 首都医科大学宣武医院 A kind of generation method of brain template
CN105427237B (en) * 2015-11-25 2018-04-24 长春乙天科技有限公司 A kind of steel mesh image registration of large format optical measuring system and detection method
CN112967236B (en) * 2018-12-29 2024-02-27 上海联影智能医疗科技有限公司 Image registration method, device, computer equipment and storage medium
CN110322491B (en) * 2019-06-11 2022-03-04 大连理工大学 Algorithm for registering deformable mouse whole-body atlas and mouse image
CN111091567B (en) * 2020-03-23 2020-06-23 南京景三医疗科技有限公司 Medical image registration method, medical device and storage medium
CN112116642A (en) * 2020-09-27 2020-12-22 上海联影医疗科技股份有限公司 Registration method of medical image, electronic device and storage medium
CN112995528B (en) * 2021-05-06 2021-09-21 中国工程物理研究院流体物理研究所 Method for registering images among channels of photoelectric framing camera
CN113298856B (en) * 2021-05-28 2023-10-20 上海联影医疗科技股份有限公司 Image registration method, device, equipment and medium
CN114742869B (en) * 2022-06-15 2022-08-16 西安交通大学医学院第一附属医院 Brain neurosurgery registration method based on pattern recognition and electronic equipment

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103077522A (en) * 2013-01-22 2013-05-01 南方医科大学 Improved registration method for image with obvious local deformation
CN104599268A (en) * 2015-01-06 2015-05-06 广州医科大学附属肿瘤医院 Local area accurate deformation registration algorithm combining point registration

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9524552B2 (en) * 2011-08-03 2016-12-20 The Regents Of The University Of California 2D/3D registration of a digital mouse atlas with X-ray projection images and optical camera photos

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103077522A (en) * 2013-01-22 2013-05-01 南方医科大学 Improved registration method for image with obvious local deformation
CN104599268A (en) * 2015-01-06 2015-05-06 广州医科大学附属肿瘤医院 Local area accurate deformation registration algorithm combining point registration

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
"一种多模态医学图像数据融合方法与应用";高峰等;《中国医疗设备》;20131231;第28卷(第5期);第164-167页 *
"基于区域标识的‘粗粒度’扩散光学层析成像方法稳态测量模式";赵会娟、高峰等;《光学学报》;20130228;第33卷(第2期);第1-8页 *
"小鼠荧光层析成像***及数据提取扩展方法";赵会娟等;《光子学报》;20141030;第43卷(第10期);第1-7页 *

Also Published As

Publication number Publication date
CN104867104A (en) 2015-08-26

Similar Documents

Publication Publication Date Title
CN104867104B (en) Target mouse anatomical structure collection of illustrative plates acquisition methods based on the non-rigidity registration of XCT images
Cao et al. Deformable image registration using a cue-aware deep regression network
Claes et al. Computerized craniofacial reconstruction: conceptual framework and review
CN101339670B (en) Computer auxiliary three-dimensional craniofacial rejuvenation method
CN110459301A (en) Brain neuroblastoma surgical navigation method for registering based on thermodynamic chart and facial key point
CN106204561A (en) Prostate multi-modality images non-rigid registration method based on mixed model
CN104050666B (en) Brain MR image method for registering based on segmentation
CN102903103B (en) Migratory active contour model based stomach CT (computerized tomography) sequence image segmentation method
CN113570627B (en) Training method of deep learning segmentation network and medical image segmentation method
CN107507195A (en) The multi-modal nasopharyngeal carcinoma image partition methods of PET CT based on hypergraph model
CN109118455B (en) Ancient human skull craniofacial interactive restoration method based on modern soft tissue distribution
Xie et al. Feature‐based rectal contour propagation from planning CT to cone beam CT
CN107220965A (en) A kind of image partition method and system
Desvignes et al. 3D semi-landmarks based statistical face reconstruction
CN114792326A (en) Surgical navigation point cloud segmentation and registration method based on structured light
Wang et al. Deformable torso phantoms of Chinese adults for personalized anatomy modelling
CN106709867A (en) Medical image registration method based on improved SURF and improved mutual information
Markel et al. A 4D biomechanical lung phantom for joint segmentation/registration evaluation
Wang et al. Using optimal transport to improve spherical harmonic quantification of complex biological shapes
CN111798500B (en) Differential synblast non-rigid registration algorithm based on hierarchical neighborhood spectral features
CN110322491B (en) Algorithm for registering deformable mouse whole-body atlas and mouse image
CN112598669A (en) Lung lobe segmentation method based on digital human technology
Villard et al. ISACHI: integrated segmentation and alignment correction for heart images
Sun et al. Stepwise local synthetic pseudo-CT imaging based on anatomical semantic guidance
CN114913149B (en) Head deformable statistical map construction method based on CT image

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
EXSB Decision made by sipo to initiate substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20171215

Termination date: 20200520