[Insight-users] Accumulate a dimension

[Emiliano Beronich]

Thanks James for your suggestion. Starting from 
itkGetAverageSliceImageFilter, I have made a new filter that does what I 
wanted. It was almost completely rewritten, that's why I called it with 
a new name. I called itkAccumulateImageFilter because that is the main 
idea: Each pixel is the cumulative sum of the pixels along the collapsed 
dimension. The filter reduce the largest size of the accumulated 
dimension to 1 (only on the accumulated). I have tested it and it's working.

I think it's not very difficult to modify the filter to accumulate more
than one dimension (i.e. reducing a 4D to a 2D Image). It isn't what I
need at the moment, so I'm fine with it.

I would be happy to contribute the code to the community of ITK. I think 
this code could close the bug 1224, opened by myself.

Regards,
Emiliano





Miller, James V (Research) wrote:
> You might want to try the BasicFilters/itkGetAverageSliceImageFilter.h
> 
> This filter isn't well documented, but I believe its intention is to take
> a ND image and output a (N-1)D image where each pixel is the average of the
> pixels along the collapsed dimension.
> 
> I don't think this filter has been tested much.  So if you have a problem
> please let us know and we'll try to fix it.
> 
> Jim
> 
> 


-------------- next part --------------
/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkAccumulateImageFilter.h,v $
  Language:  C++
  Date:      $Date: 2004/10/06 12:00:00 $
  Version:   $Revision: 1.4 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     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 __itkAccumulateImageFilter_h
#define __itkAccumulateImageFilter_h

#include "itkImageToImageFilter.h"

namespace itk
{
  
/** \class AccumulateImageFilter
 * \brief Implements an accumulation of an image along a selected direction.
 *
 *    This class accumulates an image along a dimension and reduce the size 
 * of this dimension to 1. The dimension being accumulated is set by 
 * AccumulatedDimension. 
 *
 *   Each pixel is the cumulative sum of the pixels along the collapsed
 * dimension and reduce the size of the accumulated dimension to 1 (only 
 * on the accumulated). 
 *
 *   The dimensions of the InputImage and the OutputImage must be the same.
 *
 *
 *
 * This class is parameterized over the type of the input image and
 * the type of the output image.
 *
 *
 * \author Emiliano Beronich

 * This filter was contributed by Emiliano Beronich
 *
 * \ingroup   IntensityImageFilters     Singlethreaded
 */
template <class TInputImage, class TOutputImage>
class ITK_EXPORT AccumulateImageFilter : public ImageToImageFilter<TInputImage,TOutputImage>
{
public:
  /** Standard class typedefs. */
  typedef AccumulateImageFilter  Self;
  typedef ImageToImageFilter<TInputImage,TOutputImage>  Superclass;
  typedef SmartPointer<Self>   Pointer;
  typedef SmartPointer<const Self>  ConstPointer;

  /** Method for creation through the object factory. */
  itkNewMacro(Self);

  /** Run-time type information (and related methods). */
  itkTypeMacro(AccumulateImageFilter, ImageToImageFilter);

  /** Some convenient typedefs. */
  typedef TInputImage InputImageType;
  typedef typename    InputImageType::Pointer    InputImagePointer;
  typedef typename    InputImageType::RegionType InputImageRegionType;
  typedef typename    InputImageType::PixelType  InputImagePixelType;
  typedef TOutputImage OutputImageType;
  typedef typename     OutputImageType::Pointer    OutputImagePointer;
  typedef typename     OutputImageType::RegionType OutputImageRegionType;
  typedef typename     OutputImageType::PixelType  OutputImagePixelType;


  /** ImageDimension enumeration */
  itkStaticConstMacro(InputImageDimension, unsigned int,
                      TInputImage::ImageDimension);
  itkStaticConstMacro(OutputImageDimension, unsigned int,
                      TOutputImage::ImageDimension);

  /** Set the direction in which to accumulate the data.
  It must be set before the update of the filter.*/
  itkGetMacro( AccumulatedDimension, unsigned int );
  itkSetMacro( AccumulatedDimension, unsigned int );

  /** Perform a division by the size of the accumulated dimension after the
  accumulation is done. If true, the output image is the average of the accumulated
  dimension, if false the output is the sum of the pixels along the selected direction.
  The default value is false.**/
  itkSetMacro( AverageOutput, bool );
  itkGetMacro( AverageOutput, bool );

  

protected:
  AccumulateImageFilter();
  virtual ~AccumulateImageFilter() {};
  void PrintSelf(std::ostream& os, Indent indent) const;

  /** Apply changes to the output image information. */
  virtual void GenerateOutputInformation();

  /** Apply changes to the input image requested region. */
  virtual void GenerateInputRequestedRegion();

  /** This method implements the actual accumulation of the image.
   *
   * \sa ImageToImageFilter::ThreadedGenerateData(),
   *     ImageToImageFilter::GenerateData()  */
  void GenerateData(void);

private:
  AccumulateImageFilter(const Self&); //purposely not implemented
  void operator=(const Self&); //purposely not implemented

  unsigned int m_AccumulatedDimension;
  bool m_AverageOutput;

};

} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkAccumulateImageFilter.txx"
#endif

#endif


-------------- next part --------------
/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkAccumulateImageFilter.txx,v $
  Language:  C++
  Date:      $Date: 2004/10/06 12:00:00 $
  Version:   $Revision: 1.6 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     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 _itkAccumulateImageFilter_txx
#define _itkAccumulateImageFilter_txx

#include "itkAccumulateImageFilter.h"
#include "itkImageRegionIterator.h"
#include "itkImageRegionConstIterator.h"


namespace itk
{

/**
 * Constructor
 */
template <class TInputImage, class TOutputImage >
AccumulateImageFilter<TInputImage,TOutputImage >
::AccumulateImageFilter()
{
  this->SetNumberOfRequiredInputs( 1 );
  m_AccumulatedDimension=-1;
  m_AverageOutput = false;
}


template <class TInputImage, class TOutputImage>
void
AccumulateImageFilter<TInputImage,TOutputImage>
::GenerateOutputInformation()
{
  itkDebugMacro("GenerateOutputInformation Start");

  typename TOutputImage::RegionType outputRegion;
  typename TInputImage::IndexType inputIndex;
  typename TInputImage::SizeType  inputSize;
  typename TOutputImage::SizeType  outputSize;
  typename TOutputImage::IndexType outputIndex;
  typename TInputImage::SpacingType inSpacing;
  typename TInputImage::PointType inOrigin;
  typename TOutputImage::SpacingType outSpacing;
  typename TOutputImage::PointType outOrigin;

  // Get pointers to the input and output
  typename Superclass::OutputImagePointer output = this->GetOutput();
  typename Superclass::InputImagePointer input = const_cast< TInputImage * >( this->GetInput() );

  inputIndex = input->GetLargestPossibleRegion().GetIndex();
  inputSize = input->GetLargestPossibleRegion().GetSize();
  inSpacing = input->GetSpacing();
  inOrigin = input->GetOrigin();

// Set the LargestPossibleRegion of the output.
// Reduce the size of the accumulated dimension.
  for(int i = 0; i<InputImageDimension; i++) {
    if (i != m_AccumulatedDimension) {
      outputSize[i]  = inputSize[i];
      outputIndex[i] = inputIndex[i];
      outSpacing[i] = inSpacing[i];
      outOrigin[i]  = inOrigin[i];
    } else {
      outputSize[i]  = 1;
      outputIndex[i] = 0;
      outSpacing[i] = inSpacing[i]*inputSize[i];
      outOrigin[i]  = inOrigin[i] + (i-1)*inSpacing[i]/2;
    }
  }

  outputRegion.SetSize(outputSize);
  outputRegion.SetIndex(outputIndex);
  output->SetOrigin(outOrigin);
  output->SetSpacing(outSpacing);
  output->SetLargestPossibleRegion(outputRegion);

  itkDebugMacro("GenerateOutputInformation End");
}


template <class TInputImage, class  TOutputImage>
void
AccumulateImageFilter<TInputImage,TOutputImage>
::GenerateInputRequestedRegion()
{
  itkDebugMacro("GenerateInputRequestedRegion Start");
  Superclass::GenerateInputRequestedRegion();

  if ( this->GetInput() )
    {
      typename TInputImage::RegionType RequestedRegion;
      typename TInputImage::SizeType  inputSize;
      typename TInputImage::IndexType inputIndex;
      typename TInputImage::SizeType  inputLargSize;
      typename TInputImage::IndexType inputLargIndex;
      typename TOutputImage::SizeType  outputSize;
      typename TOutputImage::IndexType outputIndex;

      outputIndex = this->GetOutput()->GetRequestedRegion().GetIndex();
      outputSize = this->GetOutput()->GetRequestedRegion().GetSize();
      inputLargSize = this->GetInput()->GetLargestPossibleRegion().GetSize();
      inputLargIndex = this->GetInput()->GetLargestPossibleRegion().GetIndex();

      for(int i=0; i<TInputImage::ImageDimension; i++)
      {
        if(i!=m_AccumulatedDimension)
        {
          inputSize[i] = outputSize[i];
          inputIndex[i] = outputIndex[i];
        }
        else
        {
          inputSize[i]=inputLargSize[i];
          inputIndex[i]=inputLargIndex[i];
        }
      }

      RequestedRegion.SetSize(inputSize);
      RequestedRegion.SetIndex(inputIndex);
      InputImagePointer input = const_cast< TInputImage * > ( this->GetInput() );
      input->SetRequestedRegion (RequestedRegion);
    }

  itkDebugMacro("GenerateInputRequestedRegion End");
}


/**
 * GenerateData Performs the accumulation
 */
template <class TInputImage, class TOutputImage >
void
AccumulateImageFilter<TInputImage,TOutputImage>
::GenerateData( void )
{
  if(m_AccumulatedDimension<0 || m_AccumulatedDimension>=TInputImage::ImageDimension)
  {
    itkExceptionMacro(<<"AccumulateImageFilter: wrong Dimension to accumulate. Quitting!\n");
  }

  typename Superclass::InputImageConstPointer  inputImage = this->GetInput();
  typename TOutputImage::Pointer outputImage = this->GetOutput();
  outputImage->SetBufferedRegion( outputImage->GetRequestedRegion() );
  outputImage->Allocate();

// Accumulate over the Nth dimension ( = m_AccumulatedDimension)
// and divide by the size of the accumulated dimension.
  typedef itk::ImageRegionIterator<TOutputImage> outputIterType;
  outputIterType outputIter(outputImage, outputImage->GetBufferedRegion());
  typedef itk::ImageRegionConstIterator<TInputImage> inputIterType;
  TInputImage::RegionType AccumulatedRegion;
  TInputImage::SizeType AccumulatedSize = inputImage->GetLargestPossibleRegion().GetSize();
  TInputImage::IndexType AccumulatedIndex = inputImage->GetLargestPossibleRegion().GetIndex();
  TOutputImage::PixelType SizeAccumulatedDimension = AccumulatedSize[m_AccumulatedDimension];
  TOutputImage::PixelType IndexAccumulatedDimension = AccumulatedIndex[m_AccumulatedDimension];
  for(int i=0; i< InputImageDimension; i++) {
    if (i != m_AccumulatedDimension ) {
      AccumulatedSize[i] = 1;
    }
  }
  AccumulatedRegion.SetSize(AccumulatedSize);
  outputIter.GoToBegin();
  while(!outputIter.IsAtEnd())
  {
    TOutputImage::IndexType OutputIndex = outputIter.GetIndex();
    for(int i=0; i<InputImageDimension; i++) {
      if (i != m_AccumulatedDimension) {
        AccumulatedIndex[i] = OutputIndex[i];
      } else {
        AccumulatedIndex[i] = IndexAccumulatedDimension;
      }
    }
    AccumulatedRegion.SetIndex(AccumulatedIndex);
    inputIterType inputIter(inputImage, AccumulatedRegion);
    inputIter.GoToBegin();
    TOutputImage::PixelType Value=0;
    while(!inputIter.IsAtEnd())
    {
      Value+=inputIter.Get();
      ++inputIter;
    }
    outputIter.Set (Value);
    ++outputIter;
  }

// Performs a division to get an averaged output
  if (m_AverageOutput) {
    outputIter.GoToBegin();
    while(!outputIter.IsAtEnd())
    {
      outputIter.Set(outputIter.Get()/SizeAccumulatedDimension);
      ++outputIter;
    }
  }


}


template <class TInputImage, class TOutputImage >
void
AccumulateImageFilter<TInputImage,TOutputImage>::
PrintSelf(std::ostream& os, Indent indent) const
{
  Superclass::PrintSelf(os,indent);

  os << indent << "AccumulatedDimension: " << m_AccumulatedDimension << std::endl;
}


} // end namespace itk


#endif