How can I apply a 2D mask in a 3D image? - SliceBySliceImageFilter

cristinamontserrat · October 17, 2022, 8:59am

Hi!

I am trying to mask a 3D image with a 2D mask with the itkSliceBySliceImageFilter. Here is my code:

// We create the filter mask
using MaskFilterType = itk::MaskImageFilter<ITKImgLabelType2D, ITKImgLabelType2D, ITKImgLabelType2D>;
MaskFilterType::Pointer maskFilter3D = MaskFilterType::New();
maskFilter3D->SetInput(image);
maskFilter3D->SetMaskImage(mask);
maskFilter3D->Update();

// Apply a 2D mask to every slice of a 3D image
typedef itk::SliceBySliceImageFilter < ITKImgInType, ITKImgLabelType, MaskFilterType, MaskFilterType, ITKImgLabelType2D, ITKImgLabelType2D > SliceBySliceImageFilter;
SliceBySliceImageFilter::Pointer sliceBySliceFilter = SliceBySliceImageFilter::New();
sliceBySliceFilter->SetFilter(maskFilter3D);
sliceBySliceFilter->SetInput(&image3D);

try
{
sliceBySliceFilter->Update();
}
catch (const itk::ExceptionObject& err)

Where 'ITKImgLabelType2D = itk::Image<unsigned char, 2>’ , ‘ITKImgInType = itk::Image<short, 3>’ and ‘ITKImgLabelType = itk::Image<unsigned char, 3>’

When I execute sliceBySliceFilter->Update(); I get the following error:

‘itk::ERROR: SliceBySliceImageFilter(000001630B273560): At least 2 of the first 2 indexed inputs are required but only 1 are specified. The required inputs are expected to be the first inputs’

And I have no idea of what could it be. How can I solve it? Could I obtain what I want in a different way?

My ITK version is: 4.8.2

Thank you very much,
Cristina

dzenanz · October 17, 2022, 1:01pm

Maybe try sliceBySliceFilter->SetInput(1, mask);? Looking at the code, SliceBySliceImageFilter should be able to work with multiple inputs.

dzenanz · October 17, 2022, 1:06pm

But it expects all the inputs to have the same size (requested region). Depending on what you are trying to accomplish, it might be easier to construct a 3D mask using TileImageFilter, and use 3D variant of MaskImageFilter.

cristinamontserrat · October 18, 2022, 10:47am

Thank you for the answer.

Firstly

I have tried the first advice:

auto mask = ITKImgLabelType2D::New();
ITKImgLabelType2D::RegionType region = a_pmask->GetLargestPossibleRegion();

mask->SetRegions(region);
mask->Allocate();

ITKImgLabelType2D::SizeType regionSize = region.GetSize();
itk::ImageRegionIterator imageIterator(mask, region);

while (!imageIterator.IsAtEnd())
{
imageIterator.Set(itk::NumericTraits::max());
++imageIterator;
}

// We create the filter mask
using MaskFilterType = itk::MaskImageFilter<ITKImgLabelType2D, ITKImgLabelType2D, ITKImgLabelType2D>;
MaskFilterType::Pointer maskFilter3D = MaskFilterType::New();
maskFilter3D->SetInput(a_pmask);
maskFilter3D->SetMaskImage(mask);

// Apply a 2D mask to every slice of a 3D image
typedef itk::SliceBySliceImageFilter< ITKImgInType, ITKImgLabelType, MaskFilterType, MaskFilterType, ITKImgLabelType2D, ITKImgLabelType2D > SliceBySliceImageFilter;
SliceBySliceImageFilter::Pointer sliceBySliceFilter = SliceBySliceImageFilter::New();
sliceBySliceFilter->SetFilter(maskFilter3D);
sliceBySliceFilter->SetInput(1, &a_pimage3D);

And I get another error, which is:

itk::ERROR: SliceBySliceImageFilter(0000019653A533C0): Input Primary is required but not set.

Secondly

I have tried to construct a 3D mask using TileImageFilter , and use 3D variant of MaskImageFilter.

Here is my code:

constexpr unsigned int InputDimension = 2;
constexpr unsigned int OutputDimension = 3;

using PixelType = unsigned char;
using InputImageType = itk::Image<PixelType, InputDimension>;
using OutputImageType = itk::Image<PixelType, OutputDimension>;

using FilterType = itk::TileImageFilter<InputImageType, OutputImageType>;
auto filter = FilterType::New();

itk::FixedArray<unsigned int, OutputDimension> layout;
layout[0] = 2;
layout[1] = 2;
layout[2] = 0;

filter->SetLayout(layout);

ITKImgInType::SizeType insize = a_pimage3D.GetRequestedRegion().GetSize();
int nImages = insize[2];

for (int ii = 0; ii <= nImages - 1; ++ii)
{
InputImageType::Pointer input = a_pmask;
filter->SetInput(ii, input);
}

constexpr PixelType defaultValue = 128;
filter->SetDefaultPixelValue(defaultValue);

filter->GetOutput()->SetSpacing(a_pimage3D.GetSpacing());
filter->GetOutput()->SetOrigin(a_pimage3D.GetOrigin());
filter->GetOutput()->SetDirection(a_pimage3D.GetDirection());

if (f_verbose)
{
a_file << "Spacing 3D mask after: " << filter->GetOutput()->GetSpacing() << std::endl;
a_file << "Spacing 3D image: " << a_pimage3D.GetSpacing() << std::endl << std::endl;
a_file << "Origin 3D mask after: " << filter->GetOutput()->GetOrigin() << std::endl;
a_file << "Origin 3D image: " << a_pimage3D.GetOrigin() << std::endl << std::endl;
}

// Apply a 3D mask to a 3D image
using MaskFilterType = itk::MaskImageFilter<ITKImgInType, ITKImgLabelType, ITKImgLabelType>;
MaskFilterType::Pointer maskFilter3D = MaskFilterType::New();
maskFilter3D->SetInput(&a_pimage3D);
maskFilter3D->SetMaskImage(filter->GetOutput());

try
{
maskFilter3D->Update();
}
catch (const itk::ExceptionObject& err)

I don’t understand what’s wrong, because if I print the origin and the spacing of the a_image3D and the filter, as you can see in my code, I get the following information:

Spacing 3D mask after: [0.847656, 0.847656, 2.98663]
Spacing 3D image: [0.847656, 0.847656, 2.98663]

Origin 3D mask after: [-208.208, -251.595, 1148.21]
Origin 3D image: [-208.208, -251.595, 1148.21]

So, the spacing and the origin is the same. But then, when I execute maskFilter3D->Update(); I get the following error:

ExceptionObject caught !Alma::SegmentationLib::ApplyMask3D
C:\Dev\ITK\Modules\Core\Common\include\itkImageToImageFilter.hxx:248:
itk::ERROR: MaskImageFilter(00000231EBE04390): Inputs do not occupy the same physical space!
InputImage Origin: [-2.0820781e+02, -2.5159518e+02, 1.1482096e+03], InputImage_1 Origin: [-2.0820781e+02, -2.5159518e+02, 0.0000000e+00]
Tolerance: 8.4765559e-07
InputImage Spacing: [8.4765559e-01, 8.4765559e-01, 2.9866300e+00], InputImage_1 Spacing: [8.4765559e-01, 8.4765559e-01, 1.0000000e+00]
Tolerance: 8.4765559e-07

What should I do?

Thank you very much,
Cristina

dzenanz · October 18, 2022, 2:57pm

From the docs: “The output image have a larger dimension than the input images. This filter can be used to create a volume from a series of inputs by specifying a layout of 1,1,0.” You are using a tile layout 2,2,0. Your mask is probably bigger than the image.

a_file << "Spacing 3D mask after: " << filter->GetOutput()->GetSpacing() << std::endl;

You should print to the same or higher precision as the comparison tolerance (7-8 decimal places) in order to spot a possible problem here.

cristinamontserrat · October 19, 2022, 7:02pm

Thanks!

I’ve changed the tile layout and I also get exactly the same error.

My code now is:

using PixelType = unsigned char;
using InputImageType = ITKImgLabelType2D;
using OutputImageType = ITKImgLabelType;

using FilterType = itk::TileImageFilter<InputImageType, OutputImageType>;
auto filter = FilterType::New();

itk::FixedArray<unsigned int, 3> layout;
layout[0] = 1;
layout[1] = 1;
layout[2] = 0;

filter->SetLayout(layout);

ITKImgInType::SizeType insize = a_pimage3D.GetRequestedRegion().GetSize();
int nImages = insize[2];

unsigned int inputImageNumber = 0;

for (int ii = 1; ii < nImages-1; ++ii)
{
InputImageType::Pointer input = a_pmask;
input->DisconnectPipeline();
filter->SetInput(inputImageNumber++, input);
}

constexpr PixelType defaultValue = 128;
filter->SetDefaultPixelValue(defaultValue);

filter->GetOutput()->SetSpacing(a_pimage3D.GetSpacing());
filter->GetOutput()->SetOrigin(a_pimage3D.GetOrigin());
filter->GetOutput()->SetDirection(a_pimage3D.GetDirection());

if (f_verbose)
{
a_file << std::fixed << std::setprecision(9) << std::endl;
a_file << "Spacing 3D mask after: " << filter->GetOutput()->GetSpacing() << std::endl;
a_file << "Spacing 3D image: " << a_pimage3D.GetSpacing() << std::endl << std::endl;
a_file << "Origin 3D mask after: " << filter->GetOutput()->GetOrigin() << std::endl;
a_file << "Origin 3D image: " << a_pimage3D.GetOrigin() << std::endl << std::endl;
}

// Apply a 3D mask to a 3D image
using MaskFilterType = itk::MaskImageFilter<ITKImgInType, ITKImgLabelType, ITKImgLabelType>;
MaskFilterType::Pointer maskFilter3D = MaskFilterType::New();
maskFilter3D->SetInput1(&a_pimage3D);
maskFilter3D->SetInput2(filter->GetOutput());

try
{
maskFilter3D->Update();
}
catch (const itk::ExceptionObject& err)

What exactly do the layout values 1,1,0 mean?

Could it be possible that the error is related to the different input pixels type of the MaskImageFilter function?

Thank you,
Cristina

dzenanz · October 19, 2022, 7:25pm

Have you tried filter->Update(); before if (f_verbose)? Are size, bufferedRegion and direction same, too? You could inspect that using a debugger, instead of printing to console.

cristinamontserrat · October 20, 2022, 4:29pm

Yes, I have tried that. And I also have checked that the inputs have the same size, bufferedRegion and direction too.

As a possible solution, I have changed the way of resampling the image.
I have replaced this:

filter->GetOutput()->SetSpacing(a_pimage3D.GetSpacing());
filter->GetOutput()->SetOrigin(a_pimage3D.GetOrigin());
filter->GetOutput()->SetDirection(a_pimage3D.GetDirection());

For this:

using ResampleFilterType = itk::ResampleImageFilter<ITKImgLabelType, ITKImgLabelType>;
auto resampleFilter = ResampleFilterType::New();

resampleFilter->SetInput(filter->GetOutput());
resampleFilter->SetOutputSpacing(a_pimage3D.GetSpacing());
resampleFilter->SetOutputOrigin(a_pimage3D.GetOrigin());
resampleFilter->SetOutputDirection(a_pimage3D.GetDirection());

This way, I don’t get the error: Inputs do not occupy the same physical space!, but I do get the following error: Requested region is (at least partially) outside the largest possible region.

Do you know the reason for this?

Thank you for your dedication.
Cristina

dzenanz · October 20, 2022, 4:37pm

You need to invoke filter->Update(), or at least filter->UpdateOutputInformation() before filter->GetOutput() will have meaningful data.

cristinamontserrat · October 25, 2022, 6:00pm

Hi!

I have tried that and it didn’t work for me. But I have also tried another thing, that consists in “resampling” the mask with the ResampleImageFilter instead of manually. This way, I have fixed the error. But now another thing is happening.

I create the 3D mask correctly, but then, when I “resample” it, the 3D mask disappears.

I’m going to explain myself with my code.

First of all I create a 3D mask from the 2D mask:

using PixelType = unsigned char;
using InputImageType = ITKImgLabelType2D;
using OutputImageType = ITKImgLabelType;

using FilterType = itk::TileImageFilter<InputImageType, OutputImageType>;
auto filter = FilterType::New();

itk::FixedArray<unsigned int, 3> layout;
layout[0] = 1;
layout[1] = 1;
layout[2] = 0;

filter->SetLayout(layout);

ITKImgInType::SizeType insize = a_pimage3D.GetRequestedRegion().GetSize();
int nImages = insize[2];

unsigned int inputImageNumber = 0;
InputImageType::Pointer input = a_pmask;

WriteImage(a_pmask, “a_pmask.jpg”, a_file);
for (int ii = 0; ii < nImages; ++ii)
{
input->DisconnectPipeline();
filter->SetInput(inputImageNumber++, a_pmask);
}

try
{
filter->UpdateOutputInformation();
filter->UpdateLargestPossibleRegion();
}
catch (const itk::ExceptionObject& err)
{
if (f_verbose)
{
a_file << “ExceptionObject caught !” << func << std::endl;
a_file << err.what() << std::endl;
}
}

That’s done correctly, because if I print the image I can see the 2D mask in every slice of the 3D mask.
Then, I “resample” the mask:

ITKImgInType::SizeType outputSize = a_pimage3D.GetLargestPossibleRegion().GetSize();
ITKImgLabelType::SpacingType outputSpacing = a_pimage3D.GetSpacing();

using ScalarType = double;

using NearestNeighborInterpolatorType = itk::NearestNeighborInterpolateImageFunction<ITKImgLabelType, ScalarType>;
auto nearestNeighborInterpolator = NearestNeighborInterpolatorType::New();

using ResampleImageFilterType = itk::ResampleImageFilter<ITKImgLabelType, ITKImgLabelType>;
ResampleImageFilterType::Pointer resample = ResampleImageFilterType::New();

resample->SetInput(filter->GetOutput());
resample->SetReferenceImage(&a_pimage3D);
resample->SetInterpolator(nearestNeighborInterpolator);
resample->UseReferenceImageOn(); // spacing, origin and direction of the reference image will be used
resample->SetSize(outputSize);

try
{
resample->UpdateOutputInformation();
resample->UpdateLargestPossibleRegion();
}
catch (const itk::ExceptionObject& err)
{
if (f_verbose)
{
a_file << “ExceptionObject caught !” << func << std::endl;
a_file << err.what() << std::endl;
}
}

Now, if I print the 3D mask, is wholly or partially black. So, when I do apply this mask to my original 3D image it does not do anything.

Am I missing something?

Thank you very much for your time,
Cristina

dzenanz · October 25, 2022, 8:18pm

Why are you resampling? And if you are resampling, you should pay attention to spatial extent of your images, their overlap etc.

You could take a look at a tutorial from earlier this year, or its video recording. Of particular interest would be ITK image grid notebook.

cristinamontserrat · November 7, 2022, 8:56am

Thank you very much!

I actually don’t know why I was resampling. I finally managed to fix the problem by changing the
ResampleImageFilter to ChangeInformationImageFilter.