How to extract a volume around a voxel in a 3D image?

Suppose I have a 3D image of 250*250*250 size.
I want to extract a volume of 25*25*25 around a particular voxel that could be defined by me.
Is there any example that I could look into ??

In ITK you’d use the ExtractImageFilter. Here’s the documentation for it:

https://itk.org/Doxygen/html/classitk_1_1ExtractImageFilter.html

and here’s an example program that uses it:

https://itk.org/Doxygen/html/SphinxExamples_2src_2Core_2Common_2CropImageBySpecifyingRegion_2Code_8cxx-example.html#_a6

In SimpleITK python you can use do something spiffy like:

img2 = img1[cx-25:cx+25, cy-25:cy+25, cz-25:cz+25]

(except you’d probably want to do some bounds checking)

Dave, how would you implement this in C++/SimpleITK ? I am asking this because seem more easily than ITK implementation …

Hi,

I’ve been using the CropImageFilter to do something similar.

https://itk.org/Doxygen/html/classitk_1_1CropImageFilter.html

You need to specify the opposites corners of you area.

cropmageFilter->SetInput(imgVessels);
cropImageFilter->SetLowerBoundaryCropSize( lowerSize);
cropImageFilter->SetUpperBoundaryCropSize( upperSize);
cropImageFilter->Update();

Jonas, thank you for your idea. In my case, I have a 3d point (taken with vtkCellPicker), and with this 3d point I want to take only a limited volume around this point …

A little detail: how can you achieved imgVessels ? Sorry for this simple questions, I am a beginner in ITK/SimpleITK.

Sorry about imgVessels. That was my input image at the time, but replace it with your input Image.

lowerSize and upperSize are itk::Size<3,int> Objects.

lowerSize[0] = PointIndexX - offsetX;
lowerSize[1] = PointIndexY - offsetY;
lowerSize[2] = PointIndexZ - offsetZ;

upperSize[0] = PointIndexX + offsetX;
upperSize[1] = PointIndexY + offsetY;
upperSize[2] = PointIndexZ + offsetZ;

this should work as long as you check that you are not out of bond.

Unfortunately, this is a SimpleITK/Python only feature. To do extract a sub-image in C++, you have to use the CropImageFilter (as @Jonas demonstrates).

I will try your idea. If I don’t push my luck, can you tell me how to check the bound limits ? (I am trying to learn ITK/SimpleITK from this forum too).

The best way to learn is to read the documentation and look at examples.

here is the documentation on regions.
https://itk.org/Doxygen/html/classitk_1_1ImageRegion.html

You also have a lot of code example such as:
https://itk.org/ITKExamples/search.html?q=region

The documentation is well made in that regard. You should abuse it as much as possible.

Based on your guiding lines, I wrote:

	typedef float InternalPixelType;
	const unsigned int nDimension = 3;
	typedef itk::Image<InternalPixelType, nDimension> InputImageType;
	typedef itk::Image<InternalPixelType, nDimension> OutputImageType;
	typedef itk::CropImageFilter<InputImageType, OutputImageType> CropImageFilterType;
	typedef itk::Image<InternalPixelType, nDimension> InternalImageType;
	InternalImageType::Pointer itkImage = dynamic_cast <InternalImageType*>(img.GetITKBase());
	CropImageFilterType::Pointer cropFilter	= CropImageFilterType::New();
	cropFilter->SetInput(itkImage);
	// do nothing with cropFilter for now
	cropFilter->Update();

where img is sitk::Image type. But I got an error at cropFilter->Update(); line:

itkExceptionMacro(<< "Input " << requiredInputName << " is required but not set.");

on ProcessObject ::VerifyPreconditions method … obviously I’ve done something wrong … but what ?

I have to mention that I seek on https://itk.org/ITKExamples/search.html?q=region and in ITK examples to figure out what happen here, but I didn’ catch the issue … can you lead guide me a little bit ?

Hello @flaviu2,
It looks like the dynamic cast failed (the example you a likely based your code on shows that you need to check for this).

Are you only converting to an ITK image in order to crop the image? If yes, then you could use the SimpleITK crop filter.

Yes you had right … the itkImage object has been NULL, I don’t understand why …

Now I understand the difference between SimpleITK and ITK. With SimpleITK I have:

	sitk::CropImageFilter crop;
	crop.SetLowerBoundaryCropSize(lower);
	crop.SetUpperBoundaryCropSize(upper);
	crop.Execute(img);

much more simpler :slight_smile:

1 Like

Please tell me why I got out of bounds on this filter:

sitk::Image img;

int cx = img.GetWidth();  // 512
int cy = img.GetHeight();  // 512
int cz = img.GetDepth();  // 44

int x = 326;
int y = 351;
int z = 26;

std::vector<unsigned int> lower;
lower.push_back(x + 2);
lower.push_back(y + 2);
lower.push_back(z + 2);
std::vector<unsigned int> upper;
upper.push_back(x - 2);
upper.push_back(y - 2);
upper.push_back(z - 2);

sitk::CropImageFilter crop;
crop.SetLowerBoundaryCropSize(lower);
crop.SetUpperBoundaryCropSize(upper);
crop.Execute(img);

But I got this error:

itkExceptionMacro("The input image's size " << input_sz << " is less than the total of the crop size!");

That point, x,y,z seem to be within img bounds … or not ? Can you help me ?

The documentation for SetLowerBondaryCropSize says:

Set/Get the cropping sizes for the upper and lower boundaries.

So it is the size to reduce the volume by. For you case it would be:

std::vector<unsigned int> crop_size={2,2,2}

Yes, it goes now. The sitk::Image is cropped:

	sitk::Image imgTemp(img);
	sitk::CropImageFilter crop;
	crop.SetLowerBoundaryCropSize(lower);
	crop.SetUpperBoundaryCropSize(upper);
	img = crop.Execute(imgTemp);

I wonder if there is some simpleitk method to crop only a small part of volume, let say a bone or a tooth …

I have tried to save the copped image into a jpg/bmp/etc image but the sitk::writer told me that “JPEG Writer can only write 2-dimensional images