[Insight-users] Cannot run the Resampling filter

Bill Lorensen bill.lorensen at gmail.com
Sat Nov 28 08:19:11 EST 2009


The valgrind output
==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.

shows 8,529,290,300 of memort allocated. This seems very large. What
is using so much memory?

Bill

On Sat, Nov 28, 2009 at 5:14 AM, motes motes <mort.motes at gmail.com> wrote:
> 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
>>>>
>>>
>>
> _____________________________________
> 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