[Insight-users] Cannot run the Resampling filter

motes motes mort.motes at gmail.com
Sun Nov 29 15:58:15 EST 2009


You are a wizard, it now works! Damn over a week spend on this problem
and then a single line fixes the whole thing, but this is great :-)

Next thing to figure out is how to set the:

itk::MultiThreader::SetGlobalMaximumNumberOfThreads(1);

only when using the resampler and not when performing the registration
(the kdtree still works fine in the registration process). I assume
that some performance is lost when setting the maximum number of
thread to 1?










On Sun, Nov 29, 2009 at 5:29 PM, Bill Lorensen <bill.lorensen at gmail.com> wrote:
> In your program that uses the kdtree in the transform, add this line
> at the beginning of the main program:
>
> itk::MultiThreader::SetGlobalMaximumNumberOfThreads(1);
>
> This will tell us if the kdtree has problems in multiple threads.
>
> On Sat, Nov 28, 2009 at 8:49 AM, motes motes <mort.motes at gmail.com> wrote:
>> Btw: sometimes I also get this error:
>>
>> ExceptionObject caught !
>>
>> itk::ExceptionObject (00A4F634)
>> Location: "void __thiscall itk::MultiThreader::SingleMethodExecute(void)"
>> File: ..\..\..\Code\Common\itkMultiThreader.cxx
>> Line: 471
>> Description: itk::ERROR: MultiThreader(02D36B18): Exception occurred
>> during SingleMethodExecute
>>
>> The SingleMethodExecute is called just before calling the transform in
>> the resampling filter. So this is actually an error that appears
>> before the transform is called.
>>
>>
>>
>> On Sat, Nov 28, 2009 at 2:47 PM, motes motes <mort.motes at gmail.com> wrote:
>>> I am loading 3D dicom images. The original size is 256*256*177. But
>>> when I load them I subsample them to 128*128*88. I am using a machine
>>> with 6GB RAM so it should not be a problem.
>>>
>>>
>>>
>>> On Sat, Nov 28, 2009 at 2:19 PM, Bill Lorensen <bill.lorensen at gmail.com> wrote:
>>>> 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