[Insight-users] Cannot run the Resampling filter
motes motes
mort.motes at gmail.com
Sat Nov 28 05:14:21 EST 2009
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
>>>
>>
>
More information about the Insight-users
mailing list