Posted on | Edited on | In

Introduction

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 FixedImageFile(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();

Posted on | Edited on | In

Work

This week I get the same result as my fellow. As what I expected, the performance of SMEIR after 400 iterations is the same as mSMEIR. Well, my fellow’s paper only show the result of 200 iterations SMEIR. And the show the result of another 200 iterations mSMEIR. This comparison looks good, but it’s meaningless,
In order to get the trajectory of tumor, I generate the lung by XCAT, the tumor by NCAT, and then combine them in MATLAB. And the result seems unusually, there only exists motion in one direction. Let me generate these all by NCAT and see the result.

New command

It needs more work for AAPM. I thought the abstract only needs words. Unfortunately, I need provide supporting document. It’s better to read instruction rather than thinking.

Posted on | Edited on | In

work

This week, I still mainly focus on mSMEIR. It’s hard to judge this code, total in disorder. I don’t like it, but I have to use it. After I finish this project, I would either throw it away or rewrite it.

new command

My supervisor want me to submit AAPM conference, which is close to deadline(March 7). I’m not feeling good now.