WO2018000652A1 - 一种非刚性多模医学图像的配准方法及系统 - Google Patents

一种非刚性多模医学图像的配准方法及系统 Download PDF

Info

Publication number
WO2018000652A1
WO2018000652A1 PCT/CN2016/101547 CN2016101547W WO2018000652A1 WO 2018000652 A1 WO2018000652 A1 WO 2018000652A1 CN 2016101547 W CN2016101547 W CN 2016101547W WO 2018000652 A1 WO2018000652 A1 WO 2018000652A1
Authority
WO
WIPO (PCT)
Prior art keywords
image
registration
order
floating
local feature
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.)
Ceased
Application number
PCT/CN2016/101547
Other languages
English (en)
French (fr)
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.)
Huazhong University of Science and Technology
Original Assignee
Huazhong University of Science and Technology
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 Huazhong University of Science and Technology filed Critical Huazhong University of Science and Technology
Priority to US16/094,473 priority Critical patent/US10853941B2/en
Publication of WO2018000652A1 publication Critical patent/WO2018000652A1/zh
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/14—Transformations for image registration, e.g. adjusting or mapping for alignment of images
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • G06T7/0014—Biomedical image inspection using an image reference approach
    • G06T7/0016—Biomedical image inspection using an image reference approach involving temporal comparison
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T12/00—Tomographic reconstruction from projections
    • G06T12/10—Image preprocessing, e.g. calibration, positioning of sources or scatter correction
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T12/00—Tomographic reconstruction from projections
    • G06T12/20—Inverse problem, i.e. transformations from projection space into object space
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/60—Rotation of whole images or parts thereof
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H30/00—ICT specially adapted for the handling or processing of medical images
    • G16H30/40—ICT specially adapted for the handling or processing of medical images for processing medical images, e.g. editing
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H50/00—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics
    • G16H50/20—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics for computer-aided diagnosis, e.g. based on medical expert systems
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H50/00—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics
    • G16H50/50—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics for simulation or modelling of medical disorders
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10088—Magnetic resonance imaging [MRI]
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10104—Positron emission tomography [PET]
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10132—Ultrasound image
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20048—Transform domain processing
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20212—Image combination
    • G06T2207/20221—Image fusion; Image merging
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30016—Brain

Definitions

  • the invention belongs to the field of image registration in image processing and analysis, and more particularly to a registration method and system for non-rigid multi-mode medical images.
  • Medical imaging technology is an important part of modern medicine and has revolutionary significance for the diagnosis and treatment of diseases.
  • the existing imaging technology is developing rapidly, due to the limitation of its imaging principle, it can only provide single and limited information.
  • ultrasound, CT and MRI can provide anatomical information of organs, but can not provide their functional information; PET can Better provide metabolic information, but can not clearly provide morphological information of organs.
  • doctors often need to combine information from different modal images to improve the accuracy of medical imaging diagnosis.
  • multi-mode image fusion multi-mode image registration is essential.
  • image registration is to find the correspondence between reference images and floating images.
  • image registration consists of three main parts: deformation model, similarity measure and optimization method.
  • deformation model has been a problem in the field of image registration due to factors such as grayscale distortion, patient breathing and body position movement.
  • the methods for solving the problem of non-rigid multi-mode image registration are generally divided into two categories: the first type is a registration method based on mutual information measure, but such methods usually do not consider the local feature structure of the image, and the calculation takes time. And it is easy to fall into local extremes, resulting in inaccurate registration results.
  • the second method simplifies multimodal image registration to single-mode image registration by image structure characterization, such as through entropy graph, Weber Local Descriptor (WLD), and modal independent neighborhood descriptors ( Modality Independent Neighborhood Descriptor (MIND) and other features characterize the image structure, and then use the Sum of Squared Difference (SSD) as the registration measure to achieve image registration.
  • image structure characterization such as through entropy graph, Weber Local Descriptor (WLD), and modal independent neighborhood descriptors ( Modality Independent Neighborhood Descriptor (MIND) and other features characterize the image structure, and then use the Sum of Squared Difference (SSD) as the registration measure to achieve image registration.
  • WLD Web
  • the method based on entropy map is The entropy value of each image block is calculated by estimating the gray probability density function of the image block, thereby obtaining the entropy map of the whole image, and the WLD-based method is to describe the local structure by using the Laplacian operation of the image. feature.
  • the MIND-based registration method evaluates the image self-similarity by using the Euclidean distance between image blocks, and uses the sum of similarities between different image blocks to describe the local structural features of the image.
  • this method only considers the self-similarity of the image.
  • the translation invariance between image blocks neglects the rotation characteristics that may exist between image blocks, so in the case where there is rotation between image blocks, this method is difficult to provide accurate registration results.
  • the present invention provides a non-rigid multi-modal medical image accurate registration method and system based on image self-similarity, thereby constructing an image by using Zernike moments through self-similarity of images. Describe features to achieve registration of non-rigid multimodal medical images.
  • a registration method for a non-rigid multi-mode medical image comprising:
  • the registration method specifically includes the following steps:
  • Step 1 obtains a local feature descriptor ZMLD of the reference image based on the 0th-order 0-weight Zernike moment Z A 00 (i) of the pixel point i in the reference image I A (i) and the 1st-order 1 Zernike moment Z A 11 (i) A (i); wherein i is an integer from 1 to M, and M is a reference image and a size of the floating image;
  • Step 2 obtains local features of the floating image according to the 0th-order 0-heavy Zernike moment Z B 00 (i) of the pixel i at the same position in the floating image I B (i) and the 1st-order 1 Zernike moment Z B 11 (i) Descriptor ZMLD B (i);
  • Step 3 based on the local feature description sub-ZMLD A (i) and ZMLD B (i) of the reference image and the floating image, establish an objective function g(T); obtain a transformation parameter according to the objective function, and transform the floating image according to the transformation parameter, And performing interpolation processing on the transformed floating image to obtain the registration image.
  • said step 2 comprises the following substeps:
  • Step 2-1 In the reference image, obtain a similar distance D A 00 (i) corresponding to the 0-order 0-heavy Zernike moment of the pixel point i and the remaining pixel points in the image block centered thereon, and the corresponding 1st-order 1-weight Zernike
  • the similarity distance of the moment D A 11 (i) at the same time, in the floating image, the similar distance D B 00 (i) corresponding to the 0th order 0 heavy Zernike moment of the pixel point i and the remaining pixel points in the image block centered thereon is obtained.
  • Step 2-2 to obtain a reference image local feature local feature descriptor ZMLD A (i) and the floating image descriptors ZMLD B (i);
  • h 00 A (i), h 11 A (i), h 00 B (i), and h 00 B (i) are attenuation coefficients.
  • MED( ⁇ ) is a median operator
  • c 1 and c 2 are constants of 0.5 to 1
  • j is any pixel of an image block centered on pixel i.
  • the step 3 specifically includes the following sub-steps:
  • Step 3-2 transforming the local feature descriptor ZMLD B (i) of the floating image according to the transformation parameter T ⁇ , performing interpolation processing on the transformed local feature descriptor, and updating the original local component by the interpolation process local feature descriptor
  • the feature descriptor ZMLD B (i), ⁇ ⁇ +1; iteratively solves the objective function g(T ⁇ ) to obtain a transformation parameter T ⁇ ;
  • Step 3-3 If the number of iterations ⁇ is greater than or equal to the threshold number of iterations ⁇ , and g(T ⁇ ) ⁇ g(T ⁇ -1 ), the floating image is transformed according to the transformation parameter T ⁇ , and the transformed floating image is performed. Interpolation process to obtain the registration image; otherwise return to step 3-2.
  • the method of the iterative solution is a limited-domain quasi-Newton method or a gradient descent method
  • the model adopted by the transform is a B-spline-based free deformation model
  • the interpolation processing method is bilinear interpolation. Method or B-spline interpolation.
  • the calculation formula of the similarity measure SSD of the local feature descriptors ZMLD A (i) and ZMLD B (i) in the step 3-1 is:
  • the transformation parameter T is a third-order B-spline function.
  • the number of iterations threshold 50 ⁇ ⁇ ⁇ 20.
  • a registration system for a non-rigid multi-mode medical image, the registration system comprising a first Zernike moment module, a first descriptive sub-module, a second a Zernike moment module, a second description submodule, and a registration module;
  • the first Zernike moment module is configured to obtain a 0th order 0th Zernike moment and a 1st order 1 Zernike moment of the reference image and output to the first description submodule, where the first description submodule is used to obtain local features of the reference image Descriptor and output to the registration module;
  • the second Zernike moment module is configured to obtain a 0th order 0 heavy Zernike moment of the floating image and a 1st order 1 Zernike moment and output to the second description submodule, the second descriptor
  • the module is configured to obtain a local feature descriptor of the floating image and output to the registration module;
  • the registration module is configured to obtain a registration image.
  • the registration module includes a solution unit for constructing an objective function according to a local feature descriptor of the reference image and the floating image, obtaining a transformation parameter and outputting, and the determination unit is configured to determine Whether the objective function conforms to the iterative stopping criterion, otherwise transforms the local feature descriptor of the floating image according to the transform parameter, performs interpolation processing on the transformed local feature descriptor, and performs local feature descriptor of the floating image after interpolation Updating the local feature descriptor of the original floating image; if yes, transforming the floating image according to the transformation parameter, and performing interpolation processing on the floating image to obtain the registration image.
  • the above technical solution proposed by the present invention constructs the local feature descriptor ZMLD by extracting the gray information and the edge feature information of the image by the 0th order zero Zernike moment and the 1st order 1 Zernike moment, respectively. It can effectively represent the local structural information of complex medical images under the condition of overcoming the noise of images and the rotation between image features, which provides an effective basis for the accurate evaluation of multimodal image similarity.
  • FIG. 1 is a schematic structural view of a registration system for a non-rigid multi-mode medical image of the present invention
  • Embodiment 2 is a flow chart of a registration method of the non-rigid multi-mode medical image of Embodiment 1;
  • Figure 3a is a reference image T2 used in Example 1 of the present invention and Comparative Example 2-4;
  • Figure 3b is a floating image T1 used in Example 1 of the present invention and Comparative Example 2-4;
  • Figure 3c is a floating image Gad used in Example 1 of the present invention and Comparative Example 2-4;
  • Figure 3d is a registration image T1-T2 obtained by the method of Embodiment 1 of the present invention.
  • Figure 3e is a registration image T1-T2 obtained by the method of Comparative Example 4.
  • Figure 3f is a registration image T1-T2 obtained by the method of Comparative Example 3 of the present invention.
  • Figure 3g is a registration image T1-T2 obtained by the method of Comparative Example 2;
  • Figure 3h is a registration image Gad-T2 obtained by the method of Embodiment 1 of the present invention.
  • Figure 3i is a registration image Gad-T2 obtained by the method of Comparative Example 4.
  • Figure 3j is a registration image Gad-T2 obtained by the method of Comparative Example 3;
  • Figure 3k is a registration image Gad-T2 obtained by the method of Comparative Example 2;
  • Figure 4a is a reference image T1 weighted MR used in Embodiment 1 of the present invention and Comparative Example 2-4;
  • Figure 4b is a floating image CT used in Example 1 of the present invention and Comparative Example 2-4;
  • Figure 4d is a registration image obtained by the method of Comparative Example 4 of the present invention.
  • Figure 4e is a registration image obtained by the method of Comparative Example 3 of the present invention.
  • Figure 4f is a registration image obtained by the method of Comparative Example 2 of the present invention.
  • the invention provides a registration system for a non-rigid multi-mode medical image, comprising a first Zernike moment module, a first description sub-module, a second Zernike moment module, a second description sub-module and a registration module, as shown in FIG.
  • the registration module further includes a solution unit and a determination unit;
  • An output end of the first Zernike moment module is connected to an input end of the first description sub-module, an output end of the first description sub-module is connected to a first input end of the solution unit, and an output of the second Zernike moment module
  • the end is connected to the input end of the second description sub-module, the output end of the second descriptive sub-module is connected to the second input end of the solution unit; the output end of the solution unit is connected
  • An input end of the determining unit, an output end of the determining unit is connected to an input end of the second Zernike moment module;
  • the first Zernike moment module is configured to obtain a 0th order 0th Zernike moment Z A 00 (i) of the reference image I A (i) and a 1st order 1 Zernike moment Z A 11 (i), and output the first
  • the description sub-module is used to obtain a local feature descriptor ZMLD A (i) of the reference image and is output
  • the second Zernike moment module is used to obtain a 0-order 0-weight Zernike moment Z B 00 of the floating image I B (i) (i And a 1st order 1 Zernike moment Z B 11 (i) for outputting and outputting a local feature descriptor ZMLD B (i) of the floating image for use according to the reference a local feature descriptor of the image and the floating image, constructing an objective function g(T ⁇ ) and obtaining a transformation parameter T ⁇ , the determination unit is configured to determine whether the objective function meets an iterative stop criterion, or otherwise transform
  • registration of the non-rigid multimodal medical image by the registration system includes the following steps:
  • N represents the side length of the first image block centered on the pixel point i, and is usually an odd number of 3 to 11.
  • the complexity, calculation efficiency, and registration precision of the image can be comprehensively selected; When the complexity is high, the registration accuracy is not required, or a high computational efficiency is required, a smaller value such as 3 or 5 may be taken, and a larger value may be taken; f(x, y) represents the pixel i as the origin.
  • the image functions of the first image block, ⁇ and ⁇ represent the polar and polar axes of the pixel point in the image function f(x, y), respectively, when the pixel point i is located at the edge of the reference image and the registration image, first The image function value in the image block that exceeds the position range of the original reference image and the floating image is filled with the image function value of the pixel adjacent thereto; x and y respectively represent any pixel in the image function f(x, y)
  • the reference image can be obtained by the above formula
  • the 0th order 0 heavy Zernike moment Z A 00 (i) and the 1st order 1 Zernike moment Z A 11 (i) at the same time, obtain the 0th order 0 Zernike of the pixel i at the same position in the floating image I B (i) Moment Z B 00 (i) and 1st order 1 Zernike moment Z B 11 (i);
  • Step 1-2 In the reference image, the similar distances D A 00 (i) and D A 11 (i) of the pixel point i and the remaining pixel points of the second image block centered thereon are obtained, and the calculation formula is as follows: Where
  • Step 1-3 obtains the local feature descriptor ZMLD, and the formula is as follows:
  • h 00 (i, r) and h 11 (i, r) are attenuation coefficients, which are calculated as follows:
  • ⁇ nm g (i) c 1 ⁇ 1.4826MED[
  • MED( ⁇ ) is the median operator
  • c 1 and c 2 are adjustment coefficients, and the value ranges from 0.5 ⁇ c 1 ⁇ 1, 0.5 ⁇ c 2 ⁇ 1;
  • the local feature descriptor ZMLD A (i) of the reference image can be obtained separately;
  • Step 2 is the same calculation method as step 1, and can be based on the 0th order 0 of the same position in the floating image I B (i), the Zernike moment Z B 00 (i) and the 1st order 1 Zernike moment Z B 11 (i) obtaining a local feature descriptor ZMLD B (i) of the floating image;
  • Step 3 constructs an objective function according to the local feature description sub-ZMLD A (i) and ZMLD B (i) of the reference image and the floating image, and finally completes the registration of the floating image and the reference image; here the transformation model is based on the B-spline
  • the Free-form Deformation (FFD) builds the objective function as an example to illustrate the registration process, including the following sub-steps:
  • ⁇ is a weight parameter, 0 ⁇ 1
  • T ⁇ is a third-order B-spline function associated with the coordinates (x, y) of the pixel point i, representing a transformation parameter of the floating image transformed into a registration image
  • the function g(T ⁇ ) is solved iterative
  • Step 3-3 If the number of iterations ⁇ is greater than or equal to the threshold number of iterations ⁇ (general 50 ⁇ ⁇ ⁇ 20), and g(T ⁇ ) ⁇ g (T ⁇ -1 ), transform according to the transformation parameter T ⁇ (x, y)
  • the floating image is subjected to interpolation processing by bilinear interpolation or B-spline interpolation on the transformed floating image to obtain the registration image I' B (i), thereby completing image registration; otherwise, returning Step 3-2.
  • This embodiment provides a registration method for a non-rigid multi-mode medical image, including the following steps, as shown in FIG. 1:
  • Step 1 calculates a reference image and a floating image, and according to formulas (1) and (2), obtains a 0th order 0 weight and a 1st order 1 weight Zernike moment of the first image block centered on each pixel point i, in this embodiment
  • Step 2 calculates a Zernike moments based local descriptor (ZMLD) based on the moment feature.
  • ZMLD Zernike moments based local descriptor
  • Step 2-1 calculates the similar distances D A 00 (i), D A 11 (i) of the pixel and the remaining pixels of the image block centered on the reference image and the floating image based on the image self-similarity.
  • D B 00 (i) and D B 11 (i) the formula is as follows:
  • Step 2-2 calculates the local feature descriptor ZMLD, and the formula is as follows:
  • ⁇ nm g (i) c 1 ⁇ 1.4826MED[
  • MED( ⁇ ) is a median operator
  • c 1 and c 2 are adjustment coefficients.
  • each reference image can be obtained local feature local feature descriptor ZMLD A (i) and the floating image descriptors ZMLD B (i);
  • M is the image size.
  • T represents a floating image.
  • Step 3-2 takes the minimum value of g(T) as the target, and uses the bounded quasi-Newton method (L-BFGS) to perform iterative solution.
  • L-BFGS bounded quasi-Newton method
  • Step 3-3 deforms the floating image by using the transformation parameter T obtained in step 3-2, and performs interpolation by bilinear interpolation to obtain a registration image corresponding to the floating image, and finally completes image registration.
  • the registration was carried out according to the ESSD method in (Med. Image Anal. 16(1) (2012) 1-17.).
  • the specific parameters are: selecting 7 ⁇ 7 image blocks, and calculating the entropy corresponding to the image block by using Gaussian weight, local normalization method and Parzen window estimation, thereby obtaining the ESSD corresponding to the entire image.
  • the registration was carried out according to the WLD method in (Sensors 13(6)(2013)7599-7613.).
  • T s represents the random deformation and is the gold standard of the evaluation
  • T c represents the deformation obtained by the registration algorithm
  • R represents the number of pixels used for image registration performance evaluation.
  • the simulated MR image was used for the registration accuracy test.
  • the simulated T1, T2 and PD weighted MR images used in Example 1 were taken from the BrainWeb database.
  • Table 1 lists the standard deviation and mean of the TRE obtained by each algorithm. It can be seen from Table 1 that when registering different MR images, the TRE provided in Example 1 has lower mean and standard deviation than other methods, which indicates that the method proposed by the present invention has the highest registration among all the comparative methods. Precision.
  • FIG. 3a is a reference image T2
  • FIG. 3b is a floating image T1
  • FIG. 3c is a floating image Gad
  • FIG. 3d is a registration image T1-T2 obtained by the method of Embodiment 1
  • FIG. 3e is a registration image T1 obtained by the method of Comparative Example 4.
  • -T2 is the registration image T1-T2 obtained by the method of Comparative Example 3
  • FIG. 3g is the registration image T1-T2 obtained by the method of Comparative Example 2
  • FIG. 3a is a reference image T2
  • FIG. 3c is a floating image Gad
  • FIG. 3d is a registration image T1-T2 obtained by the method of Embodiment 1
  • FIG. 3e is a registration image T1 obtained by the method of Comparative Example 4.
  • -T2 is the registration image T1-T2 obtained by the method of Comparative Example 3
  • FIG. 3g is the registration image T1-T2 obtained by the method of Comparative Example
  • FIG. 3h is the registration image Gad-T2 obtained by the method of Embodiment 1.
  • Figure 3i is the registration image Gad-T2 obtained by the method of Comparative Example 4, and
  • Figure 3j is the matching method obtained by the method of Comparative Example 3.
  • the quasi image Gad-T2 Fig. 3k is the registration image Gad-T2 obtained by the method of Comparative Example 2.
  • the contrast ratio 2-4 cannot effectively correct the deformation of the lowermost portion of the contour.
  • Proportion 3-4 does not provide a good correction for the deformation of the right portion of the box, while Embodiment 1 can effectively correct the above deformation. It can be seen that the registration image obtained in Example 1 is more similar to the reference image than Comparative Example 2-4.
  • the above visual comparison proves that the local feature descriptor of the present invention has rotation invariance, and the representation of the image is more accurate, so that the registration of the non-rigid multi-mode medical image is superior.
  • Table 2 shows the corresponding TRE results for each method.
  • the results in Table 2 indicate that the registration method proposed by the present invention provides a lower TRE and thus a higher registration accuracy than other methods.
  • FIG. 4 is a result of registration of real CT and MR images by the methods of Example 1 and Comparative Examples 2-4.
  • 4a is a reference image T1 weighted MR
  • FIG. 4b is a floating image CT
  • FIG. 4c is a registration image obtained by the method of Embodiment 1
  • FIG. 4d is a registration image obtained by the method of Comparative Example 4
  • FIG. 4e is a comparison example 3
  • FIG. 4f is the registration image obtained by the method of Comparative Example 2. It can be seen that the method of Comparative Example 3-4 is difficult to effectively correct the inside of the box.
  • Table 3 shows the average of the TREs for each set of real CT images and MR images, which are the real CT images and MR images used in Figure 4.
  • the method of the present invention can achieve a lower TRE by registering the CT-MR images of Groups 1 to 5 compared to other registration methods, which illustrates the method proposed by the present invention.
  • a more proportional algorithm for CT-MR image registration has higher registration accuracy.

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Medical Informatics (AREA)
  • Public Health (AREA)
  • General Health & Medical Sciences (AREA)
  • Primary Health Care (AREA)
  • Epidemiology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Biomedical Technology (AREA)
  • Pathology (AREA)
  • Databases & Information Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Quality & Reliability (AREA)
  • Image Processing (AREA)
  • Algebra (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)

Abstract

本发明公开了一种非刚性多模医学图像的配准方法及系统。该配准方法包括：根据参考图像0阶0重Zernike矩以及1阶1重Zernike矩，获得参考图像的局部特征描述子；根据浮动图像的0阶0重Zernike矩以及1阶1重Zernike矩，获得浮动图像的局部特征描述子；最后根据所述参考图像以及浮动图像的局部特征描述子，获得配准图像。本发明利用多模医学图像的自相似性，采用基于Zernike矩的局部特征描述子，从而将非刚性多模医学图像的配准问题转化为单模医学图像配准问题，大大提高了非刚性多模医学图像配准的精度与鲁棒性。

Description

一种非刚性多模医学图像的配准方法及系统 【技术领域】
本发明属于图像处理与分析中的图像配准领域，更具体地，涉及一种非刚性多模医学图像的配准方法及系统。
【背景技术】
医学影像技术是现代医学中的重要组成部分，对疾病的诊断和治疗有着革命性的意义。现有的成像技术虽然发展迅速，但由于其成像原理的限制，往往只能提供单一且有限的信息，比如超声、CT和MRI可以提供器官的解剖信息，但是不能提供他们的功能信息；PET可以较好的提供代谢信息，但不能清晰的提供器官的形态信息。为此，医生往往需要融合不同模态图像的信息来提高医学影像诊断的准确性。而要进行多模图像融合，多模图像配准必不可少。
图像配准的目标是寻找参考图像和浮动图像之间的对应关系，大体上图像配准包含三个主要部分：形变模型、相似性度量和优化方法。而非刚性多模医学图像配准由于灰度的扭曲、病人的呼吸和体位的移动等因素的影响，一直是图像配准领域的一个难题。
目前解决非刚性多模图像配准问题的方法大体上被分为两类：第一类是基于互信息测度的配准方法，然而这类方法通常都没有考虑图像的局部特征结构，计算耗时，而且容易陷入局部极值导致配准结果不准确。第二类方法通过图像结构表征方法将多模图像配准简化为单模图像配准，如通过熵图、韦伯局部特征描述子(Weber Local Descriptor，WLD)及基于模态独立邻域描述子(Modality Independent Neighborhood Descriptor，MIND)等特征对图像结构进行表征，然后利用表征结果的差值平方和(Sum of Squared Difference，SSD)作为配准测度来实现图像配准。其中，基于熵图的方法是 通过估计图像块的灰度概率密度函数计算得到每个图像块的熵值，由此获得整幅图像的熵图，而基于WLD的方法则是利用图像的拉普拉斯运算来描述其局部结构特征。这两种方法能有效克服多模图像间灰度差异的不利影响，但是基于熵图和WLD的特征都对图像中的噪声敏感，因而在图像中存在噪声的情况下，难以产生精确的配准结果。基于MIND的配准方法通过图像块间的欧式距离来评价图像自相似性，利用不同图像块间的相似性之和来刻画图像局部结构特征，但该方法在评价图像自相似性时仅考虑了图像块之间的平移不变性，忽略了图像块间可能存在的旋转特性，因此在图像块间存在旋转的情况下，该方法难以提供精确的配准结果。
【发明内容】
针对现有技术的上述缺陷或改进需求，本发明提供了一种基于图像自相似性的非刚性多模医学图像精确配准方法及系统，从而通过图像的自相似性，采用Zernike矩来构建图像描述特征，从而实现非刚性多模医学图像的配准。
为实现上述目的，按照本发明的一方面，提供了一种非刚性多模医学图像的配准方法，包括：
根据参考图像以及浮动图像的0阶0重Zernike矩以及1阶1重Zernike矩，分别获得参考图像以及浮动图像的局部特征描述子，并获得配准图像。优选地，所述配准方法具体包括以下步骤：
步骤1根据参考图像IA(i)中像素点i的0阶0重Zernike矩ZA 00(i)以及1阶1重Zernike矩ZA 11(i)，获得参考图像的局部特征描述子ZMLDA(i)；其中，i为1～M的整数，M为参考图像以及浮动图像的尺寸；
步骤2根据浮动图像IB(i)中相同位置的像素点i的0阶0重Zernike矩ZB 00(i)以及1阶1重Zernike矩ZB 11(i)，获得浮动图像的局部特征描述子ZMLDB(i)；
步骤3根据参考图像和浮动图像的局部特征描述子ZMLDA(i)以及 ZMLDB(i)，建立目标函数g(T)；根据目标函数，获得变换参数，根据变换参数变换所述浮动图像，并对变换后的浮动图像进行插值处理，获得所述配准图像。
作为进一步优选地，所述步骤2包括以下子步骤：
步骤2-1在参考图像中，获得像素点i与以其为中心的图像块内的其余像素点对应0阶0重Zernike矩的相似距离DA 00(i)，以及对应1阶1重Zernike矩的相似距离DA 11(i)；同时，在浮动图像中，获得像素点i与以其为中心的图像块内的其余像素点对应0阶0重Zernike矩的相似距离DB 00(i)，以及对应1阶1重Zernike矩的相似距离DB 11(i)；
步骤2-2获得参考图像的局部特征描述子ZMLDA(i)以及浮动图像的局部特征描述子ZMLDB(i)；
Figure PCTCN2016101547-appb-000001
其中，h00 A(i)、h11 A(i)、h00 B(i)以及h00 B(i)为衰减系数。
作为更进一步优选地，所述步骤2-1中图像块的边长为3，在所述步骤2-2中，h00 A(i)＝{σ00 A(i)+c1·1.4826MED[|σ00 A(i)-MED(|σ00 A(i)|)|]}2，
Figure PCTCN2016101547-appb-000002
h11 A(i)＝{σ11 A(i)+c1·1.4826MED[|σ11 A(i)-MED(|σ11 A(i)|)|]}2，
Figure PCTCN2016101547-appb-000003
h00 B(i)＝{σ00 B(i)+c1·1.4826MED[|σ00 B(i)-MED(|σ00 B(i)|)|]}2，
Figure PCTCN2016101547-appb-000004
h11 B(i)＝{σ11 B(i)+c1·1.4826MED[|σ11 B(i)-MED(|σ11 B(i)|)|]}2，
Figure PCTCN2016101547-appb-000005
其中，MED(·)为中值运算符，c1和c2为0.5～1的常数，j为以像素点i为中心的图像块的任一像素点。
作为进一步优选地，所述步骤3具体包括以下子步骤：
步骤3-1建立目标函数g(Tτ)＝SSD+αR(Tτ)；其中，SSD表示局部特征描述子ZMLDA(i)以及ZMLDB(i)的相似性测度，α为常数，0<α<1，R(Tτ)表示正则化项，迭代次数τ＝1；对所述目标函数g(Tτ)迭代求解，获得变换参数Tτ；
步骤3-2根据变换参数Tτ变换所述浮动图像的局部特征描述子ZMLDB(i)，对变换后的局部特征描述子进行插值处理，并以插值处理后的局部特征描述子更新原局部特征描述子ZMLDB(i)，τ＝τ+1；对所述目标函数g(Tτ)迭代求解，获得变换参数Tτ；
步骤3-3如果迭代次数τ大于等于迭代次数阈值ξ，且g(Tτ)≥g(Tτ-1)，则根据变换参数Tτ变换所述浮动图像，并对变换后的浮动图像进行插值处理，获得所述配准图像；否则返回步骤3-2。
作为更进一步优选地，所述迭代求解的方法为限域拟牛顿法或梯度下降法，所述变换采用的模型为基于B样条的自由形变模型，所述插值处理的方法为双线性插值法或B样条插值法。
作为更进一步优选地，所述步骤3-1中的局部特征描述子ZMLDA(i)以及ZMLDB(i)的相似性测度SSD的计算公式为：
Figure PCTCN2016101547-appb-000006
作为更进一步优选地，所述变换参数T为三阶B样条函数。
作为更进一步优选地，所述迭代次数阈值50≥ξ≥20。
按照本发明的另一方面，还提供了一种非刚性多模医学图像的配准系统，所述配准系统包括第一Zernike矩模块、第一描述子模块、第二 Zernike矩模块、第二描述子模块以及配准模块；
所述第一Zernike矩模块用于获得参考图像的0阶0重Zernike矩以及1阶1重Zernike矩并输出至第一描述子模块，所述第一描述子模块用于获得参考图像的局部特征描述子并输出至配准模块；所述第二Zernike矩模块用于获得浮动图像的0阶0重Zernike矩以及1阶1重Zernike矩并输出至第二描述子模块，所述第二描述子模块用于获得浮动图像的局部特征描述子并输出至配准模块；所述配准模块用于获得配准图像。
优选地，所述配准模块包括求解单元以及判定单元，所述求解单元用于根据参考图像和浮动图像的局部特征描述子，构建目标函数，获得变换参数并输出，所述判定单元用于判断所述目标函数是否符合迭代停止标准，否则根据变换参数变换所述浮动图像的局部特征描述子，对变换后的局部特征描述子进行插值处理，并以插值处理后的浮动图像的局部特征描述子更新原浮动图像的局部特征描述子；是则根据变换参数变换所述浮动图像，对浮动图像进行插值处理，获得所述配准图像。
本发明提出的上述技术方案与现有技术相比，由于通过0阶0重Zernike矩以及1阶1重Zernike矩分别提取图像的灰度信息和边缘特征信息，由此构建局部特征描述子ZMLD，可在克服图像存在噪声以及图像特征间存在旋转的情况下，有效表征复杂医学图像的局部结构信息，为多模图像相似性的准确评价提供了有效依据。
【附图说明】
图1为本发明非刚性多模医学图像的配准系统的结构示意图；
图2为实施例1非刚性多模医学图像的配准方法的流程图；
图3a为本发明实施例1以及对比例2-4所用的参考图像T2；
图3b为本发明实施例1以及对比例2-4所用的浮动图像T1；
图3c为本发明实施例1以及对比例2-4所用的浮动图像Gad；
图3d为本发明实施例1方法获得的配准图像T1-T2；
图3e为本发明对比例4方法获得的配准图像T1-T2；
图3f为本发明对比例3方法获得的配准图像T1-T2；
图3g为本发明对比例2方法获得的配准图像T1-T2；
图3h为本发明实施例1方法获得的配准图像Gad-T2；
图3i为本发明对比例4方法获得的配准图像Gad-T2；
图3j为本发明对比例3方法获得的配准图像Gad-T2；
图3k为本发明对比例2方法获得的配准图像Gad-T2；
图4a为本发明实施例1以及对比例2-4所用的参考图像T1加权MR；
图4b为本发明实施例1以及对比例2-4所用的浮动图像CT；
图4c为本发明实施例1方法获得的配准图像；
图4d为本发明对比例4方法获得的配准图像；
图4e为本发明对比例3方法获得的配准图像；
图4f为本发明对比例2方法获得的配准图像。
【具体实施方式】
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
本发明提供了一种非刚性多模医学图像的配准系统，包括第一Zernike矩模块、第一描述子模块、第二Zernike矩模块、第二描述子模块以及配准模块，如图1所示，其中，配准模块又包括求解单元以及判定单元；
所述第一Zernike矩模块的输出端连接第一描述子模块的输入端，所述第一描述子模块的输出端连接所述求解单元的第一输入端，所述第二Zernike矩模块的输出端连接第二描述子模块的输入端，所述第二描述子模块的输出端连接所述求解单元的第二输入端；所述求解单元的输出端连接 所述判定单元的输入端，所述判定单元的输出端连接所述第二Zernike矩模块的输入端；
所述第一Zernike矩模块用于获得参考图像IA(i)的0阶0重Zernike矩ZA 00(i)以及1阶1重Zernike矩ZA 11(i)并输出，所述第一描述子模块用于获得参考图像的局部特征描述子ZMLDA(i)并输出；所述第二Zernike矩模块用于获得浮动图像IB(i)的0阶0重Zernike矩ZB 00(i)以及1阶1重Zernike矩ZB 11(i)并输出，所述第二描述子模块用于获得浮动图像的局部特征描述子ZMLDB(i)并输出，所述求解单元用于根据参考图像和浮动图像的局部特征描述子，构建目标函数g(Tτ)并获得变换参数Tτ，所述判定单元用于判断所述目标函数是否符合迭代停止标准，否则根据变换参数变换所述浮动图像的局部特征描述子，对变换后的局部特征描述子进行插值处理，并以插值处理后的浮动图像的局部特征描述子更新原浮动图像的局部特征描述子；是则根据变换参数变换所述浮动图像，对浮动图像进行插值处理，获得所述配准图像。
具体地，该配准系统进行非刚性多模医学图像的配准包括以下步骤：
步骤1根据参考图像IA(i)中像素点i的0阶0重Zernike矩ZA 00(i)以及1阶1重Zernike矩ZA 11(i)，获得参考图像的局部特征描述子ZMLDA(i)；其中，i为1～M的整数，M为参考图像以及浮动图像的尺寸，即当图像的长为X，宽为Y时，M＝X×Y；
步骤1-1Zernike矩的计算公式如下：
Figure PCTCN2016101547-appb-000007
Figure PCTCN2016101547-appb-000008
Figure PCTCN2016101547-appb-000009
其中，N表示以像素点i为中心的第一图像块的边长，通常为3～11的奇数，在应用中，可综合考虑图像的复杂程度、计算效率和配准精度来选择；当图像的复杂程度较高、配准精度要求不大或者需要较高的计算效率时，可取3、5等较小值，反之可取较大值；f(x,y)表示以像素点i为原点的第一图像块的图像函数，θ和ρ分别代表像素点在图像函数f(x,y)中的极角和极轴，当像素点i位于参考图像以及配准图像的边缘处时，第一图像块中超出原参考图像以及浮动图像的位置范围的图像函数值则以与其相邻的像素点的图像函数值进行填充；x、y分别表示图像函数f(x,y)中任意一像素点在第一图像块中与第一图像块的中心点i的相对坐标；
Figure PCTCN2016101547-appb-000010
λN为归一化系数，λN＝N2；m表示Zernike矩的重数，n表示Zernike矩的阶数，s＝0－(n－|m|)/2，在本发明中，由于n＝m＝1或n＝m＝0，因此s＝0，R00(ρ)＝1，R11(ρ)＝ρ；
则当n＝m＝0时，则获得了图像的0阶0重Zernike矩，当n＝m＝1时，则获得了图像的1阶1重Zernike矩；因此，由以上公式可以获得参考图像的0阶0重Zernike矩ZA 00(i)以及1阶1重Zernike矩ZA 11(i)，同时，获得浮动图像IB(i)中相同位置的像素点i的0阶0重Zernike矩ZB 00(i)以及1阶1重Zernike矩ZB 11(i)；
步骤1-2在参考图像中，获得像素点i与以其为中心的第二图像块的其余像素点的相似距离DA 00(i)和DA 11(i)，其计算公式如下：
Figure PCTCN2016101547-appb-000011
其中，|r|表示第二图像块的边长，通常为大于等于3的奇数，j为第二图像块中不为像素点i的任意一像素点，|Znm(i)|和|Znm(j)|是分别为像素点i和像素点j的n阶m重Zernike矩的模，当像素点i位于参考图像的边缘处时，第二图像块中超出参考图像位置范围的Zernike矩的模则以与其相邻的像素点的Zernike矩的模进行填充；由此，可获得DA 00(i)、DA 11(i)、DB 00(i)以及DB 11(i)；
步骤1-3获得局部特征描述子ZMLD，公式如下：
Figure PCTCN2016101547-appb-000012
其中，h00(i,r)和h11(i,r)为衰减系数，其计算公式如下：
hnm(i)＝[σnm l(i)+σnm g(i)]2,n＝m＝0 or n＝m＝1
例如，当|r|＝3时，σnm l(i)和σnm g(i)分别为：
σnm g(i)＝c1·1.4826MED[|σnm l(i)-MED(|σnm l(i)|)|]，
Figure PCTCN2016101547-appb-000013
其中，MED(·)为中值运算符，c1和c2为调节系数，其取值范围为0.5＜c1≤1，0.5＜c2≤1；
根据上述公式，可分别获得参考图像的局部特征描述子ZMLDA(i)；
步骤2以和步骤1相同的计算方法，可根据浮动图像IB(i)中相同位置的像素点i的0阶0重Zernike矩ZB 00(i)以及1阶1重Zernike矩ZB 11(i)，获得浮动图像的局部特征描述子ZMLDB(i)；
步骤3根据参考图像和浮动图像的局部特征描述子ZMLDA(i)以及ZMLDB(i)，构建目标函数，并最终完成浮动图像与参考图像的配准；这里以变换模型为基于B样条自由形变模型(Free-form Deformation,FFD)构建目标函数为例，说明该配准过程，具体包括以下子步骤：
步骤3-1建立目标函数g(Tτ)＝SSD+αR(Tτ)，g(Tτ)的值越小，则浮动图像与参考图像的相似性越大，SSD表示局部特征描述子ZMLDA(i)以及ZMLDB(i)的相似性测度，
Figure PCTCN2016101547-appb-000014
其中，α为权重参数，0<α<1，Tτ为与像素点i的坐标(x，y)相关的三阶B样条函数，表示浮动图像变换为配准图像的变换参数；R(Tτ)为正则化项，其计算公式为：
Figure PCTCN2016101547-appb-000015
其中，x和y分别表示像素点i在浮动图像中的横坐标和纵坐标，X和Y分别表示浮动图像的长和宽，则X×Y＝M；迭代次数τ＝1；对所述目标函数g(Tτ)用限域拟牛顿法或梯度下降法迭代求解，获得变换参数Tτ(x,y)；
步骤3-2根据变换参数Tτ(x,y)变换所述浮动图像的局部特征描述子ZMLDB(i)，对变换后的局部特征描述子通过双线性插值法或B样条插值法插值处理，并以插值处理后的局部特征描述子更新原局部特征描述子ZMLDB(i)，τ＝τ+1；对所述目标函数g(Tτ)迭代求解，获得变换参数Tτ；
步骤3-3如果迭代次数τ大于等于迭代次数阈值ξ(一般50≥ξ≥20)，且g(Tτ)≥g(Tτ-1)，则根据变换参数Tτ(x,y)变换所述浮动图像，并对变换后的浮动图像通过双线性插值法或B样条插值法进行插值处理，获得所述配准图像I’B(i)，由此完成图像配准；否则返回步骤3-2。
实施例1
本实施例提供了一种非刚性多模医学图像的配准方法，包括以下步骤，如图1所示：
步骤1计算参考图像和浮动图像中，根据公式(1)和(2)，获得以每个像素点i为中心的第一图像块的0阶0重和1阶1重Zernike矩，在本实施例中取N＝5；最终得到参考图像IA的矩特征ZA 00(i)以及ZA 11(i)，以及浮动图像IB的矩特征ZB 00(i)以及ZB 11(i)；
步骤2根据矩特征计算得到一种基于Zernike矩的局部特征描述子(Zernike moments based local descriptor，简称ZMLD)。
步骤2-1基于图像自相似性，计算参考图像和浮动图像中获得像素点i与以其为中心的图像块的其余像素点的相似距离DA 00(i)、DA 11(i)、DB 00(i)以及DB 11(i)，公式如下：
Figure PCTCN2016101547-appb-000016
其中|Znm(i)|和|Znm(j)|是像素点i和j处的n阶m重Zernike矩的模；|r|表示第二图像块的边长，|r|＝3以保证较高的配准精度；在本实施例中，当|r|＞3时，配准精度没有提高而增加了时间损耗，|r|＜3时，明显配准精度受到影响；j为第二图像块中不为像素点i的任意一像素点；
步骤2-2计算局部特征描述子ZMLD，公式如下：
Figure PCTCN2016101547-appb-000017
其中，h00(i,r)和h11(i,r)为衰减系数，其计算公式为
hnm(i)＝[σnm l(i)+σnm g(i)]2,n＝m＝0 or n＝m＝1
由于在本实施例中，|r|＝3，σnm l(i)和σnm g(i)分别为：
σnm g(i)＝c1·1.4826MED[|σnm l(i)-MED(|σnm l(i)|)|]，
Figure PCTCN2016101547-appb-000018
其中，MED(·)为中值运算符，c1和c2为调节系数，在本实施例中，c1＝c2＝0.8。
根据上述公式，可分别获得参考图像的局部特征描述子ZMLDA(i)以及浮动图像的局部特征描述子ZMLDB(i)；
步骤3-1采用基于B样条自由形变模型(Free-form Deformation,FFD)作为变换模型，建立目标函数g(T)＝SSD+αR(T)；其中，SSD表示局部特征描述子ZMLDA(i)以及ZMLDB(i)的相似性测度，
Figure PCTCN2016101547-appb-000019
M为图像尺寸，在本实施例中由于参考图像和浮动图像的长X和宽Y都为256，因此M＝X×Y＝256×256；α为权重参数，α＝0.015，T表示浮动图像变换为配准图像的变换参数；R(T)为正则化项，其计算公式为：
Figure PCTCN2016101547-appb-000020
步骤3-2以g(T)的值最小为目标，使用限域拟牛顿法(L-BFGS)进行迭代求解，当g(T)获得最小值，且迭代次数大于等于迭代次数阈值ξ时，迭代停止，获得变换参数T；
步骤3-3通过步骤3-2中获得的变换参数T，对浮动图像变形，并通过双线性插值法进行插值，获得浮动图像对应的配准图像，并最终完成图像配准。
对比例1
按照(Pattern Recognit.32(1)(1999)71-86.)里的NMI方法实现配准。
对比例2
按照(Med.Image Anal.16(1)(2012)1-17.)里的ESSD方法实现配准。其中，具体参数为：选择7×7的图像块，并利用高斯权重、局部归一化方法和Parzen窗估计来计算图像块对应的熵，由此获得整个图像对应的ESSD。
对比例3
按照(Sensors 13(6)(2013)7599-7613.)里的WLD方法实现配准。具体参数为：WLD的计算选择半径R＝1和R＝2，构建相似性测度的块大小为7×7，权重项γ＝0.01。
对比例4
按照(Med.Image Anal.16(7)(2012)1423-1435.)里的MIND方法实现配准。具体参数为：图像块大小选择3×3。
实验结果分析
为了进一步体现本发明的优点，我们将实施例1与对比例1-4的配准精度进行比较。配准精度采用目标配准误差TRE进行评价，这里TRE定 义为：
Figure PCTCN2016101547-appb-000021
其中Ts表示随机形变，也是评价的金标准，Tc表示由配准算法得到的形变，R表示用于图像配准性能评价的像素个数。
采用仿真MR图像进行配准精度测试，实施例1所用的仿真T1、T2和PD加权MR图像均取自BrainWeb数据库，表1列出了各算法所得TRE的标准差和均值。从表1可看出，对不同MR图像进行配准时，实施例1提供的TRE其均值和标准差皆低于其它方法，这说明本发明提出的方法在所有比较的方法中具有最高的配准精度。
表1 各方法在T1-T2，PD-T2和T1-PD图像配准时的TRE(mm)对比
Figure PCTCN2016101547-appb-000022
为更直观地显示本发明相对于其余方法的优越性，我们提供了实施例1与对比例2-4对应配准图像的视觉效果图，如图3所示。图3a为参考图像T2，图3b为浮动图像T1，图3c为浮动图像Gad，图3d为实施例1方法获得的配准图像T1-T2，图3e为对比例4方法获得的配准图像T1-T2，图3f为对比例3方法获得的配准图像T1-T2，图3g为对比例2方法获得的配准图像T1-T2，图3h为实施例1方法获得的配准图像Gad-T2，图3i为对比例4方法获得的配准图像Gad-T2，图3j为对比例3方法获得的配 准图像Gad-T2，图3k为对比例2方法获得的配准图像Gad-T2。
以方框内所示的图像的部分轮廓为例，对配准图像T1-T2而言，对比例2-4无法有效更正轮廓最下方部分的变形，对配准图像Gad-T2而言，对比例3-4无法对方框内右侧部分的形变取得好的更正效果，而实施例1却能对上述变形进行有效更正。可以看出，与对比例2-4相比，实施例1得到的配准图像与参考图像相似度更高。上述视觉比较充分证明了本发明的局部特征描述子具有旋转不变性，对于图像的表征更加准确，从而在非刚性多模医学图像的配准上具有优越性。
采用Altas数据库中的T1、T2和Grad加权MR图像评价所列方法的配准精度，表2给出了各方法对应TRE结果。表2中的结果表明：与其它方法相比，本发明提出的配准方法可提供更低的TRE，因而具有更高的配准精度。
表2 各方法在T1-T2、Gad-T2、Gad-T1图像配准时的TRE(mm)对比
Figure PCTCN2016101547-appb-000023
针对CT与MR图像进行了配准精度的比较，我们对浮动CT图像进行了五种随机的变形处理。图4是实施例1以及对比例2-4的方法对真实CT和MR图像进行配准的结果。其中，图4a为参考图像T1加权MR，图4b为浮动图像CT，图4c为实施例1方法获得的配准图像，图4d为对比例4方法获得的配准图像，图4e为对比例3方法获得的配准图像，图4f为对比例2方法获得的配准图像。可看出，对比例3-4的方法难以有效更正方框内 轮廓的变形，对比例2提供的配准图像中最外部轮廓局部出现失真，而实施例1方法可获得良好的配准图像。表3显示了每组真实CT图像与MR图像配准的TRE的平均值，其中，组5图像即为图4中所用的真实CT图像与MR图像。
表3 各方法在CT-MR图像配准时的TRE(mm)对比
Figure PCTCN2016101547-appb-000024
从表3中我们可以看出，与其它配准方法相比，本发明的方法对组1～组5的CT-MR图像进行配准皆可取得更低的TRE，这说明本发明提出的方法在CT-MR图像配准方面较对比例的算法具有更高的配准精度。
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (9)

  1. 一种非刚性多模医学图像的配准方法，其特征在于，包括：
    根据参考图像以及浮动图像的0阶0重Zernike矩以及1阶1重Zernike矩，分别获得参考图像以及浮动图像的局部特征描述子，并获得配准图像。
  2. 如权利要求1所述的配准方法，其特征在于，具体包括以下步骤：
    步骤1根据参考图像IA(i)中像素点i的0阶0重Zernike矩ZA 00(i)以及1阶1重Zernike矩ZA 11(i)，获得参考图像的局部特征描述子ZMLDA(i)；其中，i为1～M的整数，M为参考图像以及浮动图像的尺寸；
    步骤2根据浮动图像IB(i)中相同位置的像素点i的0阶0重Zernike矩ZB 00(i)以及1阶1重Zernike矩ZB 11(i)，获得浮动图像的局部特征描述子ZMLDB(i)；
    步骤3根据参考图像和浮动图像的局部特征描述子ZMLDA(i)以及ZMLDB(i)，建立目标函数g(T)；根据目标函数，获得变换参数，根据变换参数变换所述浮动图像，并对变换后的浮动图像进行插值处理，获得所述配准图像。
  3. 如权利要求2所述的配准方法，其特征在于，所述步骤2包括以下子步骤：
    步骤2-1在参考图像中，获得像素点i与以其为中心的图像块内的其余像素点对应0阶0重Zernike矩的相似距离DA 00(i)，以及对应1阶1重Zernike矩的相似距离DA 11(i)；同时，在浮动图像中，获得像素点i与以其为中心的图像块内的其余像素点对应0阶0重Zernike矩的相似距离DB 00(i)，以及对应1阶1重Zernike矩的相似距离DB 11(i)；
    步骤2-2获得参考图像的局部特征描述子ZMLDA(i)以及浮动图像的局部特征描述子ZMLDB(i)；
    Figure PCTCN2016101547-appb-100001
    其中，h00 A(i)、h11 A(i)、h00 B(i)以及h00 B(i)为衰减系数。
  4. 如权利要求3所述的配准方法，其特征在于，所述步骤2-1中图像块的边长为3，在所述步骤2-2中，
    h00 A(i)＝{σ00 A(i)+c1·1.4826MED[|σ00 A(i)-MED(|σ00 A(i)|)|]}2，
    Figure PCTCN2016101547-appb-100002
    h11 A(i)＝{σ11 A(i)+c1·1.4826MED[|σ11 A(i)-MED(|σ11 A(i)|)|]}2，
    Figure PCTCN2016101547-appb-100003
    h00 B(i)＝{σ00 B(i)+c1·1.4826MED[|σ00 B(i)-MED(|σ00 B(i)|)|]}2，
    Figure PCTCN2016101547-appb-100004
    h11 B(i)＝{σ11 B(i)+c1·1.4826MED[|σ11 B(i)-MED(|σ11 B(i)|)|]}2，
    Figure PCTCN2016101547-appb-100005
    其中，MED(·)为中值运算符，c1和c2为0.5～1的常数，j为以像素点i为中心的图像块内的任一像素点。
  5. 如权利要求2所述的配准方法，其特征在于，所述步骤3具体包括以下子步骤：
    步骤3-1建立目标函数g(Tτ)＝SSD+αR(Tτ)；其中，SSD表示局部特征描述子ZMLDA(i)以及ZMLDB(i)的相似性测度，α为常数，0<α<1，R(Tτ)表示正则化项，迭代次数τ＝1；对所述目标函数g(Tτ)迭代求解，获得变换参数Tτ；
    步骤3-2根据变换参数Tτ变换所述浮动图像的局部特征描述子ZMLDB(i)，对变换后的局部特征描述子进行插值处理，并以插值处理后的局部特征描述子更新原局部特征描述子ZMLDB(i)，τ＝τ+1；对所述目标函数g(Tτ)迭代求解，获得变换参数Tτ；
    步骤3-3如果迭代次数τ大于等于迭代次数阈值ξ，且g(Tτ)≥g(Tτ-1)，则根据变换参数Tτ变换所述浮动图像，并对变换后的浮动图像进行插值处理，获得所述配准图像；否则返回步骤3-2。
  6. 如权利要求5所述的配准方法，其特征在于，所述步骤3-1中的相似性测度SSD的计算公式为：
    Figure PCTCN2016101547-appb-100006
  7. 如权利要求5所述的配准方法，其特征在于，所述迭代求解的方法为限域拟牛顿法或梯度下降法，所述变换采用的模型为基于B样条的自由形变模型，所述插值处理的方法为双线性插值法或B样条插值法。
  8. 基于权利要求1-7中任意一项所述的配准方法的配准系统，其特征在于，所述配准系统包括第一Zernike矩模块、第一描述子模块、第二Zernike矩模块、第二描述子模块以及配准模块；
    所述第一Zernike矩模块用于获得参考图像的0阶0重Zernike矩以及1阶1重Zernike矩并输出至第一描述子模块，所述第一描述子模块用于获得参考图像的局部特征描述子并输出至配准模块；所述第二Zernike矩模块用于获得浮动图像的0阶0重Zernike矩以及1阶1重Zernike矩并输出至第二描述子模块，所述第二描述子模块用于获得浮动图像的局部特征描述子并输出至配准模块；所述配准模块用于获得配准图像。
  9. 如权利要求8所述的配准系统，其特征在于，所述配准模块包括求解单元以及判定单元，所述求解单元用于根据参考图像和浮动图像的局部 特征描述子，构建目标函数，获得变换参数并输出，所述判定单元用于判断所述目标函数是否符合迭代停止标准，否则根据变换参数变换所述浮动图像的局部特征描述子，对变换后的局部特征描述子进行插值处理，并以插值处理后的浮动图像的局部特征描述子更新原浮动图像的局部特征描述子；是则根据变换参数变换所述浮动图像，对浮动图像进行插值处理，获得所述配准图像。
PCT/CN2016/101547 2016-06-30 2016-10-09 一种非刚性多模医学图像的配准方法及系统 Ceased WO2018000652A1 (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US16/094,473 US10853941B2 (en) 2016-06-30 2016-10-09 Registration method and system for non-rigid multi-modal medical image

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN201610506619.9 2016-06-30
CN201610506619.9A CN106204550B (zh) 2016-06-30 2016-06-30 一种非刚性多模医学图像的配准方法及系统

Publications (1)

Publication Number Publication Date
WO2018000652A1 true WO2018000652A1 (zh) 2018-01-04

Family

ID=57463835

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2016/101547 Ceased WO2018000652A1 (zh) 2016-06-30 2016-10-09 一种非刚性多模医学图像的配准方法及系统

Country Status (3)

Country Link
US (1) US10853941B2 (zh)
CN (1) CN106204550B (zh)
WO (1) WO2018000652A1 (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111754554A (zh) * 2020-06-28 2020-10-09 上海应用技术大学 颅脑多模态医学图像配准方法
CN114581497A (zh) * 2022-04-16 2022-06-03 合肥学院 一种用于大气湍流图像畸变校正的改进b样条非刚性配准方法
CN115496785A (zh) * 2022-08-17 2022-12-20 中科超精(南京)科技有限公司 一种自适应约束模型的医学影像配准方法
CN115546268A (zh) * 2022-09-23 2022-12-30 中国人民解放军国防科技大学 多模态遥感图像配准方法、系统、终端设备及存储介质

Families Citing this family (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10398382B2 (en) * 2016-11-03 2019-09-03 Siemens Medical Solutions Usa, Inc. Respiratory motion estimation in projection domain in nuclear medical imaging
CN108416802B (zh) * 2018-03-05 2020-09-18 华中科技大学 一种基于深度学习的多模医学图像非刚性配准方法及系统
CN108711168A (zh) * 2018-06-04 2018-10-26 中北大学 基于zmld与gc离散优化的非刚性多模态医学图像配准方法
CN110517299B (zh) * 2019-07-15 2021-10-26 温州医科大学附属眼视光医院 基于局部特征熵的弹性图像配准算法
CN110533641A (zh) * 2019-08-20 2019-12-03 东软医疗系统股份有限公司 一种多模态医学图像配准方法和装置
CN110599527A (zh) * 2019-08-23 2019-12-20 首都医科大学宣武医院 一种mra影像数据的配准方法及装置
US20230059132A1 (en) * 2019-09-26 2023-02-23 The General Hospital Corporation System and method for deep learning for inverse problems without training data
CN111798500B (zh) * 2020-07-20 2023-06-23 陕西科技大学 一种基于层次邻域谱特征的微分同胚非刚性配准算法
CN112017221B (zh) * 2020-08-27 2022-12-20 北京理工大学 基于尺度空间的多模态图像配准方法、装置和设备
CN113450294A (zh) * 2021-06-07 2021-09-28 刘星宇 多模态医学图像配准融合方法、装置及电子设备
CN114022521B (zh) * 2021-10-13 2024-09-13 华中科技大学 一种非刚性多模医学图像的配准方法及系统
CN114255265B (zh) * 2021-12-01 2024-12-03 浙江大学 单模态医学图像的配准方法、系统及计算机可读存储介质
CN114677421B (zh) * 2022-04-12 2023-03-28 卡本(深圳)医疗器械有限公司 一种估算2d器官刚性/非刚性配准方法
CN114943753B (zh) * 2022-06-15 2024-07-23 北京理工大学 基于局部结构矢量对齐的位姿配准方法及装置
CN119600073B (zh) * 2024-11-21 2026-02-27 西安电子科技大学 具有辐射和旋转不变性的两阶段异源遥感图像配准方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004302863A (ja) * 2003-03-31 2004-10-28 Japan Research Institute Ltd データ管理装置、データ管理方法およびその方法をコンピュータに実行させるプログラム
CN101053531A (zh) * 2007-05-17 2007-10-17 上海交通大学 基于多模式增敏成像融合的早期肿瘤定位跟踪方法
CN101556692A (zh) * 2008-04-09 2009-10-14 西安盛泽电子有限公司 基于特征点邻域伪Zernike矩的图像拼接方法
CN103345741A (zh) * 2013-06-13 2013-10-09 华中科技大学 一种非刚性多模医学图像精确配准方法
CN103700101A (zh) * 2013-12-19 2014-04-02 华东师范大学 一种非刚性脑图像配准方法
CN104200460A (zh) * 2014-08-04 2014-12-10 西安电子科技大学 基于图像特征和互信息的图像配准方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7635832B2 (en) * 2005-08-31 2009-12-22 The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration Hybrid diversity method utilizing adaptive diversity function for recovering unknown aberrations in an optical system
US9122926B2 (en) * 2012-07-19 2015-09-01 Honeywell International Inc. Iris recognition using localized Zernike moments

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004302863A (ja) * 2003-03-31 2004-10-28 Japan Research Institute Ltd データ管理装置、データ管理方法およびその方法をコンピュータに実行させるプログラム
CN101053531A (zh) * 2007-05-17 2007-10-17 上海交通大学 基于多模式增敏成像融合的早期肿瘤定位跟踪方法
CN101556692A (zh) * 2008-04-09 2009-10-14 西安盛泽电子有限公司 基于特征点邻域伪Zernike矩的图像拼接方法
CN103345741A (zh) * 2013-06-13 2013-10-09 华中科技大学 一种非刚性多模医学图像精确配准方法
CN103700101A (zh) * 2013-12-19 2014-04-02 华东师范大学 一种非刚性脑图像配准方法
CN104200460A (zh) * 2014-08-04 2014-12-10 西安电子科技大学 基于图像特征和互信息的图像配准方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
LING ZHIGANG ET AL.: "A robust multi-source remote- Sensing image Registration Method Based on Feature Matching", CHINESE JOURNAL OF ELECTRONICS, vol. 38, no. 12, 31 December 2010 (2010-12-31), pages 2893 - 2894 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111754554A (zh) * 2020-06-28 2020-10-09 上海应用技术大学 颅脑多模态医学图像配准方法
CN111754554B (zh) * 2020-06-28 2023-09-15 上海应用技术大学 颅脑多模态医学图像配准方法
CN114581497A (zh) * 2022-04-16 2022-06-03 合肥学院 一种用于大气湍流图像畸变校正的改进b样条非刚性配准方法
CN114581497B (zh) * 2022-04-16 2024-09-10 合肥学院 一种用于大气湍流图像畸变校正的改进b样条非刚性配准方法
CN115496785A (zh) * 2022-08-17 2022-12-20 中科超精(南京)科技有限公司 一种自适应约束模型的医学影像配准方法
CN115546268A (zh) * 2022-09-23 2022-12-30 中国人民解放军国防科技大学 多模态遥感图像配准方法、系统、终端设备及存储介质

Also Published As

Publication number Publication date
US20190130572A1 (en) 2019-05-02
CN106204550B (zh) 2018-10-30
US10853941B2 (en) 2020-12-01
CN106204550A (zh) 2016-12-07

Similar Documents

Publication Publication Date Title
WO2018000652A1 (zh) 一种非刚性多模医学图像的配准方法及系统
CN108416802B (zh) 一种基于深度学习的多模医学图像非刚性配准方法及系统
CN109949349B (zh) 一种多模态三维图像的配准及融合显示方法
Li et al. Automated measurement network for accurate segmentation and parameter modification in fetal head ultrasound images
CN109523584B (zh) 图像处理方法、装置、多模态成像系统、存储介质及设备
CN107146228B (zh) 一种基于先验知识的大脑磁共振图像超体素生成方法
CN114972366A (zh) 基于图网络的大脑皮层表面全自动分割方法及系统
WO2024169341A1 (zh) 一种多模态影像引导放射治疗的配准方法
CN103400376B (zh) 一种乳腺动态增强磁共振图像序列的配准方法
Göçeri Fully automated liver segmentation using Sobolev gradient‐based level set evolution
CN110036409B (zh) 使用联合深度学习模型进行图像分割的系统和方法
CN115272261B (zh) 一种基于深度学习的多模态医学图像融合方法
CN115661219B (zh) 一种基于互信息及l-bfgs优化的三维ct/pet图像配准方法
CN110223331B (zh) 一种大脑mr医学图像配准方法
CN117197200A (zh) 基于深度概率图和矢量融合的腹部多器官配准方法
CN106504239A (zh) 一种提取超声图像中肝脏区域的方法
CN103345741B (zh) 一种非刚性多模医学图像精确配准方法
CN115409879A (zh) 图像配准的数据处理方法及装置、存储介质、电子设备
Perez-Gonzalez et al. Deep learning spatial compounding from multiple fetal head ultrasound acquisitions
CN116934712B (zh) 应用于肺部三维图像处理的配准方法及装置
CN110246566A (zh) 基于卷积神经网络的品行障碍确定方法、系统和存储介质
Wu et al. TPS-HAMMER: Improving HAMMER registration algorithm by soft correspondence matching and thin-plate splines based deformation interpolation
Zhang et al. A computational white matter atlas for aging with surface-based representation of fasciculi
Huang et al. [Retracted] Adoption of Snake Variable Model‐Based Method in Segmentation and Quantitative Calculation of Cardiac Ultrasound Medical Images
CN101305395A (zh) 基于点的自适应弹性图像配准

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 16907048

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 16907048

Country of ref document: EP

Kind code of ref document: A1