[Insight-users] HessianRecursiveGaussian execution speed

[Dan Mueller]

Hi insight-users,

I am using the itk::HessianRecursiveGaussianImageFilter and have noticed 
it executes quite slow for 'large' 3-D datasets. Here are some times I 
recorded:

itkTestHessian TestHessian_128x128x128.mhd 1 1.0 >> Time = 4.28sec
itkTestHessian TestHessian_256x256x256.mhd 1 1.0 >> Time = 6.2e+003sec = 
1hr 43min
itkTestHessian TestHessian_512x512x512.mhd 1 1.0 >> I gave up...

*My questions:*
1. Have other people experienced similar execution times using this 
filter, or am I doing something wrong? I understand the filter performs 
separable convolution multiple times in multiple dimensions (ie. it's 
not a simple computation), but these times seem impractical for datasets 
256x256x256 or bigger...
2. Are there any techniques (besides relying on Moore's Law) to speed up 
this filter? I am interested in the /whole/ image, so using an image 
function would not help. (NOTE: I compiled my test in Release mode with 
optimization for maximum speed ie. /O2).

Thanks for any help and/or ideas.

Cheers, Dan
d.mueller at qut.edu.au
*
**My system:*
Platform: Windows XP SP2
Computer: Intel Pentium 3.00 GHz, dual core, 1GB physical RAM
Compiler: Microsoft Visual Studio 2005 (v8.0.50727.762 SP1)
CMake: 2.4.2
ITK: 3.2.0 CVS on Tues 27 Mar 2007

*My test code:*
/*=========================================================================
  itkTestHessian.cxx
=========================================================================*/

#include <iomanip>
#include "itkImage.h"
#include "itkNumericTraits.h"
#include "itkCommand.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkHessianRecursiveGaussianImageFilter.h"
#include "itkTimeProbe.h"

// Declare general types
const unsigned int Dimension = 3;
typedef float PixelType;
typedef itk::Image< PixelType, Dimension > ImageType;
typedef itk::ImageFileReader< ImageType > ReaderType;
typedef itk::ImageFileWriter< ImageType > WriterType;
typedef itk::HessianRecursiveGaussianImageFilter< ImageType > 
HessianFilterType;

// Declare a command to display the progress in a nice block format
class ProgressCommand : public itk::Command
{
public:
  typedef  ProgressCommand   Self;
  typedef  itk::Command             Superclass;
  typedef  itk::SmartPointer<Self>  Pointer;
  itkNewMacro( Self );

protected:
  ProgressCommand() { m_LastProgress = -1; };
  float m_LastProgress;

public:

  void Execute(itk::Object *caller, const itk::EventObject & event)
  {
    Execute( (const itk::Object *)caller, event);
  }

  void Execute(const itk::Object * object, const itk::EventObject & event)
  {
      const itk::ProcessObject* process = dynamic_cast< const 
itk::ProcessObject* >( object );

      if( ! itk::ProgressEvent().CheckEvent( &event ) )
        return;

      int fprogress = (process->GetProgress() * 100.0);
      int progress = (int)(process->GetProgress() * 100.0);
      if ((int)m_LastProgress == progress)
          return;
      if ((int)m_LastProgress != (progress - 1))
      {
          std::cout << std::setfill('0') << std::setw(3) << (progress - 
1) << " ";
          if (fprogress > 0.0 && (progress - 1) % 10 == 0)
            std::cout << std::endl;
      }
      if (fprogress > 0.0 && fprogress <= 100.0)
          std::cout << std::setfill('0') << std::setw(3) << progress << " ";
      if (fprogress > 0.0 && progress % 10 == 0)
          std::cout << std::endl;
      m_LastProgress = fprogress;
  }  
};

int main(int argc, char* argv[])
{
    // Display header
    std::cout << "Hessian Test -------------------------" << std::endl;
   
    // Check arguments
    if (argc < 4)
    {
        std::cerr << "Usage: " << argv[0] <<
                                " InputFilename" <<
                                " NormalizeAcrossScale" <<
                                " Sigma" <<
                                " \n";
        return EXIT_FAILURE;
    }

    // Setup the algorithm parameters
    unsigned int argn = 1;
    char* InputFilename        = argv[argn++];
    bool NormalizeAcrossScale  = (bool)atoi( argv[argn++] );
    double Sigma               = atof( argv[argn++] );

    // Display parameters
    std::cout << "InputFilename:" << InputFilename << std::endl;
    std::cout << "NormalizeAcrossScale:" << NormalizeAcrossScale << 
std::endl;
    std::cout << "Sigma:" << Sigma << std::endl;
    std::cout << "------------------------------------" << std::endl;
   
    try
    {
        // Read input image
        std::cout << "Reading input..." << std::endl;
        ReaderType::Pointer reader = ReaderType::New();
        reader->SetFileName( InputFilename );
        reader->Update();

        // Setup Hessian
        std::cout << "Computing Hessian..." << std::endl;
        HessianFilterType::Pointer filterHessian = HessianFilterType::New();
        filterHessian->SetInput( reader->GetOutput() );
        filterHessian->SetSigma( Sigma );
        filterHessian->SetNormalizeAcrossScale( NormalizeAcrossScale );
        ProgressCommand::Pointer observerHessian = ProgressCommand::New();
        filterHessian->AddObserver( itk::ProgressEvent(), observerHessian );
       
        // Compute and time hessian
        itk::TimeProbe time;
        time.Start();
        filterHessian->Update( );
        time.Stop();
        std::cout << std::setprecision(3) << "Time: " << 
time.GetMeanTime() << std::endl;
    }

    catch (itk::ExceptionObject & err)
    {
        std::cout << "ExceptionObject caught !" << std::endl;
        std::cout << err << std::endl;
        return EXIT_FAILURE;
    }

    //Return
    return EXIT_SUCCESS;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20070510/9224e220/attachment-0001.htm