[ITK] [ITK-users] Writing from an external buffer to VectorImage
Kabelitz, Gordian
Gordian.Kabelitz at medma.uni-heidelberg.de
Mon Jul 3 07:22:49 EDT 2017
Hello,
I want to provide a solution to the question I had. In the case that someone else need an code example.
Just to give an idea how the implementation can look like.
The main reads an image and passes this image to the function that calls the kernel.
After the kernel is done the result is copied to the new vector image [1].
int main(int argc, char* argv[])
{
createImage();
auto nrrdIO = NRRDIOType::New();
auto reader = ReaderType::New();
reader->SetFileName("newImage.nrrd");
reader->SetImageIO(nrrdIO);
try
{
reader->Update();
}
catch (itk::ExceptionObject e)
{
std::cerr << "Exception caught!\n";
std::cerr << e << std::endl;
}
auto image = reader->GetOutput();
auto size = image->GetLargestPossibleRegion().GetSize();
int sz[3] = {size[0],size[1],size[2]};
auto totalVoxel = sz[0] * sz[1] * sz[2];
auto cuda = CudaTest();
auto out = new float4[totalVoxel];
cuda.runKernel(out, image->GetBufferPointer(), sz);
auto vecImage = VectorImageType::New();
VectorImageType::IndexType index;
index[0] = 0;
index[1] = 0;
index[2] = 0;
VectorImageType::RegionType region(index,size);
vecImage->SetRegions(region);
vecImage->SetVectorLength(3);
vecImage->Allocate();
// copy image buffer to vecImage, component by component
auto vecBuffer = vecImage->GetBufferPointer();
auto j = 0;
for (auto i = 0; i < totalVoxel; ++i)
{
vecBuffer[j] = out[i].x; j++;
vecBuffer[j] = out[i].y; j++;
vecBuffer[j] = out[i].z; j++;
}
auto writer = WriterType::New();
writer->SetInput(vecImage);
writer->SetFileName("out.nrrd");
writer->SetImageIO(nrrdIO);
try
{
writer->Update();
}
catch (itk::ExceptionObject e)
{
std::cerr << "Exception caught!\n";
std::cerr << e << std::endl;
}
delete[] out;
}
The function that starts the kernel looks like this:
void CudaTest::runKernel(float4* out, float* data, int* size)
{
auto totalVoxel = size[0] * size[1] * size[2];
// copy image to texture memory
auto channelDesc = cudaCreateChannelDesc<float>();
cudaArray *cuArray;
auto extent = make_cudaExtent(size[0], size[1], size[2]);
gpuErrChk(cudaMalloc3DArray(&cuArray, &channelDesc, extent));
cudaMemcpy3DParms copyparms = { 0 };
copyparms.extent = extent;
copyparms.dstArray = cuArray;
copyparms.kind = cudaMemcpyHostToDevice;
copyparms.srcPtr = make_cudaPitchedPtr((void*)data, size[0] * sizeof(float), size[0], size[1]);
gpuErrChk(cudaMemcpy3D(©parms));
// set texture configuration and bind array to texture
tex_volume.addressMode[0] = cudaAddressModeClamp;
tex_volume.addressMode[1] = cudaAddressModeClamp;
tex_volume.addressMode[2] = cudaAddressModeClamp;
tex_volume.filterMode = cudaFilterModeLinear;
tex_volume.normalized = false;
gpuErrChk(cudaBindTextureToArray(tex_volume, cuArray, channelDesc));
// allocate memory for gradient image data on the device
float4 *dev_gradientImage;
gpuErrChk(cudaMalloc(&dev_gradientImage, totalVoxel * sizeof(float4)));
// copy image dimension to device
int* dev_dimension;
gpuErrChk(cudaMalloc(&dev_dimension, 3 * sizeof(int)));
gpuErrChk(cudaMemcpy(dev_dimension, size, 3 * sizeof(int), cudaMemcpyHostToDevice));
// start computation for the gradient image
dim3 blockDim(16, 8, 8); //<- threads in block, change according to used GPU, max 1024 threads in a single block
dim3 grid((size[0] / blockDim.x) + 1, (size[1] / blockDim.y) + 1, (size[2] / blockDim.z) + 1);
computeGradientImage << <grid, blockDim >> >(dev_gradientImage, dev_dimension);
gpuErrChk(cudaGetLastError());
gpuErrChk(cudaDeviceSynchronize());
// unbind texture
gpuErrChk(cudaUnbindTexture(tex_volume));
gpuErrChk(cudaFreeArray(cuArray));
gpuErrChk(cudaMemcpy(out, dev_gradientImage, totalVoxel * sizeof(float4), cudaMemcpyDeviceToHost));
printf("Kernel is done.\n");
}
The kernel used in this example.
__global__ void computeGradientImage(float4* gradientImage, int* dimension)
{
// every thread computes the float4 voxel with theta,phi,magnitude from gradient image
int idx = blockIdx.x * blockDim.x + threadIdx.x;
int idy = blockIdx.y * blockDim.y + threadIdx.y;
int idz = blockIdx.z * blockDim.z + threadIdx.z;
if (idx < dimension[0] && idy < dimension[1] && idz < dimension[2])
{
// define sobel filter for each direction
float x_000 = tex3D(tex_volume, idx - 1, idy - 1, idz - 1);
float x_001 = tex3D(tex_volume, idx - 1, idy - 1, idz );
float x_002 = tex3D(tex_volume, idx - 1, idy - 1, idz + 1);
float x_010 = tex3D(tex_volume, idx - 1, idy, idz - 1);
float x_011 = tex3D(tex_volume, idx - 1, idy, idz );
float x_012 = tex3D(tex_volume, idx - 1, idy, idz + 1);
float x_020 = tex3D(tex_volume, idx - 1, idy + 1, idz - 1);
float x_021 = tex3D(tex_volume, idx - 1, idy + 1, idz );
float x_022 = tex3D(tex_volume, idx - 1, idy + 1, idz + 1);
float x_100 = tex3D(tex_volume, idx, idy - 1, idz - 1);
float x_101 = tex3D(tex_volume, idx, idy - 1, idz );
float x_102 = tex3D(tex_volume, idx, idy - 1, idz + 1);
float x_110 = tex3D(tex_volume, idx, idy, idz - 1);
// float x_111 = tex3D(tex_volume, idx, idy, idz );
float x_112 = tex3D(tex_volume, idx, idy, idz + 1);
float x_120 = tex3D(tex_volume, idx, idy + 1, idz - 1);
float x_121 = tex3D(tex_volume, idx, idy + 1, idz );
float x_122 = tex3D(tex_volume, idx, idy + 1, idz + 1);
float x_200 = tex3D(tex_volume, idx + 1, idy - 1, idz - 1);
float x_201 = tex3D(tex_volume, idx + 1, idy - 1, idz );
float x_202 = tex3D(tex_volume, idx + 1, idy - 1, idz + 1);
float x_210 = tex3D(tex_volume, idx + 1, idy, idz - 1);
float x_211 = tex3D(tex_volume, idx + 1, idy, idz );
float x_212 = tex3D(tex_volume, idx + 1, idy, idz + 1);
float x_220 = tex3D(tex_volume, idx + 1, idy + 1, idz - 1);
float x_221 = tex3D(tex_volume, idx + 1, idy + 1, idz );
float x_222 = tex3D(tex_volume, idx + 1, idy + 1, idz + 1);
// derivative in x direction
auto centerX = - 1.f*x_000 - 3.f*x_010 - 1.f*x_020
- 3.f*x_001 - 6.f*x_011 - 3.f*x_021
- 1.f*x_002 - 3.f*x_012 - 1.f*x_022
+ 1.f*x_200 + 3.f*x_210 + 1.f*x_220
+ 3.f*x_201 + 6.f*x_211 + 3.f*x_221
+ 1.f*x_202 + 3.f*x_212 + 1.f*x_222;
// derivative in y direction
auto centerY = - 1.f*x_000 - 3.f*x_100 - 1.f*x_200
- 3.f*x_001 - 6.f*x_101 - 3.f*x_201
- 1.f*x_002 - 3.f*x_102 - 1.f*x_202
+ 1.f*x_020 + 3.f*x_120 + 1.f*x_220
+ 3.f*x_021 + 6.f*x_121 + 3.f*x_221
+ 1.f*x_022 + 3.f*x_122 + 1.f*x_222;
// derivative in z direction
auto centerZ = - 1.f*x_000 - 3.f*x_100 - 1.f*x_200
- 3.f*x_010 - 6.f*x_110 - 3.f*x_210
- 1.f*x_020 - 3.f*x_120 - 1.f*x_220
+ 1.f*x_002 + 3.f*x_102 + 1.f*x_202
+ 3.f*x_012 + 6.f*x_112 + 3.f*x_212
+ 1.f*x_022 + 3.f*x_122 + 1.f*x_222;
// magnitude of each voxels gradient
auto magn = sqrtf(centerX*centerX + centerY*centerY + centerZ*centerZ);
gradientImage[idx + dimension[0] * (idy + idz * dimension[1])] = make_float4(centerX, centerY, centerZ, magn);
}
}
with kind regards,
Gordian
[1]: https://itk.org/Doxygen/html/classitk_1_1VectorImage.html
Von: Dženan Zukić [mailto:dzenanz at gmail.com]
Gesendet: Mittwoch, 28. Juni 2017 17:17
An: Kabelitz, Gordian
Cc: insight-users at itk.org
Betreff: Re: [ITK-users] Writing from an external buffer to VectorImage
Hi Gordian,
this approach looks like it should work. What is wrong with it?
Regards,
Dženan Zukić, PhD, Senior R&D Engineer, Kitware (Carrboro, N.C.)
On Wed, Jun 28, 2017 at 9:51 AM, Kabelitz, Gordian <Gordian.Kabelitz at medma.uni-heidelberg.de<mailto:Gordian.Kabelitz at medma.uni-heidelberg.de>> wrote:
Hello,
i computed a gradient with my own function and as a result a pointer to an image buffer is provided. I know the size, origin and spacing of the gradient component image.
I want to copy the gradient image into an itk::VectorImage with the components for the x,y,z gradients.
The way I copied the image to the GPU is that I retrieved the buffer pointer from my input image and use the pointer to copy the image data to the GPU.
I used the way proposed in [1]. The computeGradientImage method is listed at the end of this mail.
[...]
// get float pointer to image data
ImageType::Pointer image = reader->GetOutput();
Image->Update();
float* data = image.GetBufferPointer();
// copy image data to GPU texture memory (this works)
gpu_dev->setVoxels(dimension, voxelSize, data);
[...]
computeGradientImage<<<n,m>>> (dev_gradientImage, dimension);
// copy resulting gradientImage to host variable
float4* host_gradientImage;
cudaMemcpy(host_gradient, dev_gradientImage, numberOfVoxels*sizeof(float4));
--> Pseudo Code <--
// Now I want to reverse the copy process. I have a float4 image and want to copy this into a itk::VectorImage with VariableVectorLength of 3 (skipping the magnitude value).
[...] -> size, spacing, origin, region definition
Itk::VectorImageType vecImage = VectorImageType::New();
vecImage->setRegion(region);
vecImage ->SetVectorLength(3);
vecImage->Allocate();
// copy image buffer to vecImage, component by component
auto vecBuffer = vecImage->getBufferPointer();
auto j = 0;
for (i=0; i<totalVoxel; ++i)
{
vecbuffer[j] = host_gradient[i].x; j++;
vecbuffer[j] = host_gradient[i].y; j++;
vecbuffer[j] = host_gradient[i].z; j++;
}
// save vecImage as nrrd image
[...]
I haven't found a way to achieve my idea.
Are there any suggestions or examples?
As far I can see I cannot use the itk::ImportImageFilter.
Thank you for any suggestions.
With kind regards,
Gordian
[1]: https://itk.org/CourseWare/Training/GettingStarted-V.pdf
void computeGradientImage(float4* gradientImage, int* dimension)
{
// every thread computes the float4 voxel with theta,phi,magnitude from gradient image
int idx = blockIdx.x * blockDim.x + threadIdx.x;
int idy = blockIdx.y * blockDim.y + threadIdx.y;
int idz = blockIdx.z * blockDim.z + threadIdx.z;
if (idx < dimension[0] && idy < dimension[1] && idz < dimension[2])
{
// define sobel filter for each direction
[...]
// run sobel on image in texture memory for each direction and put result into a float4 image
gradientImage[idx + dimension[0] * (idy + idz * dimension[1])] = make_float4(sobelX, sobelY, sobelZ, magn);
}
}
_____________________________________
Powered by www.kitware.com<http://www.kitware.com>
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php
Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ
Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/insight-users
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20170703/12148ad3/attachment-0001.html>
-------------- next part --------------
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php
Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ
Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/insight-users
More information about the Community
mailing list