[Insight-users] Cast & Deformable Filter

[Luis Ibanez]

This is a multi-part message in MIME format.
--------------000303000001000800030108
Content-Type: text/plain; charset=us-ascii; format=flowed
Content-Transfer-Encoding: 7bit


Hi Pieter,

I have been trying to reproduce the
problem that you reported concerning
the interaction between the itkCastImageFilter
and the DeformableMesh3DFilter.

Sorry that it took some time to start doing this.

 From the code in your message I reconstructed
the attached source file which compiles without
problem.  

However it is seg. faulting at write.

Could you please look at the parameters
used in this file and tune them to values that
you consider appropiated.

I doubt that the problem is not comming from
the itkCastImageFilter. This is one of the simplest
filters in the toolkit and use the same structure of
filters like "add" and "multiply".

The only suspicious thing I found in your code
is that you are calling "SetOutputImage()" on
the Deformable filter.  This call is deprecated.

ITK filters allocate and produce their output.

You should only to get the image from the filter.

The fact that your code is running one iteration
before producing nan's could lead to think that
memory is being corrupted during the execution
of the filter.

Could you please take a look at the attached file
and verify that it reflects correctly the set of operations
that you are doing in your code. If you can, please
select appropiated valued for the "const" parameters
at the top of  main() and send it back to me.

Once we reproduce the problem, we can look for
the reasons behind and hopefully find a solution.



  Thanks


     Luis






--------------000303000001000800030108
Content-Type: text/plain;
 name="CMakeLists.txt"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="CMakeLists.txt"

PROJECT(PieterProject)

# 
# Find ITK
#
FIND_PATH( ITK_BINARY_DIR itkConfigure.h )

# Load in the values from ITK if found
IF ( ITK_BINARY_DIR )
  LOAD_CACHE(${ITK_BINARY_DIR})
  INCLUDE (${ITK_SOURCE_DIR}/itkCMakeOptions.cmake)
ENDIF (ITK_BINARY_DIR )


INCLUDE_DIRECTORIES(
  ${ITK_SOURCE_DIR}/Code/Common
  ${ITK_SOURCE_DIR}/Code/BasicFilters
  ${ITK_SOURCE_DIR}/Code/Algorithms
  ${ITK_SOURCE_DIR}/Code/IO
)


LINK_DIRECTORIES(
${ITK_BINARY_DIR}/Code/Common
${ITK_BINARY_DIR}/Code/BasicFilters
${ITK_BINARY_DIR}/Code/Algorithms
${ITK_BINARY_DIR}/Code/IO
)


LINK_LIBRARIES (
  VXLNumerics
  ITKCommon
  ITKIO
)


ADD_EXECUTABLE(CastTest CastTest )


--------------000303000001000800030108
Content-Type: text/plain;
 name="CastTest.cpp"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="CastTest.cpp"


#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkVTKImageIO.h"
#include "itkCastImageFilter.h"
#include "itkDeformableMesh3DFilter.h"
#include "itkSphereMeshSource.h"
#include "itkMesh.h"
#include "itkGradientRecursiveGaussianImageFilter.h"



int main()
{

  const float SEEDX = 20.0;
  const float SEEDY = 20.0;
  const float SEEDZ = 20.0;

  const int   m_radius        =   20;
  const float m_stiffness_h   = 1.0f;
  const float m_stiffness_v   = 1.0f;
  const float m_stimestep     = 1.0f;
  const float m_zdistance     = 1.0f;

  const int   m_xresolution  = 1;
  const int   m_yresolution  = 1;
  const int   m_zresolution  = 1;

  const int   m_stepthreshold = 10;
  const float m_slicedistancethreshold = 1.0f;
  const float m_modeldistancetoboundarythreshold = 1.0f;

  typedef itk::Image<float,3>             inputImageType;
  typedef itk::Image<unsigned short,3>    potImageType;

  typedef itk::ImageFileReader<inputImageType>  ReaderType;
  typedef itk::ImageFileWriter<potImageType>    WriterType;

  typedef itk::Point<float,3>               OPointType;

  typedef itk::CastImageFilter< inputImageType, potImageType > myCasterType;

  typedef itk::Mesh<float> FMesh;

  typedef itk::DeformableMesh3DFilter<FMesh, FMesh> DFilter;

  typedef DFilter::GradientImageType          GradientImageType;

  typedef potImageType::IndexType myIndexType;

  typedef itk::GradientRecursiveGaussianImageFilter< 
                                          potImageType,
                                          GradientImageType > GradientFilterType;

	ReaderType::Pointer potImgReader = 
                itk::ImageFileReader<inputImageType>::New();

  itk::VTKImageIO::Pointer vtkIO = itk::VTKImageIO::New();
	potImgReader->SetImageIO(vtkIO);
	potImgReader->SetFileName("sphere.vtk");	

	try {
		potImgReader->Update();
	} catch (itk::ExceptionObject & e) {
	    std::cerr << "exception in file potential volume " << std::endl;
	    std::cerr << e.GetDescription() << std::endl;
	    std::cerr << e.GetLocation() << std::endl;
	    return EXIT_FAILURE;
	}
	
	potImageType::Pointer potImg = potImageType::New();
  
  typedef potImageType::RegionType   potImageRegionType;

  	//Get the image setting from the reader
	potImageRegionType potRegion = 
        potImgReader->GetOutput()->GetLargestPossibleRegion();

	potRegion.SetIndex( potImageRegionType::IndexType::ZeroIndex );
	potImg->SetLargestPossibleRegion( potRegion );
	potImg->SetBufferedRegion( potRegion );
	potImg->SetRequestedRegion( potRegion );
	potImg->SetOrigin(  potImgReader->GetOutput()->GetOrigin() );
	potImg->SetSpacing( potImgReader->GetOutput()->GetSpacing() );
	potImg->Allocate();
	
//	potImg->Print(std::cout);

  /* Testing purpose: */
	itk::ImageRegionIteratorWithIndex <potImageType> ptit(potImg, potRegion);
	itk::ImageRegionIteratorWithIndex<inputImageType> readerit(potImgReader->GetOutput(), potRegion);
	ptit.GoToBegin();
	readerit.GoToBegin();
	 
	while( !ptit.IsAtEnd()) {
		ptit.Set((unsigned short) readerit.Get());
		++ptit;
		++readerit;
	} 

	myCasterType::Pointer potCaster = myCasterType::New();
	potCaster->SetInput( potImgReader->GetOutput() );
	std::cout << "update"<<std::endl;
	try{
		potCaster->Update();
	  } 
  catch (itk::ExceptionObject & e) 
    {
	    std::cerr << "exception in file potential volume " << std::endl;
	    std::cerr << e.GetDescription() << std::endl;
	    std::cerr << e.GetLocation() << std::endl;
	    return EXIT_FAILURE;
  	}
//	potCaster->Print(std::cout);

	WriterType::Pointer potImgWriter = WriterType::New();
	potImgWriter->SetImageIO(vtkIO);
	potImgWriter->SetFileName("cast.vtk");
	potImgWriter->SetInput( potCaster->GetOutput() );
	std::cout <<"Write: " <<std::endl;
	try 
    {
		potImgWriter->Write();
  	}
   catch (itk::ExceptionObject & e) {
	    std::cerr << "exception in file potential volume " << std::endl;
	    std::cerr << e.GetDescription() << std::endl;
	    std::cerr << e.GetLocation() << std::endl;
	    return EXIT_FAILURE;
	}


  DFilter::Pointer dfilter = DFilter::New();
  
//----------------spheresource for testing -------------
  typedef itk::SphereMeshSource<FMesh> SphereSourceType;
  SphereSourceType::Pointer spheresource = SphereSourceType::New();
  dfilter->SetInput(spheresource->GetOutput());

  OPointType spherecenter;
  OPointType scale;
  spherecenter[0] = SEEDX;
  spherecenter[1] = SEEDY;
  spherecenter[2] = SEEDZ;
  scale[0] = 5;
  scale[1] = 5;
  scale[2] = 5;
  spheresource->SetCenter(spherecenter);
  spheresource->SetResolutionX(20);/*Moet zelfde zijn als my_res ???*/
  spheresource->SetResolutionY(20);
  spheresource->SetSquareness1(0.5);
  spheresource->SetSquareness2(0.5);
  spheresource->SetScale(scale);
  spheresource->Update();

//  PutMeshInFile(spheresource->GetOutput(), "./Images/spheresource.vtk");
//---------------- spheresource for testing -------------
  dfilter->SetInput(  spheresource->GetOutput() );

  //Setting for the user for finetuning deformable filter
  myIndexType center;
  center[0] = static_cast<unsigned long>(SEEDX);
  center[1] = static_cast<unsigned long>(SEEDY);
  center[2] = static_cast<unsigned long>(SEEDZ);

  dfilter->SetCenter(center);
  dfilter->SetNeighborRadius(m_radius);
  dfilter->SetStiffnessV(m_stiffness_v);
  dfilter->SetStiffnessH(m_stiffness_h);
  dfilter->SetTimeStep(m_stimestep);
  dfilter->SetZDistance(m_zdistance);
  
  GradientFilterType::Pointer gradFilter = GradientFilterType::New();
  gradFilter->SetInput( potCaster->GetOutput() );

  //Connect the input images and initial model
  //dfilter->SetInput( initModel );
  dfilter->SetPotential( potImg );//potCaster->GetOutput() );	//Potential mapping
  dfilter->SetGradient( gradFilter->GetOutput() );//second derative gradient 
//  dfilter->SetImageOutput(staticsImg);			//Statistics
  dfilter->SetXResolution(m_xresolution);
  dfilter->SetYResolution(m_yresolution);
  dfilter->SetZResolution(m_zresolution);
  dfilter->SetStepThreshold(m_stepthreshold);
  dfilter->SetSliceDistanceThreshold(m_slicedistancethreshold);

dfilter->SetModelDistanceToBoundaryThreshold(m_modeldistancetoboundarythreshold);

 dfilter->Update();


return 0;
}




--------------000303000001000800030108--