Filter a 3D image slice by slice and store results

Beginner Questions
1

Dear all,
I have a 3D volume, I want to iterate over all 2D slices, apply a convolutional filter on each slice and then put them again together as a 3D volume and store this 3D volume.

Somehow, I am failing to do so and I do not know what might be the problem. I wanted to iterate over the image slice by slice and consequently update the image. But it fails at the second slice/loop Thats what I have so far:

//create the convolutional kernel (that works, should be OK, it gives me a 3D image with size (5x5x1)
ImageType::Pointer kernel = ImageType::New();
CreateKernel(kernel)
using PasteFilterType = itk::PasteImageFilter<ImageType>;
//extract an image slice
using ExtractFilterType = itk::ExtractImageFilter<ImageType, ImageType>;
ExtractFilterType::Pointer extractFilter = ExtractFilterType::New();
extractFilter->SetDirectionCollapseToSubmatrix();
//image to change
ImageType *     inputImage = image;
const typename ImageType::RegionType& imageRegion = image->GetLargestPossibleRegion();
const typename ImageType::SizeType& imageRegionSize = imageRegion.GetSize();
ImageType::RegionType inputRegion = inputImage->GetBufferedRegion();
ImageType::SizeType   size = inputRegion.GetSize();
size[2] = 1; // we extract along z direction
ImageType::IndexType start = inputRegion.GetIndex();
//start at slice 0
unsigned int   sliceNumber = 0;
//use convolutional filter
using FilterType = itk::ConvolutionImageFilter<ImageType>;
ImageType::RegionType desiredRegion;
FilterType::Pointer convolutionFilter = FilterType::New();
//iterate over all slices
for (int depth = 0; depth < imageRegionSize[2]; depth++) {
//set slice nr
			sliceNumber = depth;
			start[2] = sliceNumber;
			//get 2D slice
			desiredRegion.SetSize(size);
			desiredRegion.SetIndex(start);

			extractFilter->SetExtractionRegion(desiredRegion);
			//paste filtered slice
			PasteFilterType::Pointer pasteFilter = PasteFilterType::New();
			extractFilter->SetInput(inputImage);
			//convolve image
			convolutionFilter->SetInput(extractFilter->GetOutput());
			convolutionFilter->SetKernelImage(kernel);
//insert filtered slice
			pasteFilter->SetSourceImage(convolutionFilter->GetOutput());
			pasteFilter->SetDestinationImage(inputImage);
			
			pasteFilter->SetDestinationIndex(start);
			
			pasteFilter->SetSourceRegion(convolutionFilter->GetOutput()->GetBufferedRegion());
			pasteFilter->Update();

			
			inputImage= pasteFilter->GetOutput();
		}
		writer->SetFileName("E:/IBSI2/ibsi_2_digital_phantom/gabor1.nii");
		writer->SetInput(inputImage);
		writer->Update();
		
	}

Thanks a lot for the help!

2

Hello,

With SimpleITK 2.0 we just added some infrastructure and examples for this exact problem.

The current ITK PasteImageFilter does not allow a 2D Image to be pasted into a 3D Image. But I extended it into an NPasteImageFilter which is available in ITK from the SimpleITKFitlers remote module. I was just starting to patch the ITK PasteImageFilter with the added functionality.

I am having problem following your posted code. You can edit the most to reformat it, or post a link to your code.

4

Thanks already for your answer. Sorry, I see, the code looks very messy. I tried to make t more readable:

//create the convolutional kernel (that works, should be OK, it gives me a 3D image with size (5x5x1)):

ImageType::Pointer kernel = ImageType::New();
CreateKernel(kernel)

//define the image I want to change (imageType is a 3D image)

ImageType * inputImage = image;
const typename ImageType::RegionType& imageRegion = image->GetLargestPossibleRegion();
const typename ImageType::SizeType& imageRegionSize = imageRegion.GetSize();
ImageType::RegionType inputRegion = inputImage->GetBufferedRegion();
ImageType::SizeType size = inputRegion.GetSize();

//define extract image filter to extract an image slice

using ExtractFilterType = itk::ExtractImageFilter<ImageType, ImageType>;
ExtractFilterType::Pointer extractFilter = ExtractFilterType::New();
extractFilter->SetDirectionCollapseToSubmatrix();
size[2] = 1; // we extract along z direction
ImageType::IndexType start = inputRegion.GetIndex();

// now start at slice 0

unsigned int sliceNumber = 0;

//define convolutional and paste filter
using PasteFilterType = itk::PasteImageFilter;

using FilterType = itk::ConvolutionImageFilter;
ImageType::RegionType desiredRegion;
FilterType::Pointer convolutionFilter = FilterType::New();

//iterate over all slices

for (int depth = 0; depth < imageRegionSize[2]; depth++) {
//set slice nr
sliceNumber = depth;
start[2] = sliceNumber;
//get 2D slice
desiredRegion.SetSize(size);
desiredRegion.SetIndex(start);

extractFilter->SetExtractionRegion(desiredRegion);
//paste filtered slice
PasteFilterType::Pointer pasteFilter = PasteFilterType::New();
extractFilter->SetInput(inputImage);
//convolve image
convolutionFilter->SetInput(extractFilter->GetOutput());
convolutionFilter->SetKernelImage(kernel);

//insert filtered slice
pasteFilter->SetSourceImage(convolutionFilter->GetOutput());
pasteFilter->SetDestinationImage(inputImage);
pasteFilter->SetDestinationIndex(start);
pasteFilter->SetSourceRegion(convolutionFilter->GetOutput()->GetBufferedRegion());
pasteFilter->Update();

inputImage= pasteFilter->GetOutput();
}

5

To help debug, I would construct and configure all filters inside the for loop.

Also, your radius of [5,5,1] is equal to a diameter of [11,11,3] so the boundary condition is being used in the Z direction and the efficiency is being reduced.

1 Like
6

@Elli Please use triple backticks when including source code (see more options here).

2 Likes