ITK库---GPU类梳理:

GPU类梳理:

ITK库---GPU类梳理:

1.核心GPU并行:

1GPU公共类:

  • 基于GPU的核管理器:itk::GPUKernelManager使用OpenCL实现的GPU内核管理器,此类负责管理GPU内核和命令队列。
  • 基于GPU的上下文管理器:itk::GPUContextManager
  • 基于GPU的数据管理器:itk::GPUDataManager 使用OpenCL实现GPU内存管理,是GPUImage类必须;
  • GPUFunctor图像滤波的基本Functor类:itk::GPUFunctorBase;其中的成员函数SetGPUKernelArguments()可以设置GPU核函数的参数;
  • itk::GPUImageDataManagerGPUImage的数据管理,将负责CPU映像和GPU映像之间的数据同步。
  • itk::GPUImage n-维图像模板化GPU执行类;GPU图像滤波一起使用,隐式管理CPU和GPU内存。
  • itk::GPUImageOps:提供一些基本的GPU图像操作的内核;
  • itk::GPUInPlaceImageFilter:GPU滤波器的基类;父类为InPlaceImageFilter;
  • itk::GPUImageToImageFilter:接受CPU和GPU图像作为输入和输出,并相应地使用滤波器。 如果可以使用GPU,则将调用GPUGenerateData()。 否则,将调用父类(即ImageToImageFilter)中的GenerateData()。
  • itk::GPUUnaryFunctorImageFilter:使用GPU在一张图像上实现逐像素操作,一元函子图像滤波器的GPU版本;
  • itk::GPUReduction:GPU版本的降维算法;
  • itk::OpenCLUtil

2GPU有限差分类:

  • itk::GPUFiniteDifferenceFunction:GPU有限差分函数基类,与大多数GPU类不同,GPUFiniteDifferenceFunction没有相应的CPU父类,而仅以CPU类FiniteDifferenceFunction作为其基类。
  • itk::GPUFiniteDifferenceImageFilter:GPU版本的有限差分滤波
  • itk::GPUDenseFiniteDifferenceImageFilter:GPU版密集有限差分图像过滤器

2、配准GPU并行:

1GPU公共类:

  • itk::GPUPDEDeformableRegistrationFunction执行形变配准算法的所有偏微分方程(PDE)函数的抽象基类, PDEDeformationRegistrationFilter子类使用它来计算输出变形字段,该字段将运动图像映射到固定图像。
  • itk::GPUPDEDeformableRegistrationFilter:GPUPDEDeformableRegistrationFilter算法通过计算变形场来两副图像配准,该变形场会将运动图像映射到固定图像上。

2GPUPDE形变类:

  • itk::GPUDemonsRegistrationFunction封装了Demons配准算法,DemonsRegistrationFilter使用它来计算输出变形字段,该字段将运动图像映射到固定图像。
  • itk::GPUDemonsRegistrationFilter:此类封装了demons配准PDE算法,DemonsRegistrationFilter使用它来计算输出变形场,将运动图像映射到固定图像。

3、滤波GPU并行:

1GPUAnisotropic平滑滤波:

  • itk::GPUScalarAnisotropicDiffusionFunction
  • itk::GPUAnisotropicDiffusionFunction:基于GPU的Anisotropic扩撒函数类,提供GPU版的有限差分函数的计算更新方法和计算平均梯度幅度方差;
  • GPUAnisotropic滤波tk::GPUAnisotropicDiffusionImageFilter类也是AnisotropicDiffusionImageFilter的GPU基类;可以通过InitializeIteration调用GPU计算平均梯度幅度平方函数(GPU Calculate Average Gradient Magnitude Squared())
  • itk::GPUGradientNDAnisotropicDiffusionFunction:为GPU上的标量值图像实现了经典的Perona-Malik各向异性扩散方程的N维版本。
  • itk::GPUGradientAnisotropicDiffusionImageFilter该类利用GPUGradientNDAnisotropicDiffusionFunction实现的经典Perona-Malik基于梯度幅度的方程式,对标量itk :: Image进行各向异性扩散。

2GPU图像滤波基类:

  • itk::GPUBoxImageFilterBox邻域中GOU滤波的基类,和itk::GPUBoxImageFilter类相对应。
  • itk::GPUNeighborhoodOperatorImageFilter利用GPU将单个NeighborhoodOperator用于图像滤波,通过GPU计算NeighborhoodOperator和NeighborhoodIterator之间的连续内积,该乘积扫过图像区域中的每个像素。
  • itk::GPUcastImageFilter:

3GPU平滑滤波:

  • Itk::GPUDiscreteGaussianImageFilter:通过GPU实现离散高斯图像滤波;
  • Itk::GPUMeanImageFilter:通过GPU实现均值滤波,全局内存中读取邻域像素。

4GPU阈值化:

  • itk::Functor::GPUBinaryThreshold:成员函数:SetGPUKernelArguments():设置该functor的GPU的核参数;SetInsideValue ()、SetOutsideValue ():设置内部外部参数;SetLowerThreshold ()、SetUpperThreshold ()设置上下限阈值;
  •  itk::GPUBinaryThresholdImageFilter:其中包含了GPU、CPU超类
  • itk::GPUBinaryThresholdImageFilterFactory类:是对GPUBinaryThresholdImageFilter类的实现;

4.CPU和GPU并行模块比较分析

公共类

(1) Image与GPUImage

1)头文件:增加

  • #include "itkGPUImageDataManager.h"
  • #include "itkVersion.h"
  • #include "itkObjectFactoryBase.h"

2)设置存取器类型

  • using AccessorFunctorType = DefaultPixelAccessorFunctor< Self >;
  • using NeighborhoodAccessorFunctorType = NeighborhoodAccessorFunctor< Self >;

3)分配CPU和GPU内存空间

4)显式同步CPU / GPU缓冲区

  • void UpdateBuffers();

5)返回Pixel Accessor对象(GPU与CPU)

AccessorType GetPixelAccessor()

{     m_DataManager->SetGPUBufferDirty();

      return Superclass::GetPixelAccessor();}

6)返回一个指向容器的指针。

PixelContainer * GetPixelContainer()

{     m_DataManager->SetGPUBufferDirty(); return Superclass::GetPixelContainer();}

7)将数据和信息从一个GPUImage移植到另一个GPUImage。

virtual void Graft(const Self *data);

8)添加一个工厂实现类

 

(2ImageToImageFilterGPUImageToImageFilter

1) 继承:

  • GPUImageToImageFilter继承ImageToImageFilter类

2) 宏设置是否使用GPU

  • itkSetMacro(GPUEnabled, bool);
  • itkGetConstMacro(GPUEnabled, bool);
  • itkBooleanMacro(GPUEnabled);

3)数据生成:

  • void GenerateData() override;

4)GPU核管理:#include "itkGPUKernelManager.h"

  • typename GPUKernelManager::Pointer m_GPUKernelManager;

(3)UnaryFunctorImageFilter

1头文件:

  • 继承UnaryFunctorImageFilter并包含:
  • #include "itkGPUFunctorBase.h"
  • #include "itkGPUInPlaceImageFilter.h"

2) 类模板定义:

  • template< typename TInputImage, typename TOutputImage, typename TFunction >
  • template< typename TInputImage, typename TOutputImage, typename TFunction,
  • typename TParentImageFilter = InPlaceImageFilter< TInputImage, TOutputImage > >
  • void GenerateOutputInformation() override;
  • void GPUGenerateData() override;
  • int m_UnaryFunctorImageFilterGPUKernelHandle;

(4)GPUFiniteDifferenceFunction

1)头文件:

  • 继承FiniteDifferenceFunction;
  • #include "itkGPUDataManager.h"
  • #include "itkGPUKernelManager.h"

2)GPU函数计算更新缓冲区

  • virtual void GPUComputeUpdate( const typename TImageType::Pointer output, ypename TImageType::Pointer update, void *gd) = 0;

3)分配GPU缓冲区以计算度量统计量

  • virtual void GPUAllocateMetricData(unsigned int itkNotUsed(numPixels)) {}

4)释放GPU缓冲区以计算度量统计量

  • virtual void GPUReleaseMetricData() {}

5)GPUFiniteDifferenceFunction类的GPU内核管理器

  • typename GPUKernelManager::Pointer m_GPUKernelManager;

6)GPUComputeUpdate()的GPU内核句柄*

  • int m_ComputeUpdateGPUKernelHandle;

(5) GPUFiniteDifferenceFilter

1)类模板定义

2)由子类定义更改应用于更新缓冲区和时间步长值“dt”的输出。

  • virtual void GPUApplyUpdate(const TimeStepType& dt) = 0;

3)由子类定义使用输出中像素的变化填充更新缓冲区,它返回用于更新的时间步长值

  • virtual TimeStepType GPUCalculateChange() = 0;

4)计算有限差分解决方案的默认高级算法

  • void GPUGenerateData() override;

图像滤波

(1)AnisotropicSmoothing

内部函数

  • itk::ScalarAnisotropicDiffusionFunction
  • itk::GPUScalarAnisotropicDiffusionFunction

1)为GPUScalarAnisotropicDiffusionFunction创建一个辅助GPU内核类

  • itkGPUKernelClassMacro(GPUScalarAnisotropicDiffusionFunctionKernel)

1)GPU数据指针类型

  • using GPUDataPointer = GPUDataManager::Pointer;

2)将OpenCL内核源获取为字符串,创建一个GetOpenCLSource方法

  • itkGetOpenCLSourceFromKernelMacro(GPUScalarAnisotropicDiffusionFunctionKernel)

3)计算平均梯度幅度平方

  • void CalculateAverageGradientMagnitudeSquared(TImage *) override;
  • void GPUCalculateAverageGradientMagnitudeSquared(TImage *) override;

滤波器

  • itkGradientAnisotropicDiffusionImageFilter 和 itkGPUGradientAnisotropicDiffusionImageFilter

1)头函数:

  • #include "itkOpenCLUtil.h"
  • 其次,包含GPU版本的对应头文件;包含非GPU的滤波类头文件

2)类模板定义:

  • CPU:template< typename TFixedImage, typename TMovingImage, typename TDisplacementField >
  • GPU:template< typename TFixedImage, typename TMovingImage, typename TDisplacementField, typename TParentImageFilter = itk::GradientAnisotropicDiffusionImageFilter<TFixedImage, TMovingImage, TDisplacementField >

3)ImageToImageFilter类定义输入输出均为图像的超类;

  • using Superclass = ImageToImageFilter< TInputImage, TOutputImage >;
  • using GPUSuperclass = GPUImageToImageFilter< TInputImage, TOutputImage, TParentImageFilter >;

4)部分会加入Object Factory implemenatation

5)对GPU参数进行打印;

 

(2)GPUNeighborhoodOperatorImageFilter

1)除了上述不同以外,加入了以上两个内容;

2)设置GPU缓存区类型:

  • using NeighborhoodGPUBufferType = GPUImage< TOperatorValueType, Self::ImageDimension >;
  • 利用头文件#include "itkGPUImageToImageFilter.h"
  • void GPUGenerateData() override;

 

(3)DiscreteGaussianImageFilter

1)头文件:

  • #include "itkImageToImageFilter.h"
  • #include "itkImage.h"
  •  
  • #include "itkGPUImage.h"
  • #include "itkGPUImageToImageFilter.h"
  • #include "itkGPUNeighborhoodOperatorImageFilter.h"
  • #include "itkDiscreteGaussianImageFilter.h"

2)ImageToImageFilter类定义输入输出均为图像的超类;

3)设置滤波器类型

  • using IntermediateFilterType = GPUNeighborhoodOperatorImageFilter<RealOutputImageType,RealOutputImageType, RealOutputPixelValueType >;
  • using LastFilterType = GPUNeighborhoodOperatorImageFilter<RealOutputImageType,OutputImageType, RealOutputPixelValueType >;
  • using SingleFilterType = GPUNeighborhoodOperatorImageFilter<InputImageType,OutputImageType,RealOutputPixelValueType >;

(4)MeanImageFilter

1)多包括一个头文件

  • #include "itkBoxImageFilter.h"
  • #include "itkGPUBoxImageFilter.h"

2)创建一个辅助GPU内核类

  • itkGPUKernelClassMacro(GPUMeanImageFilterKernel);

3)将OpenCL内核源获取为字符串,创建一个GetOpenCLSource方法

  • itkGetOpenCLSourceFromKernelMacro(GPUMeanImageFilterKernel);

4)添加一个类工厂实现

(5)BinaryThresholdImageFilter

1)头文件:

  • #include "itkOpenCLUtil.h"
  • #include "itkGPUFunctorBase.h"
  • #include "itkGPUKernelManager.h"

2)设置GPU核参数:

  • int SetGPUKernelArguments(GPUKernelManager::Pointer KernelManager, int KernelHandle)
  • {
    • KernelManager->SetKernelArg(KernelHandle, 0, sizeof(TInput), &(m_LowerThreshold) );
    • KernelManager->SetKernelArg(KernelHandle, 1, sizeof(TInput), &(m_UpperThreshold) );
    • KernelManager->SetKernelArg(KernelHandle, 2, sizeof(TOutput), &(m_InsideValue) );
    • KernelManager->SetKernelArg(KernelHandle, 3, sizeof(TOutput), &(m_OutsideValue) );
    • return 4;
  • }

3)将OpenCL内核源获取为字符串,创建一个GetOpenCLSource方法

  • itkGetOpenCLSourceFromKernelMacro(GPUBinaryThresholdImageFilterKernel);

4)GPUBinaryThresholdImageFilter目标工厂实现

GPU图像配准模块

(1)GPUDeformableRegistrationFunction

1)头文件:

  • 添加相关GPU版本的类,
  • #include "itkGPUReduction.h"

2)创建一个辅助GPU内核类

  • itkGPUKernelClassMacro(GPUDemonsRegistrationFunctionKernel);

3)标准类别名定义:

4)定义GPU数据指针类型:

  • using GPUDataPointer = GPUDataManager::Pointer;

5)将OpenCL内核源获取为字符串,创建一个GetOpenCLSource方法

  • itkGetOpenCLSourceFromKernelMacro(GPUDemonsRegistrationFunctionKernel);

6)分配GPU缓冲区以计算度量统计信息

  • void GPUAllocateMetricData(unsigned int numPixels) override;

7)释放GPU缓冲区以计算度量统计信息

  • void GPUReleaseMetricData() override;

8)有限差分求解器图像过滤器在不位于数据集边界上的每个像素处调用此方法

  • virtual PixelType ComputeUpdate( const NeighborhoodType & neighborhood,void *globalData,const FloatOffsetType & offset = FloatOffsetType(0.0)) override;
  • virtual void GPUComputeUpdate( const DisplacementFieldTypePointer output, DisplacementFieldTypePointer update,void *gd) override;

9)设置GPUComputeUpdat的GPU内核句柄

  • int m_ComputeUpdateGPUKernelHandle;

10)度量:增加降维参数,计算出的固定图像和变换运动图像之间的强度均方差

  • mutable GPUReduction<int>::Pointer   m_GPUPixelCounter;
  • mutable GPUReduction<float>::Pointer m_GPUSquaredChange;
  • mutable GPUReduction<float>::Pointer m_GPUSquaredDifference;

(2)GPUDeformableRegistrationFilter

1)头文件

2)类模板定义

3)标准类别名定义:

4)设置函数类型:

  • using GPUDemonsRegistrationFunctionType = GPUDemonsRegistrationFunction<FixedImageType,MovingImageType, DisplacementFieldType >;

5)GPUDemonsRegistrationFilter的简要目标工厂实现:

  • class GPUDemonsRegistrationFilterFactory : public itk::ObjectFactoryBase{……}

(3)GPUPDEDeformableRegistrationFunction

  • 将所有相关类定义置换为GPUPDEDeformableRegistrationFunction即可;

(4)GPUPDEDeformableRegistrationFilter

1)创建一个辅助GPU内核类

  • itkGPUKernelClassMacro(GPUDemonsRegistrationFunctionKernel);

2)类模板定义:

  • template< typename TFixedImage, typename TMovingImage, typename TDisplacementField, typename TParentImageFilter = PDEDeformableRegistrationFilter< TFixedImage, TMovingImage, TDisplacementField >>

3)标准类别名定义:

4)设置函数类型:

  • using GPUPDEDeformableRegistrationFunctionType = GPUPDEDeformableRegistrationFunction<FixedImageType, MovingImageType,DisplacementFieldType >;

5)将OpenCL内核源获取为字符串,创建一个GetOpenCLSource方法

  • itkGetOpenCLSourceFromKernelMacro(GPUPDEDeformableRegistrationFilterKernel);

6)平滑位移场内核的内存缓冲区

  • typename GPUDataManager::Pointer m_GPUSmoothingKernels[ImageDimension];

7)设置GPUComputeUpdat的GPU内核句柄

  • int m_SmoothDisplacementFieldGPUKernelHandle;

GPU图像配准代码差异

itkGPUDemonsRegistrationFilterTest.cxx

1)包含有效的头文件:

  • #include "itkDemonsRegistrationFilter.h"
  • #include "itkGPUDemonsRegistrationFilter.h"
  • #include "itkGPUImage.h"
  • #include "itkGPUKernelManager.h"
  • #include "itkGPUContextManager.h"

2)计时器:

  • itk::TimeProbe类用于计时;
  • itk::TimeProbe m_GPUTime;
  • m_GPUTime.GetMean();
  • m_GPUTime.Start();
  • filter->Update();
  • m_GPUTime.Stop();

3)图像模板化:

  • using CPUDisplacementFieldType = itk::Image<VectorPixelType, ImageDimension>;
  • using GPUDisplacementFieldType = itk::GPUImage<VectorPixelType, ImageDimension>;

4)配准滤波器

  • using CPURegistrationFilterType = itk::DemonsRegistrationFilter<ImageType,ImageType,FieldType>;
  • using GPURegistrationFilterType = itk::GPUDemonsRegistrationFilter<InternalImageType,InternalImageType, DisplacementFieldType>;

GPU模块总结

1)头文件添加:

  • 继承CPU类头文件;
  • 添加GPU管理器(上下文、数据、核、图像数据等);
  • 添加itkOpenCLUtil.h;
  • 添加相关GPU版本的图像模板化类、滤波器类、滤波器函数类;
  • 可以添加计时器等类;

2)需要创建辅助GPU内核类:itkGPUKernelClassMacro(类名Kernel);

3)修改类模板定义格式:

4)更新定义GPU版本的标准类别名:通过具体的GPU模板的类;

5)利用宏定义设置是否使用GPU;

6)设置GPU核参数:SetGPUKernelArguments();

7)内存分配:CPU核GPU;

8)返回一个指向容器的指针;

9)GPU函数计算更新缓冲区GPUComputeUpdate ();

10)设置GPU缓存区类型;

11)分配与释放GPU缓冲区、计算度量统计信息量;

12)某类的GPU内核管理器、内核句柄定义;

13)利用GPU版本的评价指标计算函数执行计算;

14)GPU数据指针类型定义;

15)对GPU参数进行打印;

16)相关类的对象工厂实现定义;

5、ITK库GPU并行模块的OpenCL实现

  • ITK库主流是利用C++语言实现,并且进行封装(Wrapping);
  • GPU模块是通过开放运算语言OpenCL实现,它是第一个面向异构系统的并行编程的开放式标准;
  • CUDA核与OpenCL
    • Block:相当于 OpenCL中的work-group
    • Thread:相当于OpenCL中的work-item
    • SP:相当于OpenCL中的PE
    • SM:相当于OpenCL中的CU
    • Warp:相当于OpenCL中的wave front(简称wave),基本的调试单位
  • 计算索引,利用多线程替换CPU程序中的循环机构的计算,每个线程的index与要计算的向量的index相同;
  • CUDA写在程序内,OpenCL写在一个独立文件,并且文件后缀是.cl,由主机代码读入后执行;
  • 多个线程执行同样的核函数,每个Work-item都有一个唯一固定的ID号,一般通过这个ID号来区分需要处理的数据。

(1)GPU版的Demons配准函数的OpenCL实现

1)GPUDemonsRegistrationFunction:

  • CMakeList中设置:if (ITK_USE_GPU)endif()
  • OpenCL实现:(src\*.cl),通常是利用OpenCL实现1D\2D\3D操作;
  • 利用__kernel void ComputeUpdate ()来更新,同时通过宏定义依据输入参数决定调用那个维度的计算。
  • 其中主要涉及:1D\2D\3D的线性插值计算,1D\2D\3D的梯度计算;
  • GPUPDEDeformableRegistrationFilter:
  • 计算全局内存索引:1D设置ysize=1;2D设置zsize=1;3D
  • 平滑滤波器的Reorder
  • OpenCL实现:数据迁移从host->global->local……
    • __kernel void SmoothingFilterReorder( __global OUTPIXELTYPE *imgIn,
    • __global OUTPIXELTYPE *imgOut,
    • __constant int *imgSize, int dim,
    • __constant OUTPIXELTYPE *filter, int filterSize, int indir, int outdir,
    • __local volatile OUTPIXELTYPE *sharedFilter,
    • __local volatile OUTPIXELTYPE *sharedData)
  • Kernel关键字定义了kernel函数,返回值必须是void;
  • Host关键字定义主机内存(CPU);        Global关键字定义全局内存;
  • Constant关键字定义常量内存;              Local关键字定义局部内存;
  • Private关键字定义私有内存(每个工作项);
  • 设置参数:blocksize、filtersize、dim、id、indir等
  • 在indir维度方向加载元素,完成后释放内存;
  • 并行计算内积

(2)ImageOpsOpenCL实现

  • 实现图像的四则运算:加减乘除

(3)有限差分滤波

  • GPUDenseFiniteDifferentImageFilter
  • 必须传入图像像素维度(1D\2D\3D)
  • 正向差分:

ITK库---GPU类梳理:

  • 反向差分:

ITK库---GPU类梳理:

  • 中心差分:

ITK库---GPU类梳理:

  • OpenCL 1.0版本在内核执行期间,可能超过三级的嵌套条件语句的命令队列会无效情况。
    • bool isValid = true;
    • if(gix < 0 || gix >= width) isValid = false;
    • if(giy < 0 || giy >= height) isValid = false;
    • if(giz < 0 || giz >= depth) isValid = false;
    • for(int i=0; i<PIXELDIM; i++)
    • {
    • out[gidx] = (OUTPIXELTYPE)( (float)out[gidx] + dt*(float)(buf[gidx]) );
    • gidx ++;
    • }

(4)二进制阈值滤波

  • BinaryThresholdImageFilter:实现了1D\2D\3D操作;
  • 得到x,y,z,三个维度上的索引id;
  • unsigned int gidx = width*(giz*height + giy) + gix;
  • out[gidx] = Functor(lowerThreshold, upperThreshold, insideValue, outsideValue, in[gidx]);