专栏首页最新医学影像技术医学图像处理案例(二十一)——一致性点漂移算法(Coherent Point Drift)

医学图像处理案例(二十一)——一致性点漂移算法(Coherent Point Drift)

今天将分享使用一致性点漂移算法(Coherent Point Drift)来对点云数据进行配准。

1、常见配准变换

简单介绍一下什么是配准,其实很简单就是给定一个数据集,在另一个数据集中找到其所对应数据集的位置。

配准常见的有刚性变换配准,仿射变换配准和非刚性变换配准,其中刚性变换涉及到平移,旋转和各向同性缩放,仿射变换涉及到平移,旋转和各向异性缩放,倾斜等,非刚性变换就会很复杂有很多变换参数,例如分段式仿射变换,多项式模型等。

2、一致性点漂移算法简介

在论文《Point-Set Registration: Coherent Point Drift》中介绍了刚性和非刚性的配准概率方法,称为一致性点漂移算法。论文将点集之间的对应关系看成是一种概率关系问题,在最理想情况下真正的对应点的概率是1,错误的对应点概率值是0,所以直接避开点之间的坐标关系而转向使用概率值来描述这种对应关系,如果概率值越大,那么这种对应关系的确定性也就越大。

既然涉及到了概率,那就会考虑到概率模型,常见有均匀分布、二项分布、正态分布等,此时需要选择一个合适的概率模型来描述这种对应关系。如果针对一个点与一个点的对应关系,可以使用正态分布(高斯分布)来描述,所以当存在多个点时,就需要使用混合高斯模型(Gaussian Mixture Model)来进行描述了。

实际上两个点集的配准问题就是概率密度估计问题,因为其中一个点集合代表了高斯混合模型的质心,而另一个点集合代表了数据点,所以两个点集的配准问题就转化成一个概率密度估计的问题,也就转换成了求解混合高斯模型的参数问题了。

使用迭代最大化期望(最大似然法)来拟合高斯混合模型的质心,并找到质心的后验概率值。因为每个点就有一个高斯模型,因此会有多个高斯混合模型的质心需要处理,此时要对所有的高斯模型的质心的贡献度进行后验概率值归一化处理。最后一旦知道对应后验概率值后就可以求解相应变换参数了。

论文算法的核心是在于使高斯混合模型的质心做为一个整体来连贯运行,这就保留了点集的拓扑结构。具体数学公式推导可以阅读论文。

3、使用pycpd来实现一致性漂移算法的点云配准

原始论文的作者提供了matlab版本的代码实现,也有大佬用纯python复现了该算法,称为pycpd,具体项目可以点击原文链接。

安装过程很简单,用下面的命令安装即可。

pip install pycpd

作者也提供了使用案例,这里我就照搬使用,实现了刚性配准,仿射配准与非刚性配准。

from functools import partial
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from pycpd import deformable_registration
import numpy as np
import time


def visualize(iteration, error, X, Y, ax):
    plt.cla()
    ax.scatter(X[:,0],  X[:,1], X[:,2], color='red', label='Target')
    ax.scatter(Y[:,0],  Y[:,1], Y[:,2], color='blue', label='Source')
    ax.text2D(0.87, 0.92, 'Iteration: {:d}\nError: {:06.4f}'.format(iteration, error), horizontalalignment='center', verticalalignment='center', transform=ax.transAxes, fontsize='x-large')
    ax.legend(loc='upper left', fontsize='x-large')
    plt.draw()
    plt.pause(0.001)


def main():
    fish_target = np.loadtxt('data/fish_target.txt')
    X1 = np.zeros((fish_target.shape[0], fish_target.shape[1] + 1))
    X1[:,:-1] = fish_target
    X2 = np.ones((fish_target.shape[0], fish_target.shape[1] + 1))
    X2[:,:-1] = fish_target
    X = np.vstack((X1, X2))
    fish_source = np.loadtxt('data/fish_source.txt')
    Y1 = np.zeros((fish_source.shape[0], fish_source.shape[1] + 1))
    Y1[:,:-1] = fish_source
    Y2 = np.ones((fish_source.shape[0], fish_source.shape[1] + 1))
    Y2[:,:-1] = fish_source
    Y = np.vstack((Y1, Y2))

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    callback = partial(visualize, ax=ax)
    reg = deformable_registration(**{ 'X': X, 'Y': Y })
    reg.register(callback)
    plt.show()


if __name__ == '__main__':
    main()

4、配准结果

原始数据是人脸点云数据,通过对原始人脸点云数据进行变换(平移,旋转,缩放,扭曲等),如下图所示显示了原始人脸数据和变换后的人脸数据。

刚性变换配准迭代的最后结果如下所示。

仿射变换配准结果如下图所示。

非刚性变换配准结果如下图所示。

上面三种配准变换也可以结合起来使用,如下图是术前CT和术中DSA主动脉血管的中心线,CT的中心线是3d,DSA的中心线是2d的。首先采用刚性变换对3d中心线进行旋转,然后再采用仿射变换将3d中心线与2d中心线进行粗配准,最后采用非刚性变换将3d中心线与2d中心线进行精配准。

首先,对3d中心线进行刚性变换配准,结果如下图所示。

然后,在刚性变换配准的结果上进行仿射变换配准,结果如下图所示。

最后,在仿射变换配准结果上进行非刚性变换配准,结果如下图所示。

如果碰到任何问题,随时留言,我会尽量去回答的。

本文分享自微信公众号 - 最新医学影像技术(MedicalHealthNews),作者:最新医学影像技术

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-07-10

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • Luna16——肺结节检测和良恶性分类挑战赛(四)

    在luna16肺结节良恶性分类得例子中,有位细心的朋友提出一个很好的问题。今天首先分析上一篇中存在的问题,然后给出解决方案。

    用户7498388
  • Tensorflow和Pytorch深度学习框架安装教程

    目前主流深度学习框架有Tensorflow和pytorch,由于一些原因我只在windows10下安装了以上两个深度学习框架。Tensorflow在16年底就出...

    用户7498388
  • Hololens开发学习(五)——场景共享

    这一篇主要讲解Hololens场景共享。HoloTolkit5.8Test中提供了场景共享案例,已经实现了同步世界坐标系功能,我们只需要在此基础上进行代...

    用户7498388
  • python︱用asyncio、aiohttp实现异步及相关案例

    Asyncio 是并发(concurrency)的一种方式。对 Python 来说,并发还可以通过线程(threading)和多进程(multiprocessi...

    素质
  • 使用分布式缓存对ASP.Net Core性能提升?

    程序你好
  • LeetCode12|两个数组的交集

    这道题没有什么难度,就是一次集合的简单操作,不过每个人的解题思路都不一样,看自己的操作吧,今天写了这部分文章之后也觉得文章的风格没怎么变化,简单,给与示例程序,...

    后端Coder
  • 【程序源代码】开源学校教务管理系统

    后台:权限控制,支持多个管理员,学生管理,学生成绩,教师管理,文章管理,站点管理,网站布局自动化,多导航模式,友情链接,站点工具。

    程序源代码
  • Python + Flask程序调试技巧(2):导入自定义模块失败

    当导入字定义模块文件:model/rtindex中自定义的函数:pay_sum_query()时,提示:"unresolved reference"。具体现象如...

    数据饕餮
  • 07-Shell编程-数值运算符号使用

    [root@node1 ~]# IPADDR=$(ifconfig eth0 | grep "inet" | awk '{print $2}')

    泽阳
  • 映射真实世界后,百度地图开始用人工智能优化这个世界

    今天,百度地图与北京交管局有一个很有意思的战略合作,这将直接影响北京市民的出行:百度地图为北京交警量身打造了一个城市灯控路口路况监测平台——百度地图智慧信号灯研...

    罗超频道

扫码关注云+社区

领取腾讯云代金券