[Insight-users] variogram / cancel update method

Mon Aug 2 08:06:15 EDT 2004

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: itkImageToVariogram.h
Type: text/x-chdr
Size: 5675 bytes
Desc: not available
Url : http://public.kitware.com/pipermail/insight-users/attachments/20040802/94016506/itkImageToVariogram-0001.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;
}

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 third to 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;
  }
}

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

  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(currentPixelValue == this->background_intensity)
      {
        progress.CompletedPixel();
        continue;
      }

      // 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(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;
       }

      progress.CompletedPixel();
   }
}

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