[Insight-users] Cannot run the Resampling filter
Bill Lorensen
bill.lorensen at gmail.com
Sun Nov 29 16:36:07 EST 2009
On any filter, you can filter->SetNumberOfThreads(1);
See if that works.
Even though we fixed the problem, we should look at the possibility of
making itk::KdTree thread safe.
Also,
I'm no wizard, after all of this time you sent the full stack trace
that gave me the clue. In fact, after over 25 e-mails, multiple e-mail
threads, that one e-mail (suggested by Karthik in another thread) was
the magic.
I hope we all learn something from this,
Bill
On Sun, Nov 29, 2009 at 3:58 PM, motes motes <mort.motes at gmail.com> wrote:
> 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