[Insight-users] derivative computation of a 3D image: runtime error when image size gets bigger

[Gomez Herrero, Alberto]

Hello all,

I am trying to compute the first and second order derivatives of a 3D image using some smoothing, so I am using itk::RecursiveGaussianImageFilter. When I use images of, say, 200x200x200 voxels it works fine, but when I try with, for example 255x255x255 or bigger, it crashes at run time. I follow a simmililar scheme as in the ITK user guide for computing second order derivatives, and apparently the execution crashes while doing one of the duplicator->Update(); lines (not necessarily the first occurrence of that)
As the image size gets bigger, it crashes earlier (i.e. with 255x255x255 it crashes while computing Ixz, and with 260x260x260 it crashes at Iyz which is the previous computation).

Anyone know what I may be doing wrong? It looks like memory allocation issues to me, but as I use smart pointers I don't see where it can be.

Thank you very much,

Alberto

CODE:
typedef itk::RecursiveGaussianImageFilter<Image3DType,DoubleImage3DType > ShortFilterType;
typedef itk::RecursiveGaussianImageFilter<DoubleImage3DType,DoubleImage3DType > FilterType;
typedef itk::ImageDuplicator< DoubleImage3DType > DuplicatorType;
typedef itk::ImageDuplicator< Image3DType > ShortDuplicatorType;
DuplicatorType::Pointer duplicator = DuplicatorType::New();

//************** 1st order derivatives ***************************
ShortFilterType::Pointer dx = ShortFilterType::New();
dx->SetDirection( 0 );//dx works along x
if (sigmaG != 0)
      dx->SetSigma( sigmaG );
// Specify the derivative order in each direction
dx->SetFirstOrder();
// build the processing pipe
dx->SetInput( itkImage );

//update and save the output
//*** Ix
dx->Update();
duplicator->SetInputImage( dx->GetOutput() );
duplicator->Update();
DoubleImage3DType::Pointer Ix = duplicator->GetOutput(); //first order derivative along z
//*** Iy
dx->SetDirection( 1 );//dx works along y
dx->Update();
duplicator->Update();
DoubleImage3DType::Pointer Iy = duplicator->GetOutput();
//*** Iz
dx->SetDirection( 2 );//dx works along z
dx->Update();
duplicator->Update();
DoubleImage3DType::Pointer Iz = duplicator->GetOutput();
//************** 2nd order derivatives **************
//All derivatives are computed indepentdently
ShortFilterType::Pointer ga = ShortFilterType::New();
FilterType::Pointer gb = FilterType::New();
FilterType::Pointer gc = FilterType::New();


//DuplicatorType::Pointer duplicator = DuplicatorType::New();
ga->SetDirection( 0 );//ga works along x
gb->SetDirection( 1 );//gb works along y
gc->SetDirection( 2 );//gc works along z
if (sigma != 0){
      ga->SetSigma( sigma );
      gb->SetSigma( sigma );
      gc->SetSigma( sigma );
}
// Specify the derivative order in each direction
ga->SetZeroOrder();//No derivation in x (only smoothing)
gb->SetZeroOrder();// no derivation in y (only smoothing)
gc->SetSecondOrder();// 2nd order derivation in z (plus smoothing)
// build the processing pipe
ga->SetInput( itkImage );
gb->SetInput( ga->GetOutput() );
gc->SetInput( gb->GetOutput() );
duplicator->SetInputImage( gc->GetOutput() );
//update and save the output
//*** Izz
gc->Update();
duplicator->Update();
DoubleImage3DType::Pointer Izz = duplicator->GetOutput(); //second order derivative along z
//*** Iyy
gc->SetDirection( 1 ); // gc now works along y
gb->SetDirection( 2 ); // gb now works along z
gc->Update();
duplicator->Update();
DoubleImage3DType::Pointer Iyy = duplicator->GetOutput();
//*** Ixx
gc->SetDirection( 0 ); // gc now works along X
ga->SetDirection( 1 ); // ga now works along Y
gc->Update();
duplicator->Update();
DoubleImage3DType::Pointer Ixx = duplicator->GetOutput();
// Now we compute cross derivates
//*** Iyz=Izy
ga->SetDirection( 0 );
gb->SetDirection( 1 );
gc->SetDirection( 2 );
ga->SetZeroOrder();
gb->SetFirstOrder();
gc->SetFirstOrder();
gc->Update();
duplicator->Update();
DoubleImage3DType::Pointer Iyz = duplicator->GetOutput();
//*** Ixz=Izx
ga->SetDirection( 1 );
gb->SetDirection( 0 );
//gc->SetDirection( 2 );
//ga->SetZeroOrder();
//gb->SetFirstOrder();
//gc->SetFirstOrder();
gc->Update();
duplicator->Update();
DoubleImage3DType::Pointer Ixz = duplicator->GetOutput();
//*** Ixy = Iyx
ga->SetDirection( 2 );
//gb->SetDirection( 0 );
gc->SetDirection( 1 );
//ga->SetZeroOrder();
//gb->SetFirstOrder();
//gc->SetFirstOrder();
gc->Update();
duplicator->Update();
DoubleImage3DType::Pointer Ixy = duplicator->GetOutput();

________________________________
The information contained in this message may be confidential and legally protected under applicable law. The message is intended solely for the addressee(s). If you are not the intended recipient, you are hereby notified that any use, forwarding, dissemination, or reproduction of this message is strictly prohibited and may be unlawful. If you are not the intended recipient, please contact the sender by return e-mail and destroy all copies of the original message.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090526/706cc250/attachment-0001.htm>