[Insight-users] Problem with offset calculation in VOI

[Mikhail Kirillov]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20080420/202de02e/attachment-0001.htm>