CN106384350B - 基于cuda加速的神经元活动图像动态配准方法及装置 - Google Patents

基于cuda加速的神经元活动图像动态配准方法及装置 Download PDF

Info

Publication number
CN106384350B
CN106384350B CN201610861004.8A CN201610861004A CN106384350B CN 106384350 B CN106384350 B CN 106384350B CN 201610861004 A CN201610861004 A CN 201610861004A CN 106384350 B CN106384350 B CN 106384350B
Authority
CN
China
Prior art keywords
image
translational movement
rigid transformation
pixel
cuda
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.)
Active
Application number
CN201610861004.8A
Other languages
English (en)
Other versions
CN106384350A (zh
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.)
Institute of Automation of Chinese Academy of Science
Original Assignee
Institute of Automation of Chinese Academy of Science
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 Institute of Automation of Chinese Academy of Science filed Critical Institute of Automation of Chinese Academy of Science
Priority to CN201610861004.8A priority Critical patent/CN106384350B/zh
Publication of CN106384350A publication Critical patent/CN106384350A/zh
Application granted granted Critical
Publication of CN106384350B publication Critical patent/CN106384350B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10056Microscopic image
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20048Transform domain processing
    • G06T2207/20056Discrete and fast Fourier transform, [DFT, FFT]

Landscapes

  • Image Analysis (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种基于CUDA加速的神经元活动图像动态配准方法,包括:对输入图像做均值拉伸处理,得到均值拉伸处理后的图像;利用傅里叶变换计算均值拉伸处理后的图像整体的平移量;利用所述平移量和模板图像对均值拉伸处理后的图像做刚性变换,得到第一次刚性变换后的图像;将第一次刚性变换后的图像和模板图像分成若干个小块区域,利用傅里叶变换计算每个小块区域的平移量;利用所述小块区域的平移量计算第一次刚性变换后的图像中每个像素点的平移量,并利用平移量和模板图像更新所述第一次刚性变换后的图像中的每个像素点,得到配准图像。本发明基于CUDA架构实现神经元活动图像配准的并行加速计算,实现神经元活动图像的实时配准。

Description

基于CUDA加速的神经元活动图像动态配准方法及装置
技术领域
本发明涉及图像配准技术,尤其涉及一种基于CUDA加速的神经元活动图像动态配准方法及装置。
背景技术
随着光片显微成像技术的发展,目前可以快速采集单细胞分辨率的活体斑马鱼神经图像。由于斑马鱼受刺激或者自发的运动,导致采集的神经元活动图像出现一定的形变。由此,对于斑马鱼神经元活动图像需要进行配准处理,以对齐神经元活动图像中的神经元,便于后续神经元活动的统计分析和进一步处理。
目前常用的神经元活动图像配准方法包括刚性变换和非刚性变换,刚性变换只能实现图像的平移,旋转,缩放变换,不能解决弹性形变,非刚性变换可以解决弹性形变,但是算法较复杂,不能实现实时配准。为了准确的分析神经元的活动情况,特别是在接受某种刺激时的活动情况,需要研究快速的神经元活动图像配准方法。考虑到采集的斑马鱼图像只是在局部区域发生了弹性形变,所以可以将图像分成若干个重叠的区域实施刚性变换,从而近似实现弹性变换。考虑到GPU的高效计算能力,构建了基于CUDA加速的神经元活动图像配准框架,像素点的处理通过GPU的线程实现并行计算。对于高分辨率的神经元活动图像,单个GPU的平均配准时间达到了毫秒级别,通过扩展多个GPU处理,可以进一步降低配准时间,从而实现实时配准。
发明内容
基于上述技术问题,本发明的目的在于提供一种基于CUDA加速的神经元活动图像动态配准方法,从而对神经元活动图像中的神经元进行位置动态配准,便于后续神经功能回路机理的解析。
为了实现上述目的,根据本发明的一个方面,本发明提供了一种基于CUDA加速的神经元活动图像动态配准方法,包括:
步骤S1:对输入图像做均值拉伸处理,得到均值拉伸处理后的图像;
步骤S2:利用傅里叶变换计算均值拉伸处理后的图像整体的平移量(Δx,Δy);
步骤S3:利用所述平移量(Δx,Δy)和模板图像对均值拉伸处理后的图像做刚性变换,得到第一次刚性变换后的图像;
步骤S4:将第一次刚性变换后的图像和模板图像分成若干个小块区域,利用傅里叶变换计算每个小块区域的平移量(Δxi,Δyi);
步骤S5:利用所述小块区域的平移量(Δxi,Δyi)计算第一次刚性变换后的图像中每个像素点的平移量,并利用平移量和模板图像更新所述第一次刚性变换后的图像中的每个像素点,得到配准图像。
根据本发明第二方面,提供了一种基于CUDA加速的神经元活动图像动态配准装置,其特征在于,包括:
均值拉伸模块,用于对输入图像做均值拉伸处理,得到均值拉伸处理后的图像;
整体平移量计算模块,用于利用傅里叶变换计算均值拉伸处理后的图像整体的平移量(Δx,Δy);
刚性变换模块,用于利用所述平移量(Δx,Δy)和模板图像对均值拉伸处理后的图像做刚性变换,得到第一次刚性变换后的图像;
块平移量计算模块,用于将第一次刚性变换后的图像和模板图像分成若干个小块区域,利用傅里叶变换计算每个小块区域的平移量(Δxi,Δyi);
配准模块,用于利用所述小块区域的平移量(Δxi,Δyi)计算第一次刚性变换后的图像中每个像素点的平移量,并利用平移量和模板图像更新所述第一次刚性变换后的图像中的每个像素点,得到配准图像。
基于上述技术方案可知,本发明的神经元活动图像配准方法利用傅里叶变换的性质计算待配准图像的整体平移量,利用平移量对待配准图像做刚性变换,然后将待配准图像分成若干个互相重叠的小区域,利用傅里叶变换的性质计算每个小区域的平移量,进而得到每个像素点的平移量,利用平移量对待配准图像做第二次变换,得到输出的配准图像。实验证明该方法可以快速的实现神经元活动图像的配准,有效对齐神经元活动图像中的神经元,便于后续神经元活动的统计分析和进一步处理。本发明从刚性变换着手,构建近似弹性变换,区别于常规的刚性变换和弹性变换方法,能够快速的实现神经元活动图像的配准,对神经元的配准效果较好。
附图说明
图1为作为本发明一实施例的基于CUDA加速的神经元活动图像动态配准方法的流程图;
图2为作为本发明一实施例的中心变换示意图;
图3为作为本发明一实施例的S4步骤中的图像分块示意图。
图4为作为本发明一实施例的基于CUDA加速的神经元活动图像动态配准批处理流程图。
图5为作为本发明一实施例的基于CUDA加速的神经元活动图像配准的CUDA架构。
图6为作为本发明一实施例的基于CUDA加速的神经元活动图像配准的运行时间结果示意图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚明白,以下结合具体实施例,并参照附图,对本发明作进一步的详细说明。
下面结合附图对本发明的基于CUDA(Compute Unified Device Architecture,统一计算设备架构)加速的神经元活动图像配准方法进行详细描述。对表示符号做如下声明:
模板图像:uref,输入图像:u0,均值拉伸处理后的图像:u1,第一次刚性变换后的图像:u2,输出配准图像:u3
待配准图像的长度为M,宽度为N,图像中(i,j)位置的像素值为v(i,j),图像的灰度均值为m,待配准图像整体的平移量为(Δx,Δy),待配准图像第i个小区域的平移量为(Δxi,Δyi)。
图1为作为本发明一个实施例的基于CUDA加速的神经元活动图像配准方法的流程图。
参照图1,在步骤S101,对输入图像u0做均值拉伸处理,得到均值拉伸处理后的图像u1(t)。
首先计算输入图像u0的灰度均值,其中,m为输入图像u0的灰度均值,M为输入图像的长度,N为输入图像的宽度,i、j分别表示图像中的像素位置,v(i,j)表示输入图像中像素(i,j)位置的像素值;然后对输入图像中的每一像素点做均值拉伸处理,v′(i,j)=v(i,j)-m(i=1,2..M,j=1,2..N)。其中v(i,j)为输入图像的像素值,v′(i,j)为均值拉伸处理后的像素值。在计算输入图像的灰度均值时,首先计算输入图像的像素值之和,该过程中利用CUDA归约核函数实现,时间复杂度为○(log2MN),均值拉伸利用CUDA减法核函数实现并行操作。
在步骤S102,利用傅里叶变换的性质计算均值拉伸处理后的图像u1整体的平移量(Δx,Δy),其中,计算均值拉伸处理后的图像u1整体的平移量(Δx,Δy)包括:
计算均值拉伸处理后的图像u1与模板图像uref之间的平移量(Δx,Δy)。
计算过程根据下式进行:
F1(i,j)=fft(uref(i,j)) (1);
F2(i,j)=fft(u1(i,j)) (2);
F(i,j)=F2(i,j)*F1 *(i,j)/(F1(i,j)*F1 *(i,j)) (3);
f(i,j)=ifft(F(i,j)) (4);
f′(i,j)=fftshift(f(i,j)) (5);
(x0,y0)=argmax(f′(i,j)) (6);
其中,fft()为快速傅里叶变换函数,F1(i,j)为对模板图像uref进行快速傅里叶变换后的得到的结果;F2(i,j)为对均值拉伸处理后的图像u1进行快速傅里叶变换后的得到的结果;F1 *(i,j)为F1(i,j)的共轭,f(i,j)为对F(i,j)进行快速傅里叶变换得到的结果;ifft()为快速反傅里叶变换函数。f′(i,j)为对f(i,j)进行中心变换后得到的结果,fftshift()为中心变换函数,变换过程可参照图2,其中每个方格代表一个像素,其中的数值表示像素值,中心变换是指以图像中心为原点将图像划分为四个大小相同的区域,交换左上角区域和右下角区域的像素值,交换右上角区域和左下角区域的像素值。(x0,y0)取f′(i,j)最大值处的坐标值,最后需要对(x0,y0)做矫正处理,即公式(7)中的处理,最终得到整体的平移量(Δx,Δy)。fft()和ifft()函数利用CUDA的cufft库中的函数实现,fftshift()函数分情况设计了多个核函数实现变换处理,向量加法,乘法分别利用CUDA加法,乘法核函数实现,考虑到公式(6)中利用的仅仅是最大值所在的位置,与具体的数值无关,因此忽略了公式(3)中的模长计算,减少了计算量,提高了计算速度。
在步骤S103,利用平移量(Δx,Δy)和模板图像uref对均值拉伸处理后的图像u1做刚性变换,得到第一次刚性变换后的图像u2,刚性变换过程包括:对于均值拉伸处理后的图像u1中的每个像素点,利用下式分情况处理:
上述公式(8)中,CUDA架构中设计了像素值变换核函数,首先判断平移量类别,然后计算对应像素点的位置,利用对应像素点的位置更新该位置的像素值,核函数应用于图像中的每一个像素,实现并行处理。
在步骤S104,将第一次刚性变换后的图像u2和模板图像uref分成若干个小块区域pi,利用傅里叶变换的性质计算每个小块区域pi的平移量(xi,yi),图像分块情况可参照图3,图3将图像划分成了5×5个方格块,其中斜线区域为不同块之间相互重叠的区域,黑色圆点表示每个网格块的中心点,图中标明了点密度和方格宽度参数的指示情况,通过调节这两个参数可以更改图像分块情况,通过比较不同参数设置下的图像配准效果确定最佳的图像分块情况,神经元活动图像u2和模板图像uref的分块情况相同。以小块区域pi为单位执行步骤S102中的操作,计算神经元活动图像u2中每个小块区域相对模板图像uref中对应方格区域的平移量(Δxi,Δyi)。若Δxi或Δyi的值大于设定的最大平移量,则将其置为0。
在步骤S105,利用平移量(Δxi,Δyi)计算每个像素点的平移量,具体操作为遍历图像u2中的每个像素点,比较其所属方格区域的平移量,选取最大的平移量作为其平移量。针对图像u2中的每个像素点,将其坐标值减去平移量得到模板图像中对应像素点的坐标,利用模板图像中对应点的像素值替换图像u2中像素点的像素值,如果像素点在模板中的对应点出界了,则保持其像素值不变,从而得到结果图像u3
图5为作为本发明一实施例的基于神经元活动图像配准的CUDA架构。以上所述操作均在GPU上运行,最大程度的优化线程并行处理,在运行时间上达到毫秒级别,可以实现实时配准,该CUDA架构分别测试了100张,500张,1000张斑马鱼神经活动图像配准时间,并与串行处理的matlab程序运行时间做了对比,图6给出了相应的测试结果。
以上所述的具体实施例,对本发明的目的、技术方案和有益效果进行了进一步详细说明,应理解的是,以上所述仅为本发明的具体实施例而已,并不用于限制本发明,凡在本发明的精神和原则之内,所做的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于CUDA加速的神经元活动图像动态配准方法,其特征在于,包括:
步骤S1:基于CUDA减法核函数,对输入图像做均值拉伸处理进行并行操作,得到均值拉伸处理后的图像;
步骤S2:基于CUDA的cufft库中的函数、CUDA加法、乘法核函数,利用傅里叶变换计算均值拉伸处理后的图像整体的平移量(Δx,Δy);
步骤S3:在CUDA中设计像素值变换核函数,利用所述平移量(Δx,Δy)和模板图像对均值拉伸处理后的图像做刚性变换,得到第一次刚性变换后的图像;
步骤S4:将第一次刚性变换后的图像和模板图像分成若干个小块区域,利用傅里叶变换计算每个小块区域的平移量(xi,yi),再计算第一次刚性变换后的图像中每个小块区域相对模板图像中对应方格区域的平移量(Δxi,Δyi);
步骤S5:利用所述第一次刚性变换后的图像中每个小块区域相对模板图像中对应方格区域的平移量(Δxi,Δyi)计算第一次刚性变换后的图像中每个像素点的平移量,并利用平移量和模板图像更新所述第一次刚性变换后的图像中的每个像素点,得到配准图像。
2.如权利要求1所述的方法,其特征在于,步骤S1中对输入图像做均值拉伸处理包括:
基于CUDA归约核函数,计算输入图像的灰度均值其中,m为输入图像u0的灰度均值,M为输入图像的长度,N为输入图像的宽度,i、j分别表示输入图像中的像素位置,v(i,j)表示输入图像中像素点(i,j)位置处的像素值;
对输入图像中的每一像素点做均值拉伸处理,v′(i,j)=v(i,j)-m(i=1,2..M,j=1,2..N);其中v(i,j)为输入图像的像素值,v′(i,j)为均值拉伸处理后的像素值。
3.如权利要求1所述的方法,其特征在于,步骤S2中均值拉伸处理后的图像整体的平移量(Δx,Δy)如下计算:
F1(i,j)=fft(uref(i,j)) (1);
F2(i,j)=fft(u1(i,j)) (2);
f(i,j)=ifft(F(i,j)) (4);
f′(i,j)=fftshift(f(i,j)) (5);
(x0,y0)=argmax(f′(i,j)) (6);
其中,fft()为快速傅里叶变换函数,F1(i,j)为对模板图像uref进行快速傅里叶变换后的得到的结果;F2(i,j)为对均值拉伸处理后的图像u1进行快速傅里叶变换后的得到的结果;F(i,j)为F2(i,j)除以F1(i,j)的结果,为F1(i,j)的共轭,f(i,j)为对F(i,j)进行快速傅里叶变换得到的结果;ifft()为快速反傅里叶变换函数;f′(i,j)为对f(i,j)进行中心变换后得到的结果,fftshift()为中心变换函数;(x0,y0)取f′(i,j)最大值处的坐标值。
4.如权利要求1所述的方法,其特征在于,步骤S3中利用下式对均值拉伸处理后的图像做刚性变换:
其中,u2为第一次刚性变换后的图像;uref为模板图像;M为输入图像的长度,N为输入图像的宽度。
5.如权利要求1所述的方法,其特征在于,步骤S4所述第一次刚性变换后的图像中每个小块区域相对模板图像中对应方格区域的平移量(Δxi,Δyi)如下计算:
以所述小块区域为单位,执行步骤S2中的操作,计算第一次刚性变换后的图像中每个小块区域相对模板图像uref中对应方格区域的平移量(Δxi,Δyi);若Δxi或Δyi的值大于设定的最大平移量,则将其置为0。
6.如权利要求1所述的方法,其特征在于,步骤S5中第一次刚性变换后的图像中每个像素点的平移量如下计算:
遍历第一次刚性变换后的图像中的每个像素点,比较其所属小块区域的平移量,选取最大的平移量作为其平移量;
针对第一次刚性变换后的图像中的每个像素点,将其坐标值减去其平移量得到模板图像中对应像素点的坐标,利用模板图像中对应点的像素值更新第一次刚性变换后的图像中像素点的像素值,如果像素点在模板中的对应点出界了,则保持其像素值不变,从而得到配准图像。
7.一种基于CUDA加速的神经元活动图像动态配准装置,其特征在于,包括:
均值拉伸模块,基于CUDA减法核函数,对输入图像做均值拉伸处理,得到均值拉伸处理后的图像;
整体平移量计算模块,基于CUDA的cufft库中的函数、CUDA加法和乘法核函数,利用傅里叶变换计算均值拉伸处理后的图像整体的平移量(Δx,Δy);
刚性变换模块,在CUDA中设计像素值变换核函数,利用所述平移量(Δx,Δy)和模板图像对均值拉伸处理后的图像做刚性变换,得到第一次刚性变换后的图像;
块平移量计算模块,用于将第一次刚性变换后的图像和模板图像分成若干个小块区域,利用傅里叶变换计算每个小块区域的平移量(xi,yi),再计算第一次刚性变换后的图像中每个小块区域相对模板图像中对应方格区域的平移量(Δxi,Δyi);
配准模块,用于利用所述第一次刚性变换后的图像中每个小块区域相对模板图像中对应方格区域的平移量(Δxi,Δyi)计算第一次刚性变换后的图像中每个像素点的平移量,并利用平移量和模板图像更新所述第一次刚性变换后的图像中的每个像素点,得到配准图像。
CN201610861004.8A 2016-09-28 2016-09-28 基于cuda加速的神经元活动图像动态配准方法及装置 Active CN106384350B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610861004.8A CN106384350B (zh) 2016-09-28 2016-09-28 基于cuda加速的神经元活动图像动态配准方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610861004.8A CN106384350B (zh) 2016-09-28 2016-09-28 基于cuda加速的神经元活动图像动态配准方法及装置

Publications (2)

Publication Number Publication Date
CN106384350A CN106384350A (zh) 2017-02-08
CN106384350B true CN106384350B (zh) 2019-06-07

Family

ID=57937051

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610861004.8A Active CN106384350B (zh) 2016-09-28 2016-09-28 基于cuda加速的神经元活动图像动态配准方法及装置

Country Status (1)

Country Link
CN (1) CN106384350B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109540946A (zh) * 2018-10-25 2019-03-29 华东师范大学 基于人工智能算法的透射电子显微镜样品漂移矫正方法
CN111539997B (zh) * 2020-04-23 2022-06-10 中国科学院自动化研究所 基于gpu计算平台的图像并行配准方法、系统、装置
CN111709979B (zh) * 2020-05-15 2023-08-25 北京百度网讯科技有限公司 图像对齐的方法、装置、电子设备和存储介质

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6373970B1 (en) * 1998-12-29 2002-04-16 General Electric Company Image registration using fourier phase matching
CN103310458A (zh) * 2013-06-19 2013-09-18 北京理工大学 结合凸包匹配和多尺度分级策略的医学图像弹性配准方法
CN103761750A (zh) * 2014-02-14 2014-04-30 华中科技大学 一种心肌质点运动图像与心肌纤维走向图像配准方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2015009877A1 (en) * 2013-07-17 2015-01-22 The Regents Of The University Of California Method for focused recording and stimulation electrode array

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6373970B1 (en) * 1998-12-29 2002-04-16 General Electric Company Image registration using fourier phase matching
CN103310458A (zh) * 2013-06-19 2013-09-18 北京理工大学 结合凸包匹配和多尺度分级策略的医学图像弹性配准方法
CN103761750A (zh) * 2014-02-14 2014-04-30 华中科技大学 一种心肌质点运动图像与心肌纤维走向图像配准方法

Also Published As

Publication number Publication date
CN106384350A (zh) 2017-02-08

Similar Documents

Publication Publication Date Title
Kobler et al. Variational networks: connecting variational methods and deep learning
CN104680491B (zh) 一种基于深度神经网络的图像非均匀运动模糊去除方法
CN106384350B (zh) 基于cuda加速的神经元活动图像动态配准方法及装置
CN102073982B (zh) 用gpu实现超大sar图像各向异性扩散滤波加速方法
CN109191476A (zh) 基于U-net网络结构的生物医学图像自动分割新方法
Ugolotti et al. Particle swarm optimization and differential evolution for model-based object detection
Wu et al. Dynamic filtering with large sampling field for convnets
CN104156994A (zh) 一种压缩感知磁共振成像的重建方法
CN110781976B (zh) 训练图像的扩充方法、训练方法及相关装置
Li et al. Carn: Convolutional anchored regression network for fast and accurate single image super-resolution
RU2013134636A (ru) Способ, компьютерная программа и устройство для определения места захвата
CN110363760A (zh) 用于识别医学图像的计算机系统
CN104392409A (zh) 一种图像美容的加速方法
CN105809173A (zh) 一种基于仿生物视觉变换的图像rstn不变属性特征提取及识别方法
CN104751485A (zh) 一种基于gpu自适应的前景提取方法
Valle et al. Cascade of encoder-decoder CNNs with learned coordinates regressor for robust facial landmarks detection
CN102592135A (zh) 融合目标空间分布和时序分布特征子空间的视觉跟踪方法
An et al. Patch loss: A generic multi-scale perceptual loss for single image super-resolution
CN103839233A (zh) 一种相机抖动造成的模糊图像复原方法
CN108960203A (zh) 一种基于fpga异构计算的车辆检测方法
Menchón-Lara et al. Efficient convolution-based pairwise elastic image registration on three multimodal similarity metrics
Choe et al. Improved techniques for weakly-supervised object localization
CN111445503B (zh) 基于gpu集群上并行编程模型的金字塔互信息图像配准方法
CN106778831A (zh) 基于高斯混合模型的刚体目标在线特征分类与跟踪方法
Peppert et al. On the sufficient condition for solving the gap-filling problem using deep convolutional neural networks

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant