[Insight-users] HessianRecursiveGaussian execution speed

Iván Macía imacia at vicomtech.es
Thu May 10 10:28:40 EDT 2007


Hi Dan,
 
I think the problem could be the amount of memory the filter is consuming in
the process. For an output with PixelType float the output data in the first
case is consuming :
 
128x128x128 tensors x 6 pixels/tensor x 2 bytes / pixel = 25.165.824 bytes
 
next case is 8 times the number of pixels so you have
 
256x256x256 tensors x 6 pixels/tensor x 2 bytes / pixel = 201.326.592 bytes
512x512x512 tensors x 6 pixels/tensor x 2 bytes / pixel = 1.610.612.736
bytes
 
but probably the filter is consuming more memory in the process itself. So
maybe in the second case you have already run out of physical memory. If you
are using double precision then multiply these numbers by 2 so in the second
case only with the output you are consuming 400 Mb. I would check out this
first before considering a problem in the filter itself.
 
Answering your second questions, for computing gaussian derivatives probably
recursive filters are among the fastest. Implementation in ITK follows
Deriche's work :
HYPERLINK
"http://citeseer.ist.psu.edu/deriche93recursively.html"http://citeseer.ist.p
su.edu/deriche93recursively.html
 
You may find useful information and references here
HYPERLINK
"http://en.wikipedia.org/wiki/Scale_space_implementation"http://en.wikipedia
.org/wiki/Scale_space_implementation
 
Regards
 
Ivan Macía
 
   _____  

De: insight-users-bounces+imacia=vicomtech.es at itk.org
[mailto:insight-users-bounces+imacia=vicomtech.es at itk.org] En nombre de Dan
Mueller
Enviado el: jueves, 10 de mayo de 2007 3:53
Para: insight-users
Asunto: [Insight-users] HessianRecursiveGaussian execution speed


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
HYPERLINK "mailto:d.mueller at qut.edu.au"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;
}



No virus found in this incoming message.
Checked by AVG Free Edition.
Version: 7.5.467 / Virus Database: 269.6.6/794 - Release Date: 08/05/2007
14:23



No virus found in this outgoing message.
Checked by AVG Free Edition. 
Version: 7.5.467 / Virus Database: 269.6.6/795 - Release Date: 09/05/2007
15:07
 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20070510/a08ab170/attachment.html


More information about the Insight-users mailing list