首页 > itarticle > applying gpu in demons registration by using itk library Read image file GPU registration Output

applying gpu in demons registration by using itk library Read image file GPU registration Output

admin 11月 13, 2020 0

I currently process CBCT image, mainly focusing on reconstruction. I use ITK as motion estimation tool. And to accelerate the calculation procedure, I also use GPU part in ITK

Read image file

My files are stored as raw binary format.Since ITK does not provide the raw bin file reader, I have to convert them to itk::Image type by my own. Luckily, ITK provides ImportImageFilter class to import data from a standard C array into an itk::Image.

image file

For image file, each voxel has one float number representing attenuation coefficient.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
 
typedef float   PixelType;
const unsigned int Dimension = 3;
typedef itk::Vector< PixelType, Dimension > VectorPixelType;
unsigned int nump = width * height * slices;
PixelType *FixedImageArray = new PixelType[nump];
std::ifstream (argv[1], std::ios::binary | std::ios::in);
FixedImageFile.read((char*)FixedImageArray, nump * sizeof(float));
FixedImageFile.close();
typedef itk::ImportImageFilter<PixelType, Dimension> ImageImportFilterType;    
ImageImportFilterType::Pointer FixedImportFilter = ImageImportFilterType::New();    
ImageImportFilterType::IndexType start;
start[0] = 0;
start[1] = 0;
start[2] = 0;
ImageImportFilterType::SizeType size;
size[0] = width;
size[1] = height;
size[2] = slices;
ImageImportFilterType::RegionType region;
region.SetSize(size);
region.SetIndex(start);
FixedImportFilter->SetRegion(region);
double origin[Dimension];
origin[0] = 0.0;
origin[1] = 0.0;
origin[2] = 0.0;
FixedImportFilter->SetOrigin(origin);
double spacing[Dimension];
spacing[0] = 1;
spacing[1] = 1;
spacing[2] = 1;
FixedImportFilter->SetSpacing(spacing);
FixedImportFilter->SetImportPointer(FixedImageArray, nump,true);

Then the output of FixedImportFilter is image of ITK type built from the source binary file.

DVF file

For DVF file, each voxel has three float numbers representing movements in three directions. So I have to combine them into a vector form. In this part I use ComposeImageFilter class

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
 
typedef itk::Vector< PixelType, Dimension >             VectorPixelType;
typedef itk::Image< PixelType, Dimension >       ImageType;
typedef itk::Image< VectorPixelType, Dimension >    DisplacementFieldType;
PixelType *DVFArrayX = new PixelType[nump];
PixelType *DVFArrayY = new PixelType[nump];
PixelType *DVFArrayZ = new PixelType[nump];
std::ifstream InitDVFFile(argv[3], std::ios::binary | std::ios::in);
InitDVFFile.read((char*)DVFArrayX, nump * sizeof(float));
InitDVFFile.read((char*)DVFArrayY, nump * sizeof(float));
InitDVFFile.read((char*)DVFArrayZ, nump * sizeof(float));
InitDVFFile.close();
ImageImportFilterType::Pointer DVFImportFilterX = ImageImportFilterType::New();
ImageImportFilterType::Pointer DVFImportFilterY = ImageImportFilterType::New();
ImageImportFilterType::Pointer DVFImportFilterZ = ImageImportFilterType::New();
DVFImportFilterX->SetRegion(region);
DVFImportFilterX->SetOrigin(origin);
DVFImportFilterX->SetSpacing(spacing);
DVFImportFilterY->SetRegion(region);
DVFImportFilterY->SetOrigin(origin);
DVFImportFilterY->SetSpacing(spacing);
DVFImportFilterZ->SetRegion(region);
DVFImportFilterZ->SetOrigin(origin);
DVFImportFilterZ->SetSpacing(spacing);
DVFImportFilterX->SetImportPointer(DVFArrayX, nump,true);
DVFImportFilterY->SetImportPointer(DVFArrayY, nump, true);
DVFImportFilterZ->SetImportPointer(DVFArrayZ, nump, true);
typedef itk::ComposeImageFilter< ImageType, DisplacementFieldType> ImageToVectorImageFilterType;
ImageToVectorImageFilterType::Pointer DVFImportFilter = ImageToVectorImageFilterType::New();
DVFImportFilter->SetInput(0, DVFImportFilterX->GetOutput());
DVFImportFilter->SetInput(1, DVFImportFilterY->GetOutput());
DVFImportFilter->SetInput(2, DVFImportFilterZ->GetOutput());
DVFImportFilter->Update();

Then the output of DVFImportFilter is image of ITK type which is image of verctors.

different types of vector images

– VectorImage : Each pixel represents k measurements, each of datatype TPixel. The memory organization of the resulting image is as follows: … Pi0 Pi1 Pi2 Pi3 P(i+1)0 P(i+1)1 P(i+1)2 P(i+1)3 P(i+2)0 … where Pi0 represents the 0th measurement of the pixel at index i.
Conceptually, a VectorImage< TPixel, 3 > is the same as a Image< VariableLengthVector< TPixel >, 3 >. The difference lies in the memory organization. The latter results in a fragmented organization with each location in the Image holding a pointer to an VariableLengthVector holding the actual pixel. The former stores the k pixels instead of a pointer reference, which apart from avoiding fragmentation of memory also avoids storing a 8 bytes of pointer reference for each pixel. The parameter k can be set using SetVectorLength.

GPU registration

In order to accelerate the calculation, it is best to use GPU.

Convert CPUImage to GPUImage

I use CastImageFilter to convert CPUImage to GPUImage

1
2
3
4
5
6
7
8
9
10
11
12
13
 
typedef float                                         InternalPixelType; 
typedef itk::GPUImage< InternalPixelType, Dimension > InternalImageType;
typedef itk::CastImageFilter< ImageType, InternalImageType >  ImageCasterType;
ImageCasterType::Pointer fixedImageCaster = ImageCasterType::New();
fixedImageCaster->SetInput(FixedImportFilter->GetOutput());
InternalImageType::Pointer GPUFixedImage = fixedImageCaster->GetOutput();

typedef itk::Vector< InternalPixelType, Dimension >    InternalVectorPixelType; 
typedef itk::GPUImage< InternalVectorPixelType, Dimension > InternalVectorImageType;
typedef itk::CastImageFilter< DisplacementFieldType, InternalVectorImageType >  VectorImageCasterType;
VectorImageCasterType::Pointer DVFCaster = VectorImageCasterType::New();
DVFCaster->SetInput(DVFImportFilter->GetOutput());
InternalVectorImageType::Pointer DVF = DVFCaster->GetOutput();

Preform GPU demons registration

ITK provides GPUDemonsRegistrationFilter class, which is almost the same as DemonsRegistrationFilter, except the imput image type.

1
2
3
4
5
6
7
8
9
10
11
12
13
 
typedef itk::GPUImage<  VectorPixelType, Dimension >   DeformationFieldType;
typedef itk::GPUDemonsRegistrationFilter<
    InternalImageType,
    InternalImageType,
    DeformationFieldType> RegistrationFilterType;

RegistrationFilterType::Pointer GPUDemonsFilter = RegistrationFilterType::New();

GPUDemonsFilter->SetFixedImage(GPUFixedImage);
GPUDemonsFilter->SetMovingImage(GPUMovingImage);
GPUDemonsFilter->SetInitialDisplacementField(DVF);
GPUDemonsFilter->SetNumberOfIterations(75);
GPUDemonsFilter->SetStandardDeviations(1);

Output

Since the output file storage format is different from vectorimage, I have to transfer first.

1
2
3
4
5
6
7
8
9
10
11
12
 
PixelType *DisplacementArray = new PixelType[Dimension * nump];
memcpy(DisplacementArray, GPUDemonsFilter->GetOutput()->GetBufferPointer(), Dimension* nump * sizeof(float));
PixelType *OutputDisplacementArray = new PixelType[Dimension * nump];
for (unsigned int i = 0; i < nump; i++) {
    OutputDisplacementArray[i] = DisplacementArray[i * 3];
    OutputDisplacementArray[i + nump] = DisplacementArray[i * 3 + 1];
    OutputDisplacementArray[i + 2 * nump] = DisplacementArray[i * 3 + 2];
}

std::ofstream DisplacementImageFile(argv[5], std::ios::binary | std::ios::out);
DisplacementImageFile.write((char*)OutputDisplacementArray, Dimension * nump * sizeof(float));
DisplacementImageFile.close();

微海报分享

近期文章

近期评论

标签

热门

文章归档

分类目录

功能