[Insight-users] Cannot run the Resampling filter
Bill Lorensen
bill.lorensen at gmail.com
Sat Nov 28 08:19:11 EST 2009
The valgrind output
==8389==
==8389== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 8 from 1)
==8389== malloc/free: in use at exit: 0 bytes in 0 blocks.
==8389== malloc/free: 13,153,309 allocs, 13,153,309 frees,
8,529,290,300 bytes allocated.
==8389== For counts of detected errors, rerun with: -v
==8389== All heap blocks were freed -- no leaks are possible.
shows 8,529,290,300 of memort allocated. This seems very large. What
is using so much memory?
Bill
On Sat, Nov 28, 2009 at 5:14 AM, motes motes <mort.motes at gmail.com> wrote:
> Ok I now tried to run Valgrind on the application that computes the
> registered image and here is the result:
>
>
> valgrind --leak-check=yes --track-origins=yes ./MyApp
> INFO:: Computing registered image
> ==8389==
> ==8389== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 8 from 1)
> ==8389== malloc/free: in use at exit: 0 bytes in 0 blocks.
> ==8389== malloc/free: 13,153,309 allocs, 13,153,309 frees,
> 8,529,290,300 bytes allocated.
> ==8389== For counts of detected errors, rerun with: -v
> ==8389== All heap blocks were freed -- no leaks are possible.
>
>
> So Valgrind cannot find any errors but I still get a segmentation
> fault when running the app without Valgrind.
>
> ./MyApp
> INFO:: Computing registered image
> Segmentation fault
>
>
> I have stripped down my transform to only making the call to the
> kd-tree and just pass the incoming points as output points.
>
>
> // Transform a point
> template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
> void
> AdaptiveBSplineDeformableTransform<TScalarType, NDimensions, VSplineOrder>
> ::TransformPoint(
> const InputPointType & point,
> OutputPointType & outputPoint,
> WeightsType & weights,
> ParameterIndexArrayType & indices,
> bool & inside ) const
> {
> unsigned int j;
> IndexType supportIndex;
> inside = true;
> InputPointType transformedPoint;
>
> if ( m_BulkTransform ) {
> transformedPoint = m_BulkTransform->TransformPoint( point );
> } else {
> transformedPoint = point;
> }
>
> outputPoint.Fill( NumericTraits<ScalarType>::Zero );
> if ( m_CoefficientImage[0] ) {
> outputPoint.Fill( NumericTraits<ScalarType>::Zero );
>
> VectorType vectorPoint;
> vectorPoint.Fill(0);
>
> for (int i=0; i<NDimensions; i++) {
> vectorPoint[i] = point[i];
> }
> NeighborsType neighbors;
> tree->Search(vectorPoint, 5.0, neighbors);
>
> ...
> ...
>
>
>
>
> I have also written the following unit-test of the transform where the test
>
> BOOST_AUTO_TEST_CASE(test_resampler)
>
> still give a segmentation error:
>
>
>
> #include "itkArray.h"
>
> #include "itkVector.h"
>
> #include "itkImage.h"
>
> #include "itkImageFileReader.h"
>
> #include "My_BSpline_deformable_transform.h"
>
> #include "My_create_regular_grid.h"
>
> #include "project_paths.h"
>
> #include "display_image.h"
>
>
>
> #define BOOST_AUTO_TEST_MAIN
>
> #include <boost/test/auto_unit_test.hpp>
>
>
>
> // Boost Test declaration and Checking macros
>
> #include <boost/test/unit_test_suite.hpp>
>
> #include <boost/test/test_tools.hpp>
>
> #include <boost/test/floating_point_comparison.hpp>
>
>
>
>
>
>
>
> static const unsigned int Dimension = 3;
>
> static const unsigned int SpaceDimension = Dimension;
>
> static const unsigned int SplineOrder = 3;
>
> typedef float
> PixelType;
>
> typedef double
> CoordinateRepType;
>
> typedef itk::Image<PixelType, Dimension>
> ImageType;
>
> typedef itk::ImageFileReader<ImageType>
> ImageReaderType;
>
> typedef itk::MyBSplineDeformableTransform<CoordinateRepType,
> SpaceDimension, SplineOrder > TransformType;
>
> typedef TransformType::RegionType
> RegionType;
>
> typedef TransformType::ControlPointContainerType
> ControlPointContainerType;
>
> typedef TransformType::ControlPointType
> ControlPointType;
>
> typedef TransformType::ControlPointIteratorType
> ControlPointIteratorType;
>
> typedef TransformType::InputPointType
> InputPointType;
>
> typedef TransformType::ParametersType
> ParametersType;
>
> typedef TransformType::VectorType
> VectorType;
>
> typedef RegionType::IndexType
> IndexType;
>
> typedef RegionType::SizeType
> SizeType;
>
> typedef itk::ResampleImageFilter< FixedImageType, FixedImageType >
> ResampleFilterType;
>
> typedef itk::LinearInterpolateImageFunction< FixedImageType, double >
> LinearInterpolatorType;
>
>
>
> std::string fixedImagePath = base_dir + "/local/images/fixed/fix.dcm";
>
> std::string movingImagePath = base_dir + "/local/images/moving/mov.dcm";
>
> FixedImageType::Pointer imageF = FixedImageType::New();
>
> FixedImageType::Pointer imageM = FixedImageType::New();
>
> FixedImageType::Pointer imageR = FixedImageType::New();
>
> FileManager filemanager;
>
>
>
>
>
> struct LoadImages {
>
> LoadImages() {
>
> std::cout << "Loading images\n";
>
> filemanager.LoadFile<CImg_Image,
> FixedImageType::Pointer>(fixedImagePath, imageF);
>
> filemanager.LoadFile<CImg_Image,
> FixedImageType::Pointer>(movingImagePath,imageM);
>
>
>
> }
>
>
>
> ~LoadImages() {
>
> std::cout << "global teardown\n";
>
> }
>
> };
>
>
>
>
>
> BOOST_GLOBAL_FIXTURE( LoadImages );
>
>
>
>
>
> BOOST_AUTO_TEST_SUITE(apply_transform);
>
> BOOST_AUTO_TEST_SUITE();
>
>
>
>
>
> BOOST_AUTO_TEST_CASE(test_image_data)
>
> {
>
> double tolerance = 0.00001;
>
> // Image dimension
>
> FixedImageType::SizeType imageSize =
> imageF->GetLargestPossibleRegion().GetSize();
>
> BOOST_CHECK_EQUAL(imageSize[0],128);
>
> BOOST_CHECK_EQUAL(imageSize[1],128);
>
> BOOST_CHECK_EQUAL(imageSize[2],88);
>
>
>
> // Image spacing
>
> VectorType imageSpacing = imageF->GetSpacing();
>
> BOOST_CHECK_CLOSE(imageSpacing[0],3.125, tolerance);
>
> BOOST_CHECK_CLOSE(imageSpacing[1],3.125, tolerance);
>
> BOOST_CHECK_CLOSE(imageSpacing[2],4.0, tolerance);
>
>
>
> // Image origin
>
> InputPointType imageOrigin = imageF->GetOrigin();
>
> BOOST_CHECK_EQUAL(imageOrigin[0],0);
>
> BOOST_CHECK_EQUAL(imageOrigin[1],0);
>
> BOOST_CHECK_EQUAL(imageOrigin[2],0);
>
>
>
> }
>
>
>
>
>
> BOOST_AUTO_TEST_CASE(test_index_operator)
>
> {
>
> int nodes = 4;
>
> ControlPointContainerType controlPointContainer;
>
> My::InitRegularControlPointGrid<ControlPointContainerType,
> FixedImageType, InputPointType>(nodes, Dimension,
> controlPointContainer, imageF);
>
> // Test the indexing operator. controlPointContainer[i] should
> return the control point with id=i.
>
> int ctrlNum = controlPointContainer.size();
>
> for (int i=0; i<ctrlNum; i++) {
>
> ControlPointType cp = controlPointContainer[i];
>
> BOOST_CHECK_EQUAL(cp.getId(),i);
>
> }
>
>
>
> }
>
>
>
>
>
>
>
> BOOST_AUTO_TEST_CASE(test_physical_coordinates)
>
> {
>
>
>
> /**
>
> * Create regular 3D cube of 4*4*4= 64 points.
>
> *
>
> */
>
> int nodes = 4;
>
> ControlPointContainerType controlPointContainer;
>
> My::InitRegularControlPointGrid<ControlPointContainerType,
> FixedImageType, InputPointType>(nodes, Dimension,
> controlPointContainer, imageF);
>
> ControlPointIteratorType it = controlPointContainer.begin();
>
>
>
> /*
>
> * We expect the control point with ID=63 to have maximum physical coordinates:
>
> *
>
> * Max index = [127, 127, 87]
>
> * Spacing = [3.125, 3.125, 4.0]
>
> * Max physical = [127, 127, 87] * [3.125, 3.125, 4.0] = [396.875,
> 396.875, 348]
>
> *
>
> */
>
> double pretty_tolerant = 3.0; // There will some rounding errors
> when computing the offset between control points.
>
> ControlPointType cp = controlPointContainer[63];
>
> BOOST_CHECK_CLOSE(cp.getLocation()[0], 396.875, pretty_tolerant);
>
> BOOST_CHECK_CLOSE(cp.getLocation()[1], 396.875, pretty_tolerant);
>
> BOOST_CHECK_CLOSE(cp.getLocation()[2], 348.0, pretty_tolerant);
>
> }
>
>
>
>
>
> BOOST_AUTO_TEST_CASE(test_resampler)
>
> {
>
> LinearInterpolatorType::Pointer linearInterpolator =
> LinearInterpolatorType::New();
>
> TransformType::Pointer MyTransform = TransformType::New();
>
>
>
> // Setup the cube
>
> int nodes = 4;
>
> ControlPointContainerType controlPointContainer;
>
> My::InitRegularControlPointGrid<ControlPointContainerType,
> FixedImageType, InputPointType>(nodes, Dimension,
> controlPointContainer, imageF);
>
> MyTransform->SetControlPoints(controlPointContainer);
>
> MyTransform->SetActiveControlPoints(controlPointContainer.size());
>
> MyTransform->SetGridSpacing(imageF->GetSpacing());
>
> MyTransform->SetGridOrigin(imageF->GetOrigin());
>
> MyTransform->BuildKdTree();
>
> FixedImageType::SizeType fixedImageSize =
> imageF->GetLargestPossibleRegion().GetSize();
>
> MyTransform->ComputeKernelSize(nodes, fixedImageSize);
>
>
>
> // Add parameters
>
> unsigned int numberOfBSplineParameters = MyTransform->GetNumberOfParameters();
> std::cout << "Number of bspline param = " <<
> numberOfBSplineParameters << std::endl;
> // Init all parameters to 0.0.
> ParametersType parameters(numberOfBSplineParameters);
> parameters.Fill(0.0);
> MyTransform->SetParameters(parameters);
>
>
>
>
> ResampleFilterType::Pointer resampler = ResampleFilterType::New();
> resampler->SetInput(imageM);
> resampler->SetTransform(MyTransform);
> resampler->SetInterpolator(linearInterpolator);
> resampler->SetOutputOrigin(imageF->GetOrigin());
> resampler->SetOutputSpacing(imageF->GetSpacing());
> resampler->SetSize(imageF->GetLargestPossibleRegion().GetSize());
> resampler->SetDefaultPixelValue(0);
>
>
> try {
>
> resampler->Update();
>
> }
>
> catch( itk::ExceptionObject & err ) {
>
> std::cerr << "ExceptionObject caught !" << std::endl;
>
> std::cerr << err << std::endl;
>
> }
> CopyImageToImage<typename FixedImageType::Pointer,
> FixedImageType>(resampler->GetOutput(), imageR);
>
> display_image<FixedImageType, CImg_Image>(imageR);
>
>
>
>
>
>
>
>
>
> }
>
>
>
>
>
>
>
> BOOST_AUTO_TEST_SUITE_END();
>
>
>
> BOOST_AUTO_TEST_SUITE_END();
>
>
>
>
>
> Any ideas on what migh be causing this segmentation fault is most
> welcome (have almost spend a week now trying to find the cause).
>
>
>
>
>
> On Thu, Nov 26, 2009 at 7:58 PM, motes motes <mort.motes at gmail.com> wrote:
>> On Thu, Nov 26, 2009 at 7:42 PM, Luis Ibanez <luis.ibanez at kitware.com> wrote:
>>> Hi Motes,
>>>
>>>
>>> 1) The Scope of the "try" block doesn't destroy any of the
>>> components of the registration framework.
>>>
>>
>> Ok that was also my understanding.
>>
>>
>>>
>>> 2) I don't see how you are passing any arguments to the
>>> function:
>>>
>>> computeRegisteredImage();
>>>
>>> are you declaring them as global ?
>>>
>>>
>>
>> Yes the registration components and images are "global" fields in a
>> wrapper class that I have written.
>>
>>
>>> 3) No offense,
>>> but your custom Transform is the main suspect.
>>>
>>> Here is the important question:
>>>
>>>
>>> Have you written a Unit Test
>>> for your new Transform ?
>>>
>>>
>>> I can anticipate that the answer is "no", :-)
>>>
>>> and therefore
>>>
>>> I will strongly,
>>> very, very, very strongly suggest
>>>
>>> that you write such unit test and fully debug
>>> your transform before you attempt to use it
>>> as part of a larger program.
>>>
>>> For examples on how to test Transforms,
>>> Please look at the files
>>>
>>> Insight/Testing/Code/Common/
>>>
>>> itkAffineTransformTest.cxx
>>> itkAzimuthElevationToCartesianTransformTest.cxx
>>> itkBSplineDeformableTransformTest2.cxx
>>> itkBSplineDeformableTransformTest.cxx
>>> itkCenteredAffineTransformTest.cxx
>>> itkCenteredEuler3DTransformTest.cxx
>>> itkCenteredRigid2DTransformTest.cxx
>>> itkCenteredTransformInitializerTest.cxx
>>> itkCenteredVersorTransformInitializerTest.cxx
>>> itkEuler2DTransformTest.cxx
>>> itkEuler3DTransformTest.cxx
>>> itkFixedCenterOfRotationAffineTransformTest.cxx
>>> itkIdentityTransformTest.cxx
>>> itkImageTransformTest.cxx
>>> itkLandmarkBasedTransformInitializerTest.cxx
>>> itkQuaternionRigidTransformTest.cxx
>>> itkRigid2DTransformTest.cxx
>>> itkRigid3DPerspectiveTransformTest.cxx
>>> itkRigid3DTransformTest.cxx
>>> itkRigid3DTransformTest.cxx.orig
>>> itkScaleLogarithmicTransformTest.cxx
>>> itkScaleSkewVersor3DTransformTest.cxx
>>> itkScaleTransformTest.cxx
>>> itkSimilarity2DTransformTest.cxx
>>> itkSimilarity3DTransformTest.cxx
>>> itkSplineKernelTransformTest.cxx
>>> itkTransformsSetParametersTest.cxx
>>> itkTransformTest.cxx
>>> itkTranslationTransformTest.cxx
>>> itkVersorRigid3DTransformTest.cxx
>>> itkVersorTransformTest.cxx
>>>
>>>
>>>
>>
>> Actually all new code I write is done through unit-tests. But thanks
>> for the above suggestions!
>>
>>
>>
>>
>>
>>
>>> Regards,
>>>
>>>
>>> Luis
>>>
>>>
>>> ------------------------------------------------------------------------------------------
>>> On Wed, Nov 25, 2009 at 8:44 PM, motes motes <mort.motes at gmail.com> wrote:
>>>> Which objects survive in the transform component when the image
>>>> registration method is finished?
>>>>
>>>> By finished I mean reaches the scope after this try/catch block:
>>>>
>>>> try {
>>>> registration->Update();
>>>> }
>>>> catch( itk::ExceptionObject & err ) {
>>>> std::cerr << "ExceptionObject caught !" << std::endl;
>>>> std::cerr << err << std::endl;
>>>> return EXIT_FAILURE;
>>>> }
>>>>
>>>> // What lives in the image registration method when arriving here?
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> The reason I ask is because I am trying to compute the registered
>>>> image AFTER the registration method is done. Something like this:
>>>>
>>>> try {
>>>>
>>>> registration->Update();
>>>> }
>>>> catch( itk::ExceptionObject & err ) {
>>>> std::cerr << "ExceptionObject caught !" << std::endl;
>>>> std::cerr << err << std::endl;
>>>> return EXIT_FAILURE;
>>>> }
>>>>
>>>> computeRegisteredImage();
>>>> return EXIT_SUCCESS;
>>>>
>>>>
>>>> where:
>>>>
>>>> void computeRegisteredImage(){
>>>> ResampleFilterType::Pointer resampler = ResampleFilterType::New();
>>>> resampler->SetInput(imageM);
>>>> resampler->SetTransform(registration->GetTransform());
>>>> resampler->SetInterpolator(registration->GetInterpolator());
>>>> resampler->SetOutputOrigin(imageF->GetOrigin());
>>>> resampler->SetOutputSpacing(imageF->GetSpacing());
>>>> resampler->SetSize(imageF->GetLargestPossibleRegion().GetSize());
>>>> resampler->SetDefaultPixelValue(default_pixel_value);
>>>> try {
>>>> resampler->Update();
>>>> }
>>>> catch( itk::ExceptionObject & err ) {
>>>> std::cerr << "ExceptionObject caught !" << std::endl;
>>>> std::cerr << err << std::endl;
>>>> }
>>>>
>>>>
>>>> imageR = resampler->GetOutput();
>>>>
>>>>
>>>> }
>>>>
>>>>
>>>>
>>>> But:
>>>>
>>>> resampler->Update();
>>>>
>>>> never returns. It thows a segmentation fault. I have located the error
>>>> to be in the TransformPoint method. I have made my own transform which
>>>> uses a kdTree in the TransformPoint function. This kdTree is stored
>>>> locally in my transform like:
>>>>
>>>>
>>>> #include "kd_tree.h" // My own itk:kdTree wrapper.
>>>> namespace itk
>>>> {
>>>> template <
>>>> class TScalarType = double, // Data type for scalars
>>>> unsigned int NDimensions = 3, // Number of dimensions
>>>> unsigned int VSplineOrder = 3 > // Spline order
>>>> class ITK_EXPORT MyDeformableTransform :
>>>> public Transform< TScalarType, NDimensions, NDimensions >
>>>> {
>>>> public:
>>>> /** Standard class typedefs. */
>>>> typedef MyDeformableTransform Self;
>>>> ...
>>>> ...
>>>> ...
>>>>
>>>> typedef KdTree<VectorType, ControlPointContainerType, NDimensions>
>>>> KdTreeType;
>>>>
>>>>
>>>> ....
>>>> .....
>>>> ...
>>>> protected:
>>>> mutable KdTreeType m_kdTree;
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> Can this have something to do with the segmentation error? Does it
>>>> need to be stored as a smart pointer instead?
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> In the TransformPoint function I do:
>>>>
>>>> // Transform a point
>>>> template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
>>>> void
>>>> MyDeformableTransform<TScalarType, NDimensions, VSplineOrder>
>>>> ::TransformPoint(
>>>> const InputPointType & point,
>>>> OutputPointType & outputPoint,
>>>> WeightsType & weights,
>>>> ParameterIndexArrayType & indices,
>>>> bool & inside ) const
>>>> {
>>>> unsigned int j;
>>>> IndexType supportIndex;
>>>> inside = true;
>>>> InputPointType transformedPoint;
>>>>
>>>> if ( m_BulkTransform ) {
>>>> transformedPoint = m_BulkTransform->TransformPoint( point );
>>>> } else {
>>>> transformedPoint = point;
>>>> }
>>>>
>>>> outputPoint.Fill( NumericTraits<ScalarType>::Zero );
>>>> if ( GetNumberOfParameters()>0) {
>>>> outputPoint.Fill( NumericTraits<ScalarType>::Zero );
>>>>
>>>>
>>>>
>>>> NeighborContainerType neighbors;
>>>>
>>>> // This results in a segmentation error
>>>> this->FindNeighbors(point, m_kernel_radius, neighbors);
>>>>
>>>>
>>>>
>>>> its the call:
>>>>
>>>> this->FindNeighbors(point, m_kernel_radius, neighbors);
>>>>
>>>> that gives the error which basically just calls the itk::kdTree class.
>>>>
>>>> Any suggestions are most welcome!
>>>> _____________________________________
>>>> Powered by www.kitware.com
>>>>
>>>> Visit other Kitware open-source projects at
>>>> http://www.kitware.com/opensource/opensource.html
>>>>
>>>> Kitware offers ITK Training Courses, for more information visit:
>>>> http://www.kitware.com/products/protraining.html
>>>>
>>>> Please keep messages on-topic and check the ITK FAQ at:
>>>> http://www.itk.org/Wiki/ITK_FAQ
>>>>
>>>> Follow this link to subscribe/unsubscribe:
>>>> http://www.itk.org/mailman/listinfo/insight-users
>>>>
>>>
>>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.html
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>
More information about the Insight-users
mailing list