[Insight-users] Problem with offset calculation in VOI

[Luis Ibanez]

Hi Mikahil,

Thanks for your detailed report.

What version of ITK are you using ?

 From your description, the problem seems to involve mainly
the vtkExtractVOI filter and the export import connection.


Could you help us reduce this to a minimal case ?


Please do the following:

1) make a backup copy of your code   :-)

2) remove the segmentation filters that follow
    the call

      importer->UpdateLargestPossibleRegion();

3) remove the vtk Flip filter.

4) Rerun your test.


and let us knonw if you still get the error message.

If you do, please go to the file

  Insight/Code/Common/itkImageBase.cxx

and in line: 351, just at the end of the method
VerifyRequestedRegion(), add the following:

if( retval == false )
{
std::cout << this->GetRequestedRegion() << std::endl;
std::cout << this->GetLargestPossibleRegion() << std::endl;
}

Then run the test again, and post to the list the output
of the two cout statements above.



     Please let us know what you find.


        Thanks



           Luis


------------------------
Mikhail Kirillov wrote:
> Hi fellas!
> I'm working on my software graduation project at the University of 
> Wollongong, and I have to present it on a trade show in 2 weeks.
> This problem has been giving me headaches for the last 2 weeks, and so 
> far I only got a possible cause, but have no idea how to fix it.
> I have a volume in DICOM format (768x768x75) which I acquire using ITK 
> classes.
> When the volume is imported from hard disk, I pipeline it to VTK in 
> order to display and apply some user interaction (such as selection of 
> ROI, placement of a seed pixel for segmentation, etc.).
> I also extract the volume of interest using vtkExtractVOI filter 
> (itk::ExtractImageFilter and itk::CropImageFilter work as well, but they 
> don't solve my problem either).
>  
> When I run the application, select a region, place a seed pixel and then 
> try to apply segmentation (threshold level set), I get the following 
> exception message:
>  
> *itk::InvalidRequestedRegionError (016DF528)
> Location: "void __thiscall 
> itk::DataObject::PropagateRequestedRegion(void) throw (class 
> itk::InvalidRequestedRegionError)"
> File: ..\..\..\InsightToolkit-3.4.0\Code\Common\itkDataObject.cxx
> Line: 397
> Description: Requested region is (at least partially) outside the 
> largest possible region.*
>  
> It seems that one of the filters involved in segmentation applies the 
> coordinates to a global space, not a local coordinate system of the 
> extracted volume. BTW when I don't apply extaction of the VOI, 
> segmentation and visualization works like a charm.
> Is there a way to fix it?
> Please check out part of my code below and tell me if I'm missing 
> something or doing anything wrong. I might have some extra unnecessary 
> function calls in there (such as excessive use of 
> UpdateLargestPossibleRegion(), I've added them for debugging purposes to 
> compare the image status after it is being passed to each filter):
>  
> // Extraction of VOI
> vtkExtractVOI *voiExtractor = vtkExtractVOI::New();
> voiExtractor->SetVOI(minX, maxX, minY, maxY, minZ, maxZ);  // the values 
> are acquired through user interaction, but even hard-coded ones don't 
> solve the problem
> voiExtractor->SetSampleRate(1,1,1);
> voiExtractor->SetInput(stencil->GetOutput());
> voiExtractor->ReleaseDataFlagOff();
>  
> // Flipthe image back so it is interpreted correctly by ITK
>  vtkImageFlip *flipper = vtkImageFlip::New();
>  flipper->SetInput(voiExtractor->GetOutput());
>  flipper->SetInformationInput(voiExtractor->GetOutput());
>  flipper->SetFilteredAxis(1);
>  flipper->ReleaseDataFlagOff();
>  flipper->UpdateWholeExtent();
>  
> // Connecting VTK and ITK pipelines
>  vtkImageExport *exporter = vtkImageExport::New();
>  exporter->SetInput(flipper->GetOutput());
>  exporter->UpdateWholeExtent();
>  
>  typedef itk::Image<signed short, 3> ImportImageType;
>  typedef itk::VTKImageImport<ImportImageType> ImporterType;
>  ImporterType::Pointer importer = ImporterType::New();
>  
>  ConnectPipelines(exporter, importer);
>  importer->UpdateLargestPossibleRegion();
>  
> // Casting pixel type to float
>  typedef itk::Image<float, 3> InputImageType;
>  typedef itk::CastImageFilter<ImportImageType, InputImageType> 
> ShortToFloatFilter;
>  
>  ImportImageType::Pointer inputImage = importer->GetOutput();
>  
>  ShortToFloatFilter::Pointer caster = ShortToFloatFilter::New();
>  caster->SetInput(inputImage);
>  caster->UpdateLargestPossibleRegion();
>  
> // This part is mostly copy/paste from 
> ThresholdSegmentationLevelSetImageFilter.cxx example
>  typedef itk::Image< signed short, 3 > OutputImageType;
>  typedef itk::BinaryThresholdImageFilter<InputImageType, 
> OutputImageType> ThresholdingFilterType;
>  ThresholdingFilterType::Pointer thresholder = 
> ThresholdingFilterType::New();
>  thresholder->SetLowerThreshold( -1000.0 );
>  thresholder->SetUpperThreshold(     0.0 );
>  thresholder->SetOutsideValue(  0  );
>  thresholder->SetInsideValue(  255 );
>  
>  typedef  itk::FastMarchingImageFilter< InputImageType, InputImageType > 
> FastMarchingFilterType;
>  FastMarchingFilterType::Pointer  fastMarching = 
> FastMarchingFilterType::New();
>  
>  typedef  itk::ThresholdSegmentationLevelSetImageFilter< InputImageType, 
> InputImageType > ThresholdSegmentationLevelSetImageFilterType;
>  ThresholdSegmentationLevelSetImageFilterType::Pointer 
> thresholdSegmentation = ThresholdSegmentationLevelSetImageFilterType::New();
> 
>  thresholdSegmentation->SetPropagationScaling( 1.0 );
>  thresholdSegmentation->SetCurvatureScaling( 1.0 );
>  thresholdSegmentation->SetMaximumRMSError( 0.02 );
>  thresholdSegmentation->SetNumberOfIterations( 1200 );
>  thresholdSegmentation->SetUpperThreshold( 800);
>  thresholdSegmentation->SetLowerThreshold( -600);
>  thresholdSegmentation->SetIsoSurfaceValue(0.0);
>  thresholdSegmentation->SetInput( fastMarching->GetOutput() );
>  thresholdSegmentation->SetFeatureImage( caster->GetOutput() );
> 
> thresholder->SetInput( thresholdSegmentation->GetOutput() );
> 
>  typedef FastMarchingFilterType::NodeContainer NodeContainer;
>  typedef FastMarchingFilterType::NodeType NodeType;
> 
>  NodeContainer::Pointer seeds = NodeContainer::New();
>  
>  InputImageType::IndexType  seedPosition;
>  seedPosition[0] = seedX;
>  seedPosition[1] = seedY;
>  seedPosition[2] = seedZ;
>  
>  const double initialDistance = 5.0;
> 
>  const double seedValue = - initialDistance;
>  
>  NodeType node;
>  node.SetValue( seedValue );
>  node.SetIndex( seedPosition );
>  
>  seeds->Initialize();
>  seeds->InsertElement( 0, node );
>  
>  fastMarching->SetTrialPoints(  seeds  );
>  fastMarching->SetSpeedConstant( 1.0 );
>  
> try
> {
>   fastMarching->SetOutputSize(caster->GetOutput()->GetBufferedRegion().GetSize() 
> );
>   thresholder->UpdateLargestPossibleRegion();   // This is where I get 
> my exception caught
> }
> catch( itk::ExceptionObject & excep )
> {
>    std::cerr << "Exception caught !" << std::endl;
>    std::cerr << excep << std::endl;
> }
>  
> Thanks in advance guys.
> If you need more info on the problem, please let me know.
> -- 
> "Then it comes to be that the soothing light at the end of your tunnel 
> is just a freight train coming your way..." No Leaf Clover, Metallica
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users