WO2012028053A1 - 一种图像配准方法 - Google Patents

一种图像配准方法 Download PDF

Info

Publication number
WO2012028053A1
WO2012028053A1 PCT/CN2011/078267 CN2011078267W WO2012028053A1 WO 2012028053 A1 WO2012028053 A1 WO 2012028053A1 CN 2011078267 W CN2011078267 W CN 2011078267W WO 2012028053 A1 WO2012028053 A1 WO 2012028053A1
Authority
WO
WIPO (PCT)
Prior art keywords
average distance
image
initial
registration
coordinate transformation
Prior art date
Application number
PCT/CN2011/078267
Other languages
English (en)
French (fr)
Inventor
孙毅勇
宫晶晶
刘道志
Original Assignee
微创医疗器械(上海)有限公司
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 微创医疗器械(上海)有限公司 filed Critical 微创医疗器械(上海)有限公司
Priority to US13/814,760 priority Critical patent/US9042614B2/en
Priority to ES11821079T priority patent/ES2739052T3/es
Priority to EP11821079.8A priority patent/EP2613293B1/en
Publication of WO2012028053A1 publication Critical patent/WO2012028053A1/zh

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/74—Image or video pattern matching; Proximity measures in feature spaces
    • G06V10/75—Organisation of the matching processes, e.g. simultaneous or sequential comparisons of image or video features; Coarse-fine approaches, e.g. multi-scale approaches; using context analysis; Selection of dictionaries
    • G06V10/754—Organisation of the matching processes, e.g. simultaneous or sequential comparisons of image or video features; Coarse-fine approaches, e.g. multi-scale approaches; using context analysis; Selection of dictionaries involving a deformation of the sample pattern or of the reference pattern; Elastic matching
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR 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; CALCULATING OR 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
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00—Scenes; Scene-specific elements
    • G06V20/60—Type of objects
    • G06V20/69—Microscopic objects, e.g. biological cells or cellular parts
    • G06V20/695—Preprocessing, e.g. image segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR 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; CALCULATING OR 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
    • G06T2207/10136—3D ultrasound image

Definitions

  • the present application relates to the field of image processing technologies, and in particular, to an image registration method based on an ICP algorithm. Background technique
  • the ICP (Iterative Closest Point) algorithm is a method for registering images.
  • the ICP algorithm is a registration method based on freeform surfaces. It is repeated to determine the corresponding point set and calculate the optimal rigid body variation. The process, until the preset registration convergence criterion is reached, the final coordinate change is the synthesis of each transformation.
  • the conventional ICP algorithm has a registration target and the initial offset range of the registration source cannot be too large, it can converge to a certain local limit, but it is usually not the global optimal (that is, it cannot break away from the local limit), so the image is matched.
  • the accuracy and success rate of the quasi-standard are not enough to meet the actual needs.
  • the embodiment of the present application provides an image registration method, so that registration of the registration source and the registration target in a large range of initial offset can be realized without manual adjustment of parameters, and Effectively get rid of the limitations of local limits, improve registration accuracy and registration success rate.
  • An image registration method includes:
  • the image ⁇ , ⁇ and the image ⁇ , ⁇ are iterated in the preset step length L to obtain the coordinate transformation between the image and the image ⁇ ' ⁇ ⁇ 3 ⁇ 4 ⁇ , 0 ⁇ 1 ⁇ L, calculate the average distance corresponding to it according to the coordinate transformation ⁇ 3 ⁇ 4 ⁇ ⁇ ⁇ ;
  • the sum is a column vector of 3*1.
  • the acquiring the closest point on the image ⁇ ⁇ corresponding to the point on the image ⁇ ⁇ ' ⁇ is specifically:
  • a KD binary tree is used to search for the closest point on the image ⁇ ⁇ corresponding to the point on the image ⁇ Pl , ⁇ .
  • the determining the average distance ⁇ ⁇ and the ideal average distance ⁇ ⁇ is specifically:
  • the determination of the average distance ⁇ and the average distance over ⁇ ⁇ particular body size as follows: Selecting a minimum average distance E mm of the average distance ⁇ ⁇ ;
  • the minimum average distance Emin and the ideal average distance E are determined.
  • the method further comprises: sequentially determining the average distance ⁇ ⁇ and an initial average distance E. Size, where: the initial average distance E. Transform H with the initial coordinates. Corresponding;
  • All values in the average distance ⁇ ⁇ are greater than the initial average distance E. Then, after the random disturbance control parameter ⁇ is attenuated according to a fixed attenuation rate, a random disturbance is applied to the image ⁇ Pl ⁇ , and the initial coordinate transformation is used to recalculate the average distance ⁇ " ⁇ , and the average is determined. Whether the distance ⁇ " ⁇ is less than or equal to the preset ideal average distance E step.
  • the method further includes: selecting a minimum average distance E mm of the average distance ⁇ ⁇ ;
  • the minimum average distance E mm and the initial average distance E are determined. Size, where: the initial average distance E. Transform H with the initial coordinates. Corresponding;
  • the coordinate transformation H min corresponding to the minimum average distance E mm is replaced by the initial coordinate transformation 3 ⁇ 4, and the random disturbance control parameter ⁇ is attenuated according to a fixed attenuation rate, according to the coordinate transformation H min and the attenuated random Disturbing the control parameter ⁇ , recalculating the average distance ⁇ , ⁇ , and determining whether the average distance ⁇ , ⁇ is less than or equal to the preset ideal average distance E step;
  • the fixed attenuation rate of the random disturbance control parameter ⁇ is .
  • the step size L ranges from 50 to 200. It can be seen from the technical solution provided by the above embodiment of the present application that the random disturbance control parameter ⁇ in the method can be automatically attenuated during the entire iterative process.
  • the preset initial transformation is the optimal starting position, there is no need to worry about transforming the model to a worse position because of the increase of random disturbance; when the initial transformation itself is not optimal, the model can be freed from local by increasing the random disturbance.
  • the extreme value is bound and transformed to a better position for iteration, that is, the global optimum can be approximated by searching.
  • the random disturbance control parameter ⁇ in the method is not set according to experience, but is comprehensively set after considering the global search performance and the iterative loop speed.
  • the random disturbance control parameters ⁇ can be set to the same value, and both have better registration effects.
  • the image registration method provided by the embodiment of the present invention can not only achieve the registration of the registration source and the registration target in a large range of initial offset, but also can effectively get rid of the situation without manually adjusting the parameters.
  • the constraint of local limits improves the registration accuracy and registration success rate.
  • Embodiment 1 is a part of the embodiments of the present application, and not all of the embodiments. All other embodiments obtained by a person of ordinary skill in the art based on the embodiments of the present application, without departing from the inventive work, should fall within the scope of the present application.
  • Embodiment 1 is a part of the embodiments of the present application, and not all of the embodiments. All other embodiments obtained by a person of ordinary skill in the art based on the embodiments of the present application, without departing from the inventive work, should fall within the scope of the present application.
  • FIG. 1 is a schematic flowchart diagram of an image registration method according to Embodiment 1 of the present application.
  • the method includes:
  • Step S101 Select the registration source image ⁇ pi ⁇ and the registration target image ⁇ qj ⁇ .
  • Step S102 Apply a random disturbance to the image ⁇ ⁇ according to the preset random disturbance control parameter ⁇ to generate a deformation, obtain an image ⁇ Pl ' ⁇ , and obtain an image corresponding to a point on the image ⁇ ⁇ ' ⁇ on the image ⁇ ⁇ Recently, the set of points obtained is the image ⁇ ' ⁇ .
  • a random perturbation is applied to the registration source image ⁇ ⁇ , specifically:
  • Step S103 According to the preset initial coordinate transformation 3 ⁇ 4, within the preset step length L, the image ⁇ Pl , ⁇ Iterative operation is performed with the image ⁇ ' ⁇ to obtain a coordinate transformation ⁇ 3 ⁇ 4 ⁇ between the image ⁇ ⁇ and the image ⁇ ' ⁇ , 0 ⁇ 1 ⁇ L, and the corresponding average distance ⁇ ⁇ is calculated according to the coordinate transformation ⁇ 3 ⁇ 4 ⁇ .
  • Step 104 Determine whether the average distance ⁇ ⁇ is greater than a preset ideal average distance E x , and if not, the registration is terminated.
  • Method 1 sequentially determine the average distance ⁇ ⁇ and the ideal average distance E; mode 2: in the average distance ⁇ ⁇ Select the minimum average distance E mm and then determine the minimum average distance E mm and the ideal average distance E.
  • Embodiment 2 :
  • FIG. 2 is a schematic flowchart diagram of an image registration method according to Embodiment 2 of the present application. As shown in FIG. 2, steps S201 to S204 in the embodiment of the present application are the same as steps S101 to S104 in the first embodiment, and will not be described again. In step S204, if the average distance ⁇ ⁇ is greater than the preset ideal average distance, E x , then the method may further comprise:
  • Step S205 It is sequentially determined whether the average distance ⁇ ⁇ is greater than the initial average distance E. , where: the initial average distance Eo corresponds to the initial coordinate transformation Ho;
  • Step S206 If no, the judgment between the average distance ⁇ ⁇ and the initial average distance Eo is terminated, and the coordinate transformation H g corresponding to the current average distance E g is replaced with the initial coordinate transformation 3 ⁇ 4, and will be random according to the fixed attenuation rate.
  • the disturbance control parameter ⁇ is attenuated, and the average distance ⁇ , ⁇ is recalculated according to the coordinate transformation H g and the attenuated random disturbance control parameter ⁇ , and it is determined whether the average distance ⁇ , ⁇ is less than or equal to the preset ideal average distance E step.
  • Step S207 If yes, attenuate the random disturbance control parameter ⁇ according to a fixed decay rate, re-apply random disturbance to the image ⁇ ⁇ , and use the initial coordinate transformation ⁇ . , Recalculate the average distance ⁇ 3 ⁇ 4" ⁇ and determine if the average distance ⁇ " ⁇ is less than or equal to the preset ideal average distance E step.
  • the global search performance of the algorithm is inversely related to the iteration speed.
  • the relationship between the global search performance and the iteration speed can be balanced by setting the range of the random disturbance control parameter ⁇ and setting the step length L.
  • the preset random disturbance control is preset.
  • the parameter ⁇ is 10% of the geometric ratio of the registration source image ⁇ ⁇ or the registration target image ⁇ ⁇ .
  • the step size L ranges from 50 to 200.
  • FIG. 3 is a schematic flowchart diagram of an image registration method according to Embodiment 3 of the present application.
  • steps S301 to S304 in the embodiment of the present application are the same as steps S101 to S104 in the first embodiment, and will not be described again.
  • step S304 the average distance ⁇ ⁇ and the preset ideal average distance are determined.
  • the method can also include:
  • Step S305 Select a minimum average distance Emin in the average distance ⁇ ⁇ ;
  • Step S306 It is judged whether the minimum average distance Emin is greater than the initial average distance E. , where: the initial average distance Eo corresponds to the initial coordinate transformation Ho;
  • Step S307 Random disturbance control parameter ⁇ , recalculate the average distance ⁇ , ⁇ , and judge Whether the average distance ⁇ , ⁇ is less than or equal to the preset ideal average distance E step;
  • Step S308 If yes, the random disturbance control parameter ⁇ is attenuated according to the fixed attenuation rate, and then the random disturbance is applied to the image ⁇ ⁇ , and the initial coordinate transformation is used. , Recalculate the average distance ⁇ 3 ⁇ 4" ⁇ and determine if the average distance ⁇ " ⁇ is less than or equal to the preset ideal average distance E step.
  • Step S305 and step S 2 05 The only difference is how to determine the average distance from the initial average distance ⁇ ⁇ E.
  • FIG. 4 is a schematic flowchart diagram of an image registration method according to Embodiment 4 of the present application.
  • FIG. 5 is a schematic flow chart of another image registration method according to Embodiment 4 of the present application. As shown in FIG. 4 and FIG. 5, after the step S206 or S207 in the second embodiment, or after the third step S307 or S308, the method may further include: the registration is terminated, and an error is reported.
  • the preset iteration loop number K is to control the iterative loop within a reasonable range, avoiding a large amount of time to perform calculations without unregistration.
  • the preset number of iteration cycles K is 100.
  • the random disturbance control parameter ⁇ in the method can be automatically attenuated during the whole iterative process.
  • the preset initial transformation is the optimal starting position, there is no need to worry because the randomization is increased. Disturbing to transform the model to a worse position; and when the initial transformation itself is not optimal, it can also be iterated by increasing the random disturbance to make the model break away from the local extremum and transform to a better position for iteration. Close to global optimum.
  • the random disturbance control parameter ⁇ in the method is not set according to experience, but is comprehensively set after considering the global search performance and the iterative cycle speed.
  • the random disturbance control parameter ⁇ can be set to the same Values, and all get better registration results.
  • the image registration method provided by the embodiment of the present application does not need to manually adjust parameters.
  • the registration source and the registration target are registered in a large range of initial offset, but also can effectively get rid of the local limit, and improve the registration accuracy and registration success rate.

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Health & Medical Sciences (AREA)
  • Multimedia (AREA)
  • General Health & Medical Sciences (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Medical Informatics (AREA)
  • Evolutionary Computation (AREA)
  • Computing Systems (AREA)
  • Artificial Intelligence (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Biomedical Technology (AREA)
  • Molecular Biology (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)

Description

说 明 书 一种图像配准方法 技术领域
本申请涉及图像处理技术领域， 特别是涉及一种基于 ICP算法的图像配 准方法。 背景技术
随着新型传感器的不断涌现， 人们获取图像的能力迅速提高， 不同物理 特性的传感器所产生的图像也不断增多。 由于不同的图像传感器获取的图像 数据存在明显的局限性和差异性， 所以仅仅利用一种图像数据往往 4艮难满足 实际需求， 为此需要通过图像融合技术将不同传感器获取的图像总和起来使 用， 以达到对图像中的目标更加全面清晰、 准确的理解和认识。 例如在医学 上， 为了通过对解剖和生理信息进行综合分析而通过对不同形式诸如计算机 断层扫描成像（CT )、 核磁共振成像（MRI )、 超声 （US )获得图像进行融合， 以实现改进诊断过程。
图像配准技术是实现图像融合的重要前提， 是图像融合首先要解决的问 题。 ICP ( Iterative Closest Point, 迭代最近点）算法是实现对图像进行配准的 一种方法， ICP算法是一种基于自由形态曲面的配准方法， 重复进行 "确定对 应点集， 计算最优刚体变化" 的过程， 直到达到预设的配准收敛准则， 最终 坐标变化是每次变换的合成。 但由于传统的 ICP算法存在配准目标和配准源 的初始偏移范围不能太大， 可以收敛到某一局部极限， 但通常不是全局最优 (即无法挣脱局部极限） 的缺点， 因此图像配准的精度和成功率较低， 无法 满足实际需求。
通过对现有技术研究， 申请人发现： 现有的基于 ICP算法的图像配准方 法通常采用 "给配准施加随机刚体扰动变换" 和 "给配准目标的三维坐标点 施加随机扰动使其产生形变" 两种方法达到挣脱传统 ICP算法局部极限， 使 配准效果更好。 但采用 "给配准施加随机刚体扰动变换" 这种方法需要在六 维自由度空间采样， 计算量大且费时； 而采用 "给配准目标的三维坐标点施 加随机扰动使其产生形变" 方法， 其参数是在实际应用中根据不同的实验模 型靠经验进行设定的， 对于不同的实验模型， 参数设定比较困难， 设置错误 就可能导致配准失败。 发明内容
有鉴于此， 本申请实施例提供了一种图像配准方法， 以实现无需手动调 节参数， 即可实现在配准源和配准目标在大范围初始偏移的情况下进行配准， 而且能够有效地摆脱局部极限的束缚， 提高了配准精度和配准成功率。
为实现上述目的， 技术方案如下：
一种图像配准方法， 包括：
选取配准源图像 {pi}和配准目标图像 其中： 0<i<M, 0<j<N, M, N 为正整数， 且 M«N;
根据预设的随机扰动控制参数 σ, 对所述图像 { }施加随机扰动使其产生 形变， 得到图像 {ρ }, 并获取所述图像 { }上与图像 上的点相对应的最近 点， 得到的点的集合即为图像 { '};
根据预设的初始坐标变换 ¾,在预设步长 L范围内，对图像 { ，}与图像 { ，} 进行迭代运算， 得到图像 与图像 { '}之间的坐标变换 {¾}, 0<1<L, 根据 坐标变换 {¾}计算与其相对应的平均距离 { }；
判断所述平均距离 { }和所述预设的理想平均距离 E 大小， 当所述平均 距离 { }小于等于预设的理想平均距离时， 配准终止。
优选地， 所述 和 为 3*1的列向量。
优选地， 所述获取所述图像 { }上与图像 {Ρι'}上的点相对应的最近点， 具 体为：
运用 KD二叉树搜索所述图像 { }上与图像 {Pl，}上的点相对应的最近点。 优选地， 所述判断所述平均距离 { }和所述理想平均距离 Εχ的大小具体 为：
依次判断所述平均距离 { }是否小于等于所述理想平均距离 Εχ。
优选地，所述判断所述平均距离 { }和所述理想平均距离 Εχ的大小具体体 为： 选择所述平均距离 { }中的最小平均距离 Emm；
判断所述最小平均距离 Emin和所述理想平均距离 E 大小。
优选地， 当所述平均距离 { }大于预设的理想平均距离时， 还包括： 依次判断所述平均距离 { }和初始平均距离 E。的大小， 其中： 所述初始平 均距离 E。 与所述初始坐标变换 H。相对应；
当所述平均距离 { }中任意一个值小于等于所述初始平均距离 E。，则终止 所述平均距离 { }与所述初始平均距离 E。之间的判断， 并将与当前平均距离 Eg 相对应的坐标变换 Hg替换初始坐标变换 H。， 且按照固定衰减速率将所述随机 扰动控制参数 σ进行衰减，根据坐标变换 Hg和所述衰减后的随机扰动控制参数 σ , 重新计算平均距离 { ，} , 并判断所述平均距离 { ，}是否小于等于所述预 设的理想平均距离 E 步骤；
当所述平均距离 { }中所有值均大于所述初始平均距离 E。时，则按照固定 衰减速率将所述随机扰动控制参数 σ衰减后， 重新对所述图像 {Pl}施加随机扰 动， 并利用初始坐标变换 ¾, 重新计算平均距离 { "} , 并判断所述平均距离 { "}是否小于等于所述预设的理想平均距离 E 步骤。
优选地， 当所述平均距离 { }大于预设的理想平均距离时， 还包括： 选择所述平均距离 { }中的最小平均距离 Emm；
判断所述最小平均距离 Emm和初始平均距离 E。的大小， 其中： 所述初始平 均距离 E。 与所述初始坐标变换 H。相对应；
当所述最小平均距离 Emm小于等于所述初始平均距离 E。， 则将与最小平均 距离 Emm相对应的坐标变换 Hmin替换初始坐标变换 ¾, 且按照固定衰减速率将 所述随机扰动控制参数 σ进行衰减，根据坐标变换 Hmin和所述衰减后的随机扰 动控制参数 σ , 重新计算平均距离 { ，} , 并判断所述平均距离 { ，}是否小于 等于所述预设的理想平均距离 E 步骤；
当所述最小平均距离 Emm大于所述初始平均距离 E。时， 则按照固定衰减速 率将所述随机扰动控制参数 σ衰减后， 重新对所述图像 { }施加随机扰动， 并 利用初始坐标变换 Η。， 重新计算平均距离 { "} , 并判断所述平均距离 { "} 是否小于等于所述预设的理想平均距离 E 步骤。
优选地， 所述重新计算平均距离 { ，}或重新计算平均距离 { "}后， 还包 括： 环次数大于预设的迭代循环次数 K时， 则配准终止， 并报错。
优选地， 所述随机扰动控制参 数 σ的固定衰减速率为 。 优选地， 所述步长 L的范围为 50〜200。 由以上本申请实施例提供的技术方案可见， 该方法中的随机扰动控制参 数 σ在整个迭代过程中可以自动衰减。 当预设的初始变换是最佳起始位置时， 无需担心因为增加了随机扰动而将模型变换到更差的位置； 当初始变换本身 不是最佳时， 还可以通过增加随机扰动使模型挣脱局部极值的束缚而变换到 更好的位置进行迭代， 即通过搜索可以接近全局最优。 同时， 该方法中随机 扰动控制参数 σ不是根据经验设定的，而是考虑了全局搜索性能与迭代循环速 度后综合设定的。对于多组心腔模型，随机扰动控制参数 σ可以设置为同一值， 且都得到较好的配准效果。
因此， 本申请实施例提供的图像配准方法， 在无需手动调节参数的情况 下， 不仅可以实现配准源和配准目标在大范围初始偏移的情况下进行配准， 而且能够有效地摆脱局部极限的束缚， 提高了配准精度和配准成功率。
附图说明
为了更清楚地说明本申请实施例或现有技术中的技术方案， 下面将对实 施例或现有技术描述中所需要使用的附图作简单地介绍， 显而易见地， 下面 描述中的附图仅仅是本申请中记载的一些实施例， 对于本领域普通技术人员 来讲， 在不付出创造性劳动的前提下， 还可以根据这些附图获得其他的附图。 具体实施方式
为了使本技术领域的人员更好地理解本申请中的技术方案， 下面将结合 本申请实施例中的附图， 对本申请实施例中的技术方案进行清楚、 完整地描 述， 显然， 所描述的实施例仅仅是本申请一部分实施例， 而不是全部的实施 例。 基于本申请中的实施例， 本领域普通技术人员在没有做出创造性劳动前 提下所获得的所有其他实施例， 都应当属于本申请保护的范围。 实施例一：
图 1为本申请实施例一提供的一种图像配准方法的流程示意图。
如图 1所示， 该方法包括：
步骤 S101: 选取配准源图像 {pi}和配准目标图像 {qj}。
配准源图像为 {pi I i=l, ...... ,M}， 配准目标图像为 {qj I j=l, ...... ,N}，
M«N, 配准源图像 {pi}和配准目标图像 {qj}处于不同的坐标系，
Figure imgf000007_0001
其中： 0<i<M, 0<j<N, M, N为正整数， 且 M«N; R为 3*3的旋转矩 阵， T为 3*1的平移向量；
步骤 S102: 根据预设的随机扰动控制参数 σ, 对图像 { }施加随机扰动使 其产生形变， 得到图像 {Pl'}， 并获取图像 { }上与图像 {Ρι'}上的点相对应的最 近点， 得到的点的集合即为图像 { '}。
首先对配准源图像 { }施加随机扰动， 具体为：
给配准源图像 {Pl}中每个点增加满足（0, σ1 )正态分布的随机扰动 nl 使得配准源图像 {Pl}产生形变，
即 Pi'= i+ni ,
其中， 为 3*1的列向量。
然后， 运用 KD二叉树搜索配准目标图像像 { =i, ...... ,Ν}与形变后的配 准源图像 {Pl'|i=l, ...... ,Μ}上相对应的最近点， 搜索得到的点的集合即为
(Φ'ϋ=ι,…… ,Ν}
并根据 =| ,''(；,'')' , 计算 3*3矩阵 X, 并根据 X=UDV对矩阵 X进行奇 异值分解， 其中： p"= p'-p , q"=q'-q , ρ = - pf , q =丄^ qf , U和 V是 3*3的正交矩 阵， D为由奇异值组成的对角矩阵。
最后， 分别计算旋转矩阵 R和平移向量 T, 其中： R=VUtrT="^-Rq 。 步骤 S103:根据预设的初始坐标变换 ¾,在预设步长 L范围内，对图像 {Pl，} 与图像 { ' }进行迭代运算， 得到图像 {ρ }与图像 { ' }之间的坐标变换 {¾} , 0<1<L, 根据坐标变换 {¾}计算与其相对应的平均距离 { }。
R 0
根据 {Ρι}=Η*{ }, JLH: , 计算 }和{ }间的坐标变换 Η,
Τ 1
并根据坐标变换 Η计算 {Pl}和 { }间的平均距离 Ε , 其中：
Figure imgf000008_0001
并且对图像 {Pl，}与图像 { ，}进行 L次迭代运算， 并记录 L次迭代过程中图 像 和图像 { }之间的坐标变换 {¾| i=l,… … ,L}， 并根据坐标变换 {¾| i=l, ...... ,L}计算得到图像 {Pl}和图像 { }之间平均距离 { |i=l, ...... ,L}， 其中：
0<1<L。
步骤 104: 判断平均距离 { }是否大于预设的理想平均距离 Ex, 如果否， 配准终止。
如果平均距离 { }的中一个值小于等于预设的理想平均距离 Ex,则与平均 距离 { }的中的这个值即为理想坐标变换， 所以配准可以终止。
在本申请实施例中，判断平均距离 { }和预设的理想平均距离 E 大小有 两种方式，方式一：依次判断平均距离 { }和理想平均距离 E 大小； 方式二： 在平均距离 { }中选择最小平均距离 Emm, 然后判断最小平均距离 Emm和理想 平均距离 E 大小。 实施例二：
图 2为本申请实施例二提供的一种图像配准方法的流程示意图。 如图 2所示， 本申请实施例中步骤 S201〜步骤 S204与实施例一中的步骤 S101〜步骤 S104相同， 再次不再贅述， 在步骤 S204 中如果平均距离 { }大于 预设的理想平均距离 Ex, 则该方法还可以包括：
步骤 S205: 依次判断平均距离 { }是否大于初始平均距离 E。， 其中： 初始 平均距离 Eo 与初始坐标变换 Ho相对应；
步骤 S206: 如果否， 则终止平均距离 { }与初始平均距离 Eo之间的判断， 并将与当前平均距离 Eg相对应的坐标变换 Hg替换初始坐标变换 ¾, 且按照固 定衰减速率将随机扰动控制参数 σ进行衰减，根据坐标变换 Hg和衰减后的随机 扰动控制参数 σ , 重新计算平均距离 { ，} , 并判断平均距离 { ，}是否小于等 于预设的理想平均距离 E 步骤。
步骤 S207: 如果是， 则按照固定衰减速率将随机扰动控制参数 σ衰减后， 重新对图像 { }施加随机扰动， 并利用初始坐标变换 Η。， 重新计算平均距离 {¾"} , 并判断平均距离 { "}是否小于等于预设的理想平均距离 E 步骤。
算法的全局搜索性能与迭代速度呈反关系通过设置随机扰动控制参数 σ 的范围和设置步长 L可以平衡全局搜索性能和迭代速度之间的关系，在本申请 实施例中，预设随机扰动控制参数 σ的为配准源图像 { }或者配准目标图像 { } 几何比例的 10%。 步长 L的范围为 50〜200。 实施例三：
图 3为本申请实施例三提供的一种图像配准方法的流程示意图。
如图 3 所示， 本申请实施例中步骤 S301〜步骤 S304与实施例一中的步骤 S101〜步骤 S104相同， 再次不再贅述， 在步骤 S304 中判断平均距离 { }和预 设的理想平均距离 E 大小之后， 该方法还可以包括：
步骤 S305： 选择平均距离 { }中的最小平均距离 Emin；
步骤 S306: 判断最小平均距离 Emin是否大于初始平均距离 E。， 其中： 初始 平均距离 Eo 与初始坐标变换 Ho相对应；
步骤 S307: 如果否， 则将与最小平均距离 Emm相对应的坐标变换 Hmin替换 初始坐标变换 ¾,且按照固定衰减速率将随机扰动控制参数 σ进行衰减，根据 坐标变换 Hmm和衰减后的随机扰动控制参数 σ , 重新计算平均距离 { ，} , 并判 断平均距离 { ，}是否小于等于预设的理想平均距离 E 步骤；
步骤 S308: 如果是， 则按照固定衰减速率将随机扰动控制参数 σ衰减后， 重新对图像 { }施加随机扰动， 并利用初始坐标变换 Η。， 重新计算平均距离 {¾"} , 并判断平均距离 { "}是否小于等于预设的理想平均距离 E 步骤。
步骤 S305和步骤 S205所不同的只是如何判断平均距离 { }与初始平均距 离 E。的大小， 步骤 S205中， 将平均距离 { }依次与初始平均距离 E。进行比较， 而步骤 S305中， 则是首先选取平均距离 { }中的最小平均距离 Emin, 然后用最 小平均距离 Emm与初始平均距离 E0进行比较。 实施例四：
图 4为本申请实施例四提供的一种图像配准方法的流程示意图。 图 5为 本申请实施例四提供的另一种图像配准方法的流程示意图。 如图 4和图 5所 示，在实施例二中步骤 S206或 S207之后，或者在实施例三步骤 S307或 S308 之后， 该方法还可以包括： 则配准终止， 并报错。
预设迭代循环次数 K, 是为了将迭代循环控制在合理的范围内， 避免无 法配准的情况下， 浪费大量的时间进行运算。 在本申请实施例中， 预设的迭 代循环次数 K选择 100。 由以上本申请实施例提供的技术方案可见， 该方法中的随机扰动控制参 数 σ在整个迭代过程中可以自动衰减， 当预设的初始变换是最佳起始位置时， 无需担心因为增加了随机扰动而将模型变换到更差的位置； 并且当初始变换 本身不是最佳时， 还可以通过增加随机扰动可以使模型挣脱局部极值的束缚 而变换到更好的位置进行迭代， 即可以搜索可以接近全局最优。 同时， 该方 法中随机扰动控制参数 σ不是根据经验设定的 ,而是考虑了全局搜索性能与迭 代循环速度后综合设定的，对于多组心腔模型， 随机扰动控制参数 σ可以设置 为同一值 , 且都得到较好的配准效果。
因此， 本申请实施例提供的图像配准方法， 在无需手动调节参数的情况 下， 不仅可以实现在配准源和配准目标在大范围初始偏移的情况下进行配准， 而且能够有效地摆脱局部极限的束缚， 提高了配准精度和配准成功率。
以上仅是本申请的优选实施方式， 使本领域技术人员能够理解或实现本 申请。 对这些实施例的多种修改对本领域的技术人员来说将是显而易见的， 本文中所定义的一般原理可以在不脱离本申请的精神或范围的情况下， 在其 它实施例中实现。 因此， 本申请将不会被限制于本文所示的这些实施例， 而 是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims

权 利 要 求 书
1、 一种图像配准方法， 其特征在于， 包括：
选取配准源图像 {pi}和配准目标图像 其中： 0<i<M, 0<j<N, M, N 为正整数， 且 M«N;
根据预设的随机扰动控制参数 σ , 对所述图像 { }施加随机扰动使其产生 形变， 得到图像 并获取所述图像 { }上与图像 上的点相对应的最近 点， 得到的点的集合即为图像 { '} ;
根据预设的初始坐标变换 ¾,在预设步长 L范围内，对图像 { ，}与图像 { ，} 进行迭代运算， 得到图像 { ，}与图像 { ，}之间的坐标变换 {¾} , 0<1<L, 根据 坐标变换 {¾}计算与其相对应的平均距离 { }；
判断所述平均距离 { }和所述预设的理想平均距离 E 大小， 当所述平均 距离 { }小于等于预设的理想平均距离时， 配准终止。
2、根据权利要求 1所述的方法，其特征在于，所述 和 为 3*1的列向量。
3、 根据权利要求 2所述的方法， 其特征在于， 所述获取所述图像 { }上 与图像 {Ρι'}上的点相对应的最近点， 具体为：
运用 KD二叉树搜索所述图像 { }上与图像 {Pl，}上的点相对应的最近点。
4、根据权利要求 3所述的方法，其特征在于，所述判断所述平均距离 { } 和所述理想平均距离 E 大小具体为：
依次判断所述平均距离 { }是否小于等于所述理想平均距离 Ex。
5、根据权利要求 3所述的方法，其特征在于，所述判断所述平均距离 { } 和所述理想平均距离 E 大小具体为：
选择所述平均距离 { }中的最小平均距离 Emm；
判断所述最小平均距离 Emin和所述理想平均距离 E 大小。
6、 根据权利要求 1所述的方法， 其特征在于， 当所述平均距离 { }大于 预设的理想平均距离时， 还包括：
依次判断所述平均距离 { }和初始平均距离 E。的大小， 其中： 所述初始平 均距离 E。 与所述初始坐标变换 H。相对应；
当所述平均距离 { }中任意一个值小于等于所述初始平均距离 E。，则终止 所述平均距离 { }与所述初始平均距离 E。之间的判断， 并将与当前平均距离 Eg 相对应的坐标变换 Hs替换初始坐标变换 HQ, 且按照固定衰减速率将所述随机 扰动控制参数 σ进行衰减，根据坐标变换 Hg和所述衰减后的随机扰动控制参数 σ , 重新计算平均距离 { ，} , 并判断所述平均距离 { ，}是否小于等于所述预 设的理想平均距离 E 步骤；
当所述平均距离 { }中所有值均大于所述初始平均距离 E。时，则按照固定 衰减速率将所述随机扰动控制参数 σ衰减后， 重新对所述图像 {Pl}施加随机扰 动， 并利用初始坐标变换 ¾, 重新计算平均距离 { "} , 并判断所述平均距离 { "}是否小于等于所述预设的理想平均距离 E 步骤。
7、 根据权利要求 1所述的方法， 其特征在于， 当所述平均距离 { }大于 预设的理想平均距离时， 还包括：
选择所述平均距离 { }中的最小平均距离 Emm；
判断所述最小平均距离 Emm和初始平均距离 E。的大小， 其中： 所述初始平 均距离 E。 与所述初始坐标变换 H。相对应；
当所述最小平均距离 Emm小于等于所述初始平均距离 E。， 则将与最小平均 距离 Emin相对应的坐标变换 Hmin替换初始坐标变换 Η。， 且按照固定衰减速率将 所述随机扰动控制参数 σ进行衰减，根据坐标变换 Hmin和所述衰减后的随机扰 动控制参数 σ , 重新计算平均距离 { ，} , 并判断所述平均距离 { ，}是否小于 等于所述预设的理想平均距离 E 步骤；
当所述最小平均距离 Emm大于所述初始平均距离 E。时， 则按照固定衰减速 率将所述随机扰动控制参数 σ衰减后， 重新对所述图像 { }施加随机扰动， 并 利用初始坐标变换 Η。， 重新计算平均距离 { "} , 并判断所述平均距离 { "} 是否小于等于所述预设的理想平均距离 E 步骤。
8、 根据权利要求 6或 7所述的方法， 其特征在于， 所述重新计算平均距 离 { ，}或重新计算平均距离 { "}后， 还包括： 环次数大于预设的迭代循环次数 K时， 则配准终止， 并报错。
9、 根据权利要求 6或 7所述的方法， 其特征在于， 所述随机扰动控制参 数 σ的固定衰减速率为 。
10、 根据权利要求 1-7任一项所述的方法， 其特征在于， 所述步长 L的 范围为 50-200„
PCT/CN2011/078267 2010-08-31 2011-08-11 一种图像配准方法 WO2012028053A1 (zh)

Priority Applications (3)

Application Number Priority Date Filing Date Title
US13/814,760 US9042614B2 (en) 2010-08-31 2011-08-11 Image registration method
ES11821079T ES2739052T3 (es) 2010-08-31 2011-08-11 Método de registro de imágenes
EP11821079.8A EP2613293B1 (en) 2010-08-31 2011-08-11 Image registration method

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN2010102701489A CN102385748B (zh) 2010-08-31 2010-08-31 一种图像配准方法
CN201010270148.9 2010-08-31

Publications (1)

Publication Number Publication Date
WO2012028053A1 true WO2012028053A1 (zh) 2012-03-08

Family

ID=45772146

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2011/078267 WO2012028053A1 (zh) 2010-08-31 2011-08-11 一种图像配准方法

Country Status (5)

Country Link
US (1) US9042614B2 (zh)
EP (1) EP2613293B1 (zh)
CN (1) CN102385748B (zh)
ES (1) ES2739052T3 (zh)
WO (1) WO2012028053A1 (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110136177A (zh) * 2018-02-08 2019-08-16 北京连心医疗科技有限公司 一种图像配准方法、设备和存储介质
CN113570648A (zh) * 2021-07-30 2021-10-29 武汉联影智融医疗科技有限公司 多骨骼影像配准方法、电子装置以及医学导航系统

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104966287B (zh) * 2015-06-08 2017-08-08 浙江大学 层次化的多片点云刚性配准方法
US20170337682A1 (en) 2016-05-18 2017-11-23 Siemens Healthcare Gmbh Method and System for Image Registration Using an Intelligent Artificial Agent
CN110368027B (zh) * 2018-04-13 2022-02-18 北京柏惠维康科技有限公司 一种图像融合方法和装置
CN112950706B (zh) * 2019-12-11 2023-11-03 京东科技信息技术有限公司 可移动设备定位数据处理方法、装置、设备及存储介质
CN111899286B (zh) * 2020-07-13 2024-01-05 江西理工大学 基于精英差分演化的图像配准方法
CN112465883B (zh) * 2020-11-23 2022-03-29 山东科技大学 一种高精度曲面非均匀图像配准方法
CN113558766B (zh) * 2021-07-19 2022-05-17 北京纳通医学研究院有限公司 图像配准方法、装置、手术机器人以及手术机器人系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1489101A (zh) * 2003-08-14 2004-04-14 中国人民解放军第一军医大学 基于相同重叠区框架下的多设备医学图像刚性配准方法
CN1818974A (zh) * 2006-03-08 2006-08-16 杭州电子科技大学 一种多模态医学体数据三维可视化方法
US20090046951A1 (en) * 2007-07-16 2009-02-19 Nikos Paragios System and Method for Dense Image Registration Using Markov Random Fields and Efficient Linear Programming

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7177486B2 (en) * 2002-04-08 2007-02-13 Rensselaer Polytechnic Institute Dual bootstrap iterative closest point method and algorithm for image registration
US7187809B2 (en) * 2004-06-10 2007-03-06 Sarnoff Corporation Method and apparatus for aligning video to three-dimensional point clouds
US7715604B2 (en) * 2005-01-18 2010-05-11 Siemens Medical Solutions Usa, Inc. System and method for automatically registering three dimensional cardiac images with electro-anatomical cardiac mapping data
US7583857B2 (en) * 2005-08-24 2009-09-01 Siemens Medical Solutions Usa, Inc. System and method for salient region feature based 3D multi modality registration of medical images
AU2006302057B2 (en) * 2005-10-11 2013-03-21 Carnegie Mellon University Sensor guided catheter navigation system
EP1780672A1 (en) * 2005-10-25 2007-05-02 Bracco Imaging, S.P.A. Method of registering images, algorithm for carrying out the method of registering images, a program for registering images using the said algorithm and a method of treating biomedical images to reduce imaging artefacts caused by object movement
US8064731B2 (en) * 2006-11-13 2011-11-22 Siemens Audiologische Technik Gmbh Generalized rigid alignment of 3D ear impression models
KR100884066B1 (ko) * 2007-03-30 2009-02-19 한국전자통신연구원 Ｓｖｄ 기반의 영상 비교시스템 및 방법
CN100559398C (zh) * 2007-06-19 2009-11-11 北京航空航天大学 自动的深度图像配准方法
US9025841B2 (en) * 2009-11-18 2015-05-05 Siemens Aktiengesellschaft Method and system for segmentation of the prostate in 3D magnetic resonance images

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1489101A (zh) * 2003-08-14 2004-04-14 中国人民解放军第一军医大学 基于相同重叠区框架下的多设备医学图像刚性配准方法
CN1818974A (zh) * 2006-03-08 2006-08-16 杭州电子科技大学 一种多模态医学体数据三维可视化方法
US20090046951A1 (en) * 2007-07-16 2009-02-19 Nikos Paragios System and Method for Dense Image Registration Using Markov Random Fields and Efficient Linear Programming

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
See also references of EP2613293A4 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110136177A (zh) * 2018-02-08 2019-08-16 北京连心医疗科技有限公司 一种图像配准方法、设备和存储介质
CN110136177B (zh) * 2018-02-08 2022-02-18 北京连心医疗科技有限公司 一种图像配准方法、设备和存储介质
CN113570648A (zh) * 2021-07-30 2021-10-29 武汉联影智融医疗科技有限公司 多骨骼影像配准方法、电子装置以及医学导航系统
CN113570648B (zh) * 2021-07-30 2023-09-26 武汉联影智融医疗科技有限公司 多骨骼影像配准方法、电子装置以及医学导航系统

Also Published As

Publication number Publication date
US9042614B2 (en) 2015-05-26
CN102385748B (zh) 2013-12-25
EP2613293B1 (en) 2019-06-12
CN102385748A (zh) 2012-03-21
EP2613293A1 (en) 2013-07-10
US20130156281A1 (en) 2013-06-20
EP2613293A4 (en) 2017-07-26
ES2739052T3 (es) 2020-01-28

Similar Documents

Publication Publication Date Title
WO2012028053A1 (zh) 一种图像配准方法
Vandemeulebroucke et al. Spatiotemporal motion estimation for respiratory‐correlated imaging of the lungs
AU2017281290A1 (en) Method for estimating at least one of shape, position and orientation of a dental restoration
Guetter et al. Efficient symmetric and inverse-consistent deformable registration through interleaved optimization
CN112614169B (zh) 基于深度学习网络的2d/3d脊椎ct层级配准方法
WO2020186371A1 (zh) 3d/2d血管配准方法及装置
JP2017509411A (ja) データ及びモデルを用いた画像再構成ならびに補正のためのシステム及び方法
CN106133789B (zh) 用于分割感兴趣区域的图像处理装置和方法
CN108198235B (zh) 一种三维超声重建方法、装置、设备及存储介质
US20140032197A1 (en) Method and apparatus for creating model of patient specified target organ based on blood vessel structure
US20150371372A1 (en) System and method for medical image quality enhancement using multiscale total variation flow
US20150254841A1 (en) Image processing device, imaging system, and image processing program
CN110288637B (zh) 多角度dsa造影图像血管匹配方法及装置
US9092666B2 (en) Method and apparatus for estimating organ deformation model and medical image system
JP2023507865A (ja) 医用画像のセグメント化
CN115601243A (zh) 图像畸变矫正方法、装置、设备及计算机可读存储介质
CN108428245B (zh) 基于自适应正则项的滑移图像配准方法
CN111402221A (zh) 一种图像处理方法、装置及电子设备
CN111161330B (zh) 非刚性图像配准方法、装置、系统、电子设备、存储介质
JP7459243B2 (ja) １以上のニューラルネットワークとしての画像形成のモデル化による画像再構成
WO2015065781A2 (en) System and method for model-based reconstruction of quantitative images
CN109919916B (zh) 一种壁面切应力优化方法及装置、存储介质
US8103075B2 (en) Method for providing extended possibilities when imaging a patient&#39;s heart
Bergvall et al. A fast and highly automated approach to myocardial motion analysis using phase contrast magnetic resonance imaging
CN115661219A (zh) 一种基于互信息及l-bfgs优化的三维ct/pet图像配准方法

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: 11821079

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

WWE Wipo information: entry into national phase

Ref document number: 13814760

Country of ref document: US

WWE Wipo information: entry into national phase

Ref document number: 2011821079

Country of ref document: EP