[Insight-users] Cannot run the Resampling filter

motes motes mort.motes at gmail.com
Sat Nov 28 05:14:21 EST 2009


Ok I now tried to run Valgrind on the application that computes the
registered image and here is the result:


valgrind --leak-check=yes --track-origins=yes ./MyApp
INFO:: Computing registered image
==8389==
==8389== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 8 from 1)
==8389== malloc/free: in use at exit: 0 bytes in 0 blocks.
==8389== malloc/free: 13,153,309 allocs, 13,153,309 frees,
8,529,290,300 bytes allocated.
==8389== For counts of detected errors, rerun with: -v
==8389== All heap blocks were freed -- no leaks are possible.


So Valgrind cannot find any errors but I still get a segmentation
fault when running the app without Valgrind.

./MyApp
INFO:: Computing registered image
Segmentation fault


I have stripped down my transform to only making the call to the
kd-tree and just pass the incoming points as output points.


// Transform a point
template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
void
AdaptiveBSplineDeformableTransform<TScalarType, NDimensions, VSplineOrder>
::TransformPoint(
  const InputPointType & point,
  OutputPointType & outputPoint,
  WeightsType & weights,
  ParameterIndexArrayType & indices,
  bool & inside ) const
{
  unsigned int j;
  IndexType supportIndex;
  inside = true;
  InputPointType transformedPoint;

  if ( m_BulkTransform ) {
    transformedPoint = m_BulkTransform->TransformPoint( point );
  } else {
    transformedPoint = point;
  }

  outputPoint.Fill( NumericTraits<ScalarType>::Zero );
  if ( m_CoefficientImage[0] ) {
      outputPoint.Fill( NumericTraits<ScalarType>::Zero );

     VectorType vectorPoint;
     vectorPoint.Fill(0);

    for (int i=0; i<NDimensions; i++) {
      vectorPoint[i] = point[i];
    }
    NeighborsType neighbors;
    tree->Search(vectorPoint, 5.0, neighbors);

...
...




I have also written the following unit-test of the transform where the test

   BOOST_AUTO_TEST_CASE(test_resampler)

still give a segmentation error:



#include "itkArray.h"

#include "itkVector.h"

#include "itkImage.h"

#include "itkImageFileReader.h"

#include "My_BSpline_deformable_transform.h"

#include "My_create_regular_grid.h"

#include "project_paths.h"

#include "display_image.h"



#define BOOST_AUTO_TEST_MAIN

#include <boost/test/auto_unit_test.hpp>



// Boost Test declaration and Checking macros

#include <boost/test/unit_test_suite.hpp>

#include <boost/test/test_tools.hpp>

#include <boost/test/floating_point_comparison.hpp>







static const unsigned int Dimension = 3;

static const unsigned int SpaceDimension = Dimension;

static const unsigned int SplineOrder = 3;

typedef float																				
        PixelType;

typedef double																									
  CoordinateRepType;

typedef itk::Image<PixelType, Dimension>					
                      ImageType;

typedef itk::ImageFileReader<ImageType>				
                        ImageReaderType;

typedef itk::MyBSplineDeformableTransform<CoordinateRepType,
SpaceDimension, SplineOrder >  TransformType;

typedef TransformType::RegionType
                           RegionType;

typedef TransformType::ControlPointContainerType
                           ControlPointContainerType;

typedef TransformType::ControlPointType
                           ControlPointType;

typedef TransformType::ControlPointIteratorType
                           ControlPointIteratorType;

typedef TransformType::InputPointType
                           InputPointType;

typedef TransformType::ParametersType
                           ParametersType;

typedef TransformType::VectorType
                           VectorType;

typedef RegionType::IndexType
                           IndexType;

typedef RegionType::SizeType
                           SizeType;

typedef itk::ResampleImageFilter< FixedImageType, FixedImageType >
                           ResampleFilterType;

typedef itk::LinearInterpolateImageFunction< FixedImageType, double >
                           LinearInterpolatorType;



std::string fixedImagePath = base_dir +  "/local/images/fixed/fix.dcm";

std::string movingImagePath = base_dir +  "/local/images/moving/mov.dcm";

FixedImageType::Pointer  imageF = FixedImageType::New();

FixedImageType::Pointer  imageM = FixedImageType::New();

FixedImageType::Pointer  imageR = FixedImageType::New();

FileManager filemanager;





struct LoadImages {

  LoadImages() {

    std::cout << "Loading images\n";

    filemanager.LoadFile<CImg_Image,
FixedImageType::Pointer>(fixedImagePath, imageF);

    filemanager.LoadFile<CImg_Image,
FixedImageType::Pointer>(movingImagePath,imageM);



  }



  ~LoadImages() {

    std::cout << "global teardown\n";

  }

};





BOOST_GLOBAL_FIXTURE( LoadImages );





BOOST_AUTO_TEST_SUITE(apply_transform);

BOOST_AUTO_TEST_SUITE();





BOOST_AUTO_TEST_CASE(test_image_data)

{

  double tolerance = 0.00001;

  // Image dimension

  FixedImageType::SizeType imageSize =
imageF->GetLargestPossibleRegion().GetSize();

	BOOST_CHECK_EQUAL(imageSize[0],128);

	BOOST_CHECK_EQUAL(imageSize[1],128);

	BOOST_CHECK_EQUAL(imageSize[2],88);



  // Image spacing

  VectorType imageSpacing = imageF->GetSpacing();

  BOOST_CHECK_CLOSE(imageSpacing[0],3.125, tolerance);

  BOOST_CHECK_CLOSE(imageSpacing[1],3.125, tolerance);

  BOOST_CHECK_CLOSE(imageSpacing[2],4.0, tolerance);



  // Image origin

  InputPointType imageOrigin = imageF->GetOrigin();

	BOOST_CHECK_EQUAL(imageOrigin[0],0);

	BOOST_CHECK_EQUAL(imageOrigin[1],0);

	BOOST_CHECK_EQUAL(imageOrigin[2],0);



}





BOOST_AUTO_TEST_CASE(test_index_operator)

{

  int nodes = 4;

  ControlPointContainerType controlPointContainer;

  My::InitRegularControlPointGrid<ControlPointContainerType,
FixedImageType, InputPointType>(nodes, Dimension,
controlPointContainer, imageF);

  // Test the indexing operator. controlPointContainer[i] should
return the control point with id=i.

  int ctrlNum = controlPointContainer.size();

  for (int i=0; i<ctrlNum; i++) {

    ControlPointType cp = controlPointContainer[i];

    BOOST_CHECK_EQUAL(cp.getId(),i);

  }



}







BOOST_AUTO_TEST_CASE(test_physical_coordinates)

{



  /**

   * Create regular 3D cube of 4*4*4= 64 points.

   *

   */

  int nodes = 4;

  ControlPointContainerType controlPointContainer;

  My::InitRegularControlPointGrid<ControlPointContainerType,
FixedImageType, InputPointType>(nodes, Dimension,
controlPointContainer, imageF);

  ControlPointIteratorType it = controlPointContainer.begin();



  /*

  * We expect the control point with ID=63 to have maximum physical coordinates:

  *

  * Max index = 	[127, 127, 87]

  * Spacing = [3.125, 3.125, 4.0]

  * Max physical = [127, 127, 87] * [3.125, 3.125, 4.0] = [396.875,
396.875, 348]

  *

  */

  double pretty_tolerant = 3.0; // There will some rounding errors
when computing the offset between control points.

  ControlPointType cp = controlPointContainer[63];

  BOOST_CHECK_CLOSE(cp.getLocation()[0], 396.875, pretty_tolerant);

  BOOST_CHECK_CLOSE(cp.getLocation()[1], 396.875, pretty_tolerant);

  BOOST_CHECK_CLOSE(cp.getLocation()[2], 348.0, pretty_tolerant);

}





BOOST_AUTO_TEST_CASE(test_resampler)

{

  LinearInterpolatorType::Pointer linearInterpolator =
LinearInterpolatorType::New();

  TransformType::Pointer MyTransform = TransformType::New();



  // Setup the cube

  int nodes = 4;

  ControlPointContainerType controlPointContainer;

  My::InitRegularControlPointGrid<ControlPointContainerType,
FixedImageType, InputPointType>(nodes, Dimension,
controlPointContainer, imageF);

  MyTransform->SetControlPoints(controlPointContainer);

  MyTransform->SetActiveControlPoints(controlPointContainer.size());

  MyTransform->SetGridSpacing(imageF->GetSpacing());

  MyTransform->SetGridOrigin(imageF->GetOrigin());

  MyTransform->BuildKdTree();

  FixedImageType::SizeType fixedImageSize =
imageF->GetLargestPossibleRegion().GetSize();

  MyTransform->ComputeKernelSize(nodes, fixedImageSize);



  // Add parameters

  unsigned int numberOfBSplineParameters = MyTransform->GetNumberOfParameters();
  std::cout << "Number of bspline param = " <<
numberOfBSplineParameters << std::endl;
  // Init all parameters to 0.0.
  ParametersType parameters(numberOfBSplineParameters);
  parameters.Fill(0.0);
  MyTransform->SetParameters(parameters);




  ResampleFilterType::Pointer resampler = ResampleFilterType::New();
  resampler->SetInput(imageM);
  resampler->SetTransform(MyTransform);
  resampler->SetInterpolator(linearInterpolator);
  resampler->SetOutputOrigin(imageF->GetOrigin());
  resampler->SetOutputSpacing(imageF->GetSpacing());
  resampler->SetSize(imageF->GetLargestPossibleRegion().GetSize());
  resampler->SetDefaultPixelValue(0);


  try {

    resampler->Update();

  }

  catch( itk::ExceptionObject & err ) {

    std::cerr << "ExceptionObject caught !" << std::endl;

    std::cerr << err << std::endl;

  }
  CopyImageToImage<typename FixedImageType::Pointer,
FixedImageType>(resampler->GetOutput(), imageR);

  display_image<FixedImageType, CImg_Image>(imageR);









}







BOOST_AUTO_TEST_SUITE_END();



BOOST_AUTO_TEST_SUITE_END();





Any ideas on what migh be causing this segmentation fault is most
welcome (have almost spend a week now trying to find the cause).





On Thu, Nov 26, 2009 at 7:58 PM, motes motes <mort.motes at gmail.com> wrote:
> On Thu, Nov 26, 2009 at 7:42 PM, Luis Ibanez <luis.ibanez at kitware.com> wrote:
>> Hi Motes,
>>
>>
>> 1) The Scope of the "try" block doesn't destroy any of the
>>     components of the registration framework.
>>
>
> Ok that was also my understanding.
>
>
>>
>> 2) I don't see how you are passing any arguments to the
>>    function:
>>
>>                     computeRegisteredImage();
>>
>>     are you declaring them as global ?
>>
>>
>
> Yes the registration components and images are "global" fields in  a
> wrapper class that I have written.
>
>
>> 3)  No offense,
>>      but your custom Transform is the main suspect.
>>
>>     Here is the important question:
>>
>>
>>                    Have you written a Unit Test
>>                    for your new Transform ?
>>
>>
>>     I can anticipate that the answer is "no",    :-)
>>
>>     and therefore
>>
>>     I will strongly,
>>     very, very, very strongly suggest
>>
>>     that you write such unit test and fully debug
>>     your transform before you attempt to use it
>>     as part of a larger program.
>>
>>     For examples on how to test Transforms,
>>     Please look at the files
>>
>>      Insight/Testing/Code/Common/
>>
>> itkAffineTransformTest.cxx
>> itkAzimuthElevationToCartesianTransformTest.cxx
>> itkBSplineDeformableTransformTest2.cxx
>> itkBSplineDeformableTransformTest.cxx
>> itkCenteredAffineTransformTest.cxx
>> itkCenteredEuler3DTransformTest.cxx
>> itkCenteredRigid2DTransformTest.cxx
>> itkCenteredTransformInitializerTest.cxx
>> itkCenteredVersorTransformInitializerTest.cxx
>> itkEuler2DTransformTest.cxx
>> itkEuler3DTransformTest.cxx
>> itkFixedCenterOfRotationAffineTransformTest.cxx
>> itkIdentityTransformTest.cxx
>> itkImageTransformTest.cxx
>> itkLandmarkBasedTransformInitializerTest.cxx
>> itkQuaternionRigidTransformTest.cxx
>> itkRigid2DTransformTest.cxx
>> itkRigid3DPerspectiveTransformTest.cxx
>> itkRigid3DTransformTest.cxx
>> itkRigid3DTransformTest.cxx.orig
>> itkScaleLogarithmicTransformTest.cxx
>> itkScaleSkewVersor3DTransformTest.cxx
>> itkScaleTransformTest.cxx
>> itkSimilarity2DTransformTest.cxx
>> itkSimilarity3DTransformTest.cxx
>> itkSplineKernelTransformTest.cxx
>> itkTransformsSetParametersTest.cxx
>> itkTransformTest.cxx
>> itkTranslationTransformTest.cxx
>> itkVersorRigid3DTransformTest.cxx
>> itkVersorTransformTest.cxx
>>
>>
>>
>
> Actually all new code I write is done through unit-tests. But thanks
> for the above suggestions!
>
>
>
>
>
>
>>   Regards,
>>
>>
>>        Luis
>>
>>
>> ------------------------------------------------------------------------------------------
>> On Wed, Nov 25, 2009 at 8:44 PM, motes motes <mort.motes at gmail.com> wrote:
>>> Which objects survive in the transform component when the image
>>> registration method is finished?
>>>
>>> By finished I mean reaches the scope after this try/catch block:
>>>
>>>    try {
>>>      registration->Update();
>>>    }
>>>    catch( itk::ExceptionObject & err ) {
>>>      std::cerr << "ExceptionObject caught !" << std::endl;
>>>      std::cerr << err << std::endl;
>>>      return EXIT_FAILURE;
>>>    }
>>>
>>> // What lives in the image registration method when arriving here?
>>>
>>>
>>>
>>>
>>>
>>> The reason I ask is because I am trying to compute the registered
>>> image AFTER the registration method is done. Something like this:
>>>
>>>    try {
>>>
>>>      registration->Update();
>>>    }
>>>    catch( itk::ExceptionObject & err ) {
>>>      std::cerr << "ExceptionObject caught !" << std::endl;
>>>      std::cerr << err << std::endl;
>>>      return EXIT_FAILURE;
>>>    }
>>>
>>>    computeRegisteredImage();
>>>    return EXIT_SUCCESS;
>>>
>>>
>>> where:
>>>
>>>  void computeRegisteredImage(){
>>>    ResampleFilterType::Pointer resampler = ResampleFilterType::New();
>>>    resampler->SetInput(imageM);
>>>    resampler->SetTransform(registration->GetTransform());
>>>    resampler->SetInterpolator(registration->GetInterpolator());
>>>    resampler->SetOutputOrigin(imageF->GetOrigin());
>>>    resampler->SetOutputSpacing(imageF->GetSpacing());
>>>    resampler->SetSize(imageF->GetLargestPossibleRegion().GetSize());
>>>    resampler->SetDefaultPixelValue(default_pixel_value);
>>>    try {
>>>      resampler->Update();
>>>    }
>>>    catch( itk::ExceptionObject & err ) {
>>>      std::cerr << "ExceptionObject caught !" << std::endl;
>>>      std::cerr << err << std::endl;
>>>    }
>>>
>>>
>>>    imageR = resampler->GetOutput();
>>>
>>>
>>>  }
>>>
>>>
>>>
>>> But:
>>>
>>> resampler->Update();
>>>
>>> never returns. It thows a segmentation fault. I have located the error
>>> to be in the TransformPoint method. I have made my own transform which
>>> uses a kdTree in the TransformPoint function. This kdTree is stored
>>> locally in my transform like:
>>>
>>>
>>> #include "kd_tree.h" // My own itk:kdTree wrapper.
>>> namespace itk
>>> {
>>> template <
>>>    class TScalarType = double,          // Data type for scalars
>>>    unsigned int NDimensions = 3,        // Number of dimensions
>>>    unsigned int VSplineOrder = 3 >      // Spline order
>>> class ITK_EXPORT MyDeformableTransform :
>>>          public Transform< TScalarType, NDimensions, NDimensions >
>>> {
>>> public:
>>>  /** Standard class typedefs. */
>>>  typedef MyDeformableTransform                                          Self;
>>> ...
>>> ...
>>> ...
>>>
>>>  typedef KdTree<VectorType, ControlPointContainerType, NDimensions>
>>>               KdTreeType;
>>>
>>>
>>> ....
>>> .....
>>> ...
>>> protected:
>>>  mutable KdTreeType m_kdTree;
>>>
>>>
>>>
>>>
>>>
>>> Can this have something to do with the segmentation error? Does it
>>> need to be stored as a smart pointer instead?
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>> In the TransformPoint function I do:
>>>
>>> // Transform a point
>>> template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
>>> void
>>> MyDeformableTransform<TScalarType, NDimensions, VSplineOrder>
>>> ::TransformPoint(
>>>  const InputPointType & point,
>>>  OutputPointType & outputPoint,
>>>  WeightsType & weights,
>>>  ParameterIndexArrayType & indices,
>>>  bool & inside ) const
>>> {
>>>  unsigned int j;
>>>  IndexType supportIndex;
>>>  inside = true;
>>>  InputPointType transformedPoint;
>>>
>>>  if ( m_BulkTransform ) {
>>>    transformedPoint = m_BulkTransform->TransformPoint( point );
>>>  } else {
>>>    transformedPoint = point;
>>>  }
>>>
>>>  outputPoint.Fill( NumericTraits<ScalarType>::Zero );
>>>  if ( GetNumberOfParameters()>0) {
>>>    outputPoint.Fill( NumericTraits<ScalarType>::Zero );
>>>
>>>
>>>
>>>     NeighborContainerType neighbors;
>>>
>>>      // This results in a segmentation error
>>>      this->FindNeighbors(point, m_kernel_radius, neighbors);
>>>
>>>
>>>
>>> its the call:
>>>
>>>      this->FindNeighbors(point, m_kernel_radius, neighbors);
>>>
>>> that gives the error which basically just calls the itk::kdTree class.
>>>
>>> Any suggestions are most welcome!
>>> _____________________________________
>>> Powered by www.kitware.com
>>>
>>> Visit other Kitware open-source projects at
>>> http://www.kitware.com/opensource/opensource.html
>>>
>>> Kitware offers ITK Training Courses, for more information visit:
>>> http://www.kitware.com/products/protraining.html
>>>
>>> Please keep messages on-topic and check the ITK FAQ at:
>>> http://www.itk.org/Wiki/ITK_FAQ
>>>
>>> Follow this link to subscribe/unsubscribe:
>>> http://www.itk.org/mailman/listinfo/insight-users
>>>
>>
>


More information about the Insight-users mailing list