[Insight-users] Cannot run the Resampling filter
motes motes
mort.motes at gmail.com
Sat Nov 28 08:47:28 EST 2009
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