[ITK Community] Fwd: Problems with slice extraction from a 3D image for directions different from z

[Luis Ibanez]

Hi Simon,

Thanks for your very clear description of the problem you are facing.

The message that you are receiving is quite common, and it is the
result of the fact that when ITK generates a 2D image from a 3D
image, the direction matrix the 3D image (which is initially a 3x3
matrix) must be reduced to a 2D direction matrix which is a 2x2
matrix).

There are multiple strategies for performing this reduction, none
of which is "optimal" nor "correct", and hence, the choice is left
to the developer.....   so... here you are   :-)

An example of when this happens,
is if you have a 3D direction matrix such as

1 0 0
0 0 -1
0 1 0

and with the simple strategy of keeping the first two dimensions,
we would end up in 2D with a matrix

1 0
0 0

which is singular....
and hence you get the exception that you are observing now.


Printing out the values of your direction matrix
will be enlightening.

std::cout << reader->GetOutput()->GetDirection() << std::endl;



When you opt for:

filter->SetDirectionCollapseToIndentity();

you are forcing the filter to produce the 2D matrix

1 0
0 1

which brings you one step closer to happiness.


What I see out of order in your code is the line

filter->InPlaceOn();

which doesn't work well in the Extract image
filter when you are extracting a 2D image out
of a 3D image.

The InPlaceOn() method is intended for saving
memory in the cases where an output image is
of the same type of an input image.

In your case, the two images are of different
dimension, and hence the output can not be
stored in the input image.



   Hope this helps.


BTW, could you please elaborate more, on what
you mean by: "the slice spacing information is lost"  ?

In principle, the reduction from 3D to 2D fully
collapses a third dimension, so the spacing is
indeed meaningless after the extraction.

One alternative is to Extract, without collapsing
dimension, and in this way the output image is
a 3D (not 2D) image, of a single slice.


    Thanks

         Luis



On Sat, Mar 8, 2014 at 3:01 PM, Simon Wimmer <wimmersimon at gmail.com> wrote:

> Hello ITK community,
>
> I've started to use ITK two weeks ago for an university project and have
> been quite satisfied with what I achieved with it so far.
> But today, I stumbled across a problem that I could not resolve by
> consulting the existing resources on the web and old mailing list threads.
> For our project we have to work with 3D MRI data that is provided as an
> .mhd file.
> We need to extract slices orthogonal to every axis from that data.
> For that matter I consulted the docs and found the example code in
> ImageReadExtractWrite.cxx which makes use of the ExtractImageFilter. This
> method works perfectly for extracting a slice in z direction as the example
> code does.
> So I tried to make the code extract a slice in x direction by simply
> changing the lines
>
> size[2] = 0;
>> InputImageType::IndexType start = inputRegion.GetIndex();
>> const unsigned int sliceNumber = atoi( argv[3] );
>> start[2] = sliceNumber;
>
>
> to
>
>
>> size[0] = 0;
>> InputImageType::IndexType start = inputRegion.GetIndex();
>> const unsigned int sliceNumber = atoi( argv[3] );
>> start[0] = sliceNumber;
>
>
> But this code throws an exception:
>
> Description: itk::ERROR: ExtractImageFilter(037F46F0): Invalid submatrix
>> extracted for collapsed direction.
>
>
> If I change
>
>> filter->SetDirectionCollapseToSubmatrix();
>
> to
>
>>   filter->SetDirectionCollapseToIndentity();
>
> the code runs cleanly, but produces an awkward result image because
> seemingly the slice spacing information is lost.
>
> Do you have any hints on how I could resolve the issue/where I have gone
> wrong or can you point out different approaches I could try to pursue?
>
> The whole code:
>
>> #include "itkImageFileReader.h"
>> #include "itkImageFileWriter.h"
>> #include "itkExtractImageFilter.h"
>> #include "itkImage.h"
>>
>> int main( int argc, char ** argv )
>> {
>>   // Verify the number of parameters in the command line
>>   if( argc < 3 )
>>     {
>>     std::cerr << "Usage: " << std::endl;
>>     std::cerr << argv[0] << " input3DImageFile  output2DImageFile " <<
>> std::endl;
>>     std::cerr << " sliceNumber " << std::endl;
>>     return EXIT_FAILURE;
>>     }
>>   typedef signed short        InputPixelType;
>>   typedef signed short        OutputPixelType;
>>   typedef itk::Image< InputPixelType,  3 >    InputImageType;
>>   typedef itk::Image< OutputPixelType, 2 >    OutputImageType;
>>   typedef itk::ImageFileReader< InputImageType  >  ReaderType;
>>   typedef itk::ImageFileWriter< OutputImageType >  WriterType;
>>   // Here we recover the file names from the command line arguments
>>   //
>>   const char * inputFilename  = argv[1];
>>   const char * outputFilename = argv[2];
>>   ReaderType::Pointer reader = ReaderType::New();
>>   WriterType::Pointer writer = WriterType::New();
>>   reader->SetFileName( inputFilename  );
>>   writer->SetFileName( outputFilename );
>>   typedef itk::ExtractImageFilter< InputImageType,
>>                                    OutputImageType > FilterType;
>>   FilterType::Pointer filter = FilterType::New();
>>   filter->InPlaceOn();
>>   filter->SetDirectionCollapseToSubmatrix();
>>   reader->UpdateOutputInformation();
>>   InputImageType::RegionType inputRegion =
>>            reader->GetOutput()->GetLargestPossibleRegion();
>>   InputImageType::SizeType size = inputRegion.GetSize();
>>   size[0] = 0;
>>   InputImageType::IndexType start = inputRegion.GetIndex();
>>   const unsigned int sliceNumber = atoi( argv[3] );
>>   start[0] = sliceNumber;
>>   InputImageType::RegionType desiredRegion;
>>   desiredRegion.SetSize(  size  );
>>   desiredRegion.SetIndex( start );
>>   filter->SetExtractionRegion( desiredRegion );
>>   filter->SetInput( reader->GetOutput() );
>>   writer->SetInput( filter->GetOutput() );
>>   try
>>     {
>>     writer->Update();
>>     }
>>   catch( itk::ExceptionObject & err )
>>     {
>>     std::cerr << "ExceptionObject caught !" << std::endl;
>>     std::cerr << err << std::endl;
>>     return EXIT_FAILURE;
>>     }
>>
>>   return EXIT_SUCCESS;
>> }
>
>
> Regards,
>
> Simon
>
>
>
> _______________________________________________
> Community mailing list
> Community at itk.org
> http://public.kitware.com/cgi-bin/mailman/listinfo/community
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20140308/24a39ad8/attachment-0002.html>