[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