[Insight-users] variogram / cancel update method

Sat Aug 28 11:06:07 EDT 2004

I would like to now submit my variogram code for inclusion.

I have included an example and a test case.

As I said before the code is quite CPU intensive, if anyone can see any
obvious optimisations I would be very interested.

Thanks

Glenn Pierce

On Tue, 2004-08-03 at 09:47 -0400, Luis Ibanez wrote:
> Hi Glenn,
> 
> The process for aborting is closely related to the Progress report
> and/or the Iteration report. When your filter is running (e.g. in
> the GenerateData() method. You should be invoking Progress events
> and/or Iteration events from time to time.
> 
> The GUI that listens to these events should offer the user a "Cancel"
> button where she/he can click in order to stop the execution of the
> filter. The Callback of such button just need to invoke in the ITK
> filter the method:
> 
>                filter->SetAbortGenerateData(true);
> 
> This will stop the execution of the filter.
> 
> Note that before attempting to run the pipeline again. You must
> invoke the method ResetPipeline() since many filters will be in
> the middle of their Update() process.
> 
> 
> Question:
> 
>     The Variogram code that you sent, is it intended as a contribution
>     to ITK, or just for the purpose of demonstrating the problem that
>     your are facing ?
> 
> 
> 
> Please let us know,
> 
> 
>     Thanks
> 
> 
> 
>          Luis
> 
> 
> 
> 
> ---------------------------------
> Glenn Pierce wrote:
> 
> > Hi
> > 
> > I have developed a variogram filter that still needs a bit of tidying
> > work. I have attached it.
> > 
> > This variogram is very processor intensive. 
> > In the worst case O(N^2), where N is the number of pixels of an image.
> > (In reality there are optimisations I have placed in means the
> > performance is quite a bit better than this). However for large images
> > it would still be slow. 
> > 
> > 
> > For my main project I am calling the update method of the variogram
> > filter in a separate thread, and wanted to allow users to cancel the
> > execution.
> > 
> > For this reason I wondered if there is an accepted way in ITK or any
> > examples to cancel the update method or in my case tell the variogram
> > filter to stop the loop in the update method ?
> > 
> > 
> > Thanks for the help
> > 
> > Glenn
> > 
> > 
> >
> 
> 
> 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: itkVariogramTest.cxx
Type: text/x-c++src
Size: 2953 bytes
Desc: not available
Url : http://public.kitware.com/pipermail/insight-users/attachments/20040828/d84e7a22/itkVariogramTest.cxx
-------------- next part --------------
A non-text attachment was scrubbed...
Name: itkVariogramExample.cxx
Type: text/x-c++src
Size: 4442 bytes
Desc: not available
Url : http://public.kitware.com/pipermail/insight-users/attachments/20040828/d84e7a22/itkVariogramExample.cxx
-------------- next part --------------
A non-text attachment was scrubbed...
Name: itkImageToVariogram.h
Type: text/x-chdr
Size: 5581 bytes
Desc: not available
Url : http://public.kitware.com/pipermail/insight-users/attachments/20040828/d84e7a22/itkImageToVariogram.h
-------------- next part --------------
/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkImageToVariogram.txx,v $
  Language:  C++
  Date:      $Date: 2004/05/16 19:03:40 $
  Version:   $Revision: 1.1 $

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#ifndef _itkImageToVariogram_txx
#define _itkImageToVariogram_txx

#include "itkProgressReporter.h"
#include "itkConstantBoundaryCondition.h"

namespace itk{ 
namespace Statistics{

template < class TImage >
ImageToVariogram< TImage >
::ImageToVariogram()
{
  this->semi_variance_maximum = 0;
  this->background_intensity = 0;
  this->has_background = false;
}

template < class TImage >
ImageToVariogram< TImage >
::~ImageToVariogram()
{
  delete lagContainer;
}

template < class TImage >
void
ImageToVariogram< TImage >
::PrintSelf(std::ostream& os, Indent indent) const
{
  Superclass::PrintSelf(os,indent);
}

template < class TImage >
void
ImageToVariogram< TImage >
::SetImage(const TImage* image)
{
  m_Image = image ;

  // get spacing between data points
  spacing_x = m_Image->GetSpacing()[0];
  spacing_y = m_Image->GetSpacing()[1];

  // Set the max lag distance to use as one quater of the image diagonal
  this->max_lag = MaximumPossibleLag() / 4;

  IndexType index = m_Image->ComputeIndex( this->Size() - 1 );
  width = index[0];
  height = index[1];
}

/** returns the number of pixels in this container */
template < class TImage >
inline unsigned int
ImageToVariogram< TImage >
::Size() const
{
  return m_Image->GetPixelContainer()->Size() ;
}

/** returns the number of lags in this image */
template < class TImage >
unsigned int
ImageToVariogram< TImage >
::TotalLags()
{
  long unsigned int size = this->Size();
  // The sum of a arithmetrix series given by (number of terms) * (average of first and last terms)
  return ((size * (size - 1)) / 2);
}

/** returns the maximum lag in this image */
template < class TImage >
inline int
ImageToVariogram< TImage >
::MinimumLag()
{
  // Get the distance between the first and the next
  return this->CalculateLag(m_Image->ComputeIndex( 0 ), m_Image->ComputeIndex( 1 ));
}

/** returns the maximum lag in this image */
template < class TImage >
inline int
ImageToVariogram< TImage >
::MaximumPossibleLag()
{
  // Get the distance between the first and last pixel
  return static_cast<int>(this->CalculateLag(m_Image->ComputeIndex( 0 ),
    m_Image->ComputeIndex( this->Size() - 1 )));
}

/** Set the max lag that the variogram shall be calculated to
    The default is 1/4 of the maximum lag of the image diagonal */
template < class TImage >
inline void
ImageToVariogram< TImage >
::SetMaxLag(int lag)
{
  // Get the maximum lag distance to operate over
  if(lag > 0 && lag <= this->MaximumPossibleLag())
  {
    this->max_lag = lag;
  }
}

/** Set the background intensity that will be ignored */
template < class TImage >
inline void
ImageToVariogram< TImage >
::SetBackgroundIntensity(PixelType value)
{
  this->has_background = true;
  this->background_intensity = value;
}

template < class TImage >
void
ImageToVariogram< TImage >
::GenerateData()
{
  unsigned int size = this->Size();
  int lag;

  lagContainer = new LagDetail[this->max_lag + 1];

  // Create a progress reporter
  ProgressReporter progress(this, 0, size);

  ConstIteratorType mainPixelIter( m_Image, m_Image->GetRequestedRegion() );
  ConstIteratorType nextPixelIter( m_Image, m_Image->GetRequestedRegion() );

  IndexType nextPixelIndex;
  IndexType index;
  int max_x, min_x, max_y;

  for ( mainPixelIter.GoToBegin(); !mainPixelIter.IsAtEnd(); ++mainPixelIter)
   {
      IndexType currentPixelIndex = mainPixelIter.GetIndex();
      PixelType currentPixelValue = mainPixelIter.Get();

      if( (this->has_background == true) && (currentPixelValue == this->background_intensity) )
      {
        try
        {
          progress.CompletedPixel();
          continue;
        }
        catch( ProcessAborted & )
        {
          throw ProcessAborted(__FILE__,__LINE__);
          return;
        }
      }

      // Calculate the maximum and minumum pixels within the max lag circle
      // So we dont have to deal with all points
      max_x = currentPixelIndex[0] + this->max_lag;
      if(max_x > this->width)
      {
        max_x = this->width;
      }

      min_x = currentPixelIndex[0] - this->max_lag;
      if(min_x < 0)
      {
       min_x = 0;
      }

      // If the next pixel is further right than the max lag chosen move on to the next line
      max_y = currentPixelIndex[1] + this->max_lag;
      if(max_y > this->height)
      {
       max_y = this->height;
      }

      nextPixelIter = mainPixelIter;
      ++nextPixelIter;

      while (!nextPixelIter.IsAtEnd())
       {  
          nextPixelIndex = nextPixelIter.GetIndex();

          if( (this->has_background == true) && (nextPixelIter.Get() == this->background_intensity))
          {
            ++nextPixelIter;
            continue;
          }

          // If the next pixel is lower than the max lag chosen move on the current pixel
          if( nextPixelIndex[1] > max_y )
          {
            break;
          }

          // If the next pixel is further right than the max lag chosen move on to the next line
          if( ( nextPixelIndex[0] > max_x ) )
          {
            // move to next line
            index[0] = 0;
            index[1] = nextPixelIndex[1] + 1;

            nextPixelIter.SetIndex(index);
            nextPixelIndex = nextPixelIter.GetIndex();

            continue;
          }

          // If the next pixel is further left than the max lag chosen move on the 
          // first pixel within range.
          if( nextPixelIndex[0] < min_x )
          {
            // move to first pixel within reach
            index[0] = min_x;
            index[1] = nextPixelIndex[1];

            nextPixelIter.SetIndex(index);
            nextPixelIndex = nextPixelIter.GetIndex();

            continue;
          }

          lag = this->CalculateLag(currentPixelIndex, nextPixelIndex);

          // Check the pixel is within the max lag
          if (lag <= this->max_lag)
          {
              lagContainer[lag].count++;
              lagContainer[lag].sum_squared_diff += pow((currentPixelValue - nextPixelIter.Get()), 2);
          }

          ++nextPixelIter;
       }

       try
        {
          progress.CompletedPixel();
          continue;
        }
        catch( ProcessAborted & )
        {
          throw ProcessAborted(__FILE__,__LINE__);

          return;
        }
   }
}

template < class TImage >
inline double
ImageToVariogram< TImage >
::GetSemiVariance(int lag)
{
    if (lag > 0 && lag <= max_lag)
    {
      double tmp = lagContainer[lag].sum_squared_diff / (2.0 * lagContainer[lag].count);

      if(tmp > this->semi_variance_maximum)
       {
         this->semi_variance_maximum = tmp;
       }

      return tmp;
    }
    else
    {
      return -1.0;
    }
}

template < class TImage >
inline int
ImageToVariogram< TImage >
::CalculateLag(ImageToVariogram< TImage >::IndexType index1, ImageToVariogram< TImage >::IndexType index2)
{
  // Difference between pixels
  typename ImageToVariogram< TImage >::IndexType::IndexValueType x_pixel_difference = abs(index2[0] - index1[0]);
  typename ImageToVariogram< TImage >::IndexType::IndexValueType y_pixel_difference = abs(index2[1] - index1[1]);

  if(!y_pixel_difference)
    return x_pixel_difference;

  if(!x_pixel_difference)
    return y_pixel_difference;

  // Always round up
  int tmp = static_cast<int>(ceil(sqrt(pow(x_pixel_difference,2) + pow(y_pixel_difference,2))));

  return tmp;
}

} // end of namespace Statistics 
} // end of namespace itk

#endif