[Insight-users] Cannot run the Resampling filter
motes motes
mort.motes at gmail.com
Sun Nov 29 16:58:59 EST 2009
Yes I should have send the stack trace as the first thing, lesson learned :-)
On Sun, Nov 29, 2009 at 10:36 PM, Bill Lorensen <bill.lorensen at gmail.com> wrote:
> 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