Sorry, forgot the files.
Here they are!
On 06/20/2013 09:20 AM, Carsten Kendziorra wrote:
Hello,
we've produced a shrink filter, that calculates the output pixel value
as the average of all corresponding input pixels.
This is basically a copy of the itkShrinkImageFilter (ITK 4.2.0), with
modifications in the ThreadedGenerateData() method, where the output
pixel values are calculated.
Author: Carsten Kendziorra, Henning Meyer, Marc Dewey. Charite Berlin,
Germany.
Maybe you want to include it in ITK.
Cheers,
Carsten
_______________________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
Kitware offers ITK Training Courses, for more information visit:
http://kitware.com/products/protraining.php
Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ
Follow this link to subscribe/unsubscribe:
http://www.itk.org/mailman/listinfo/insight-developers
#ifndef __itkShrinkAverageFilter_h
#define __itkShrinkAverageFilter_h
#include "itkImageToImageFilter.h"
namespace itk
{
/** \class ShrinkAverageFilter
* \brief Reduce the size of an image by an integer factor in each
* dimension and calculate the average value.
*
* ShrinkAverageFilter reduces the size of an image by an integer factor
* in each dimension. In contrast to itkShrinkImageFilter, this filter
* calculates the average pixel value from each corresponding input
* pixel. The algorithm implemented is a simple subsample.
*
* This is basicaly a copy of itkShrinkImageFilter, with modyfications
* in ThreadedGenerateData() method, where the output pixel values are
* calculated.
*
* \author Carsten Kendziorra, Henning Meyer, Marc Dewey. Charite Berlin, Germany.
*/
template< class TInputImage, class TOutputImage >
class ITK_EXPORT ShrinkAverageFilter:
public ImageToImageFilter< TInputImage, TOutputImage >
{
public:
/** Standard class typedefs. */
typedef ShrinkAverageFilter 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(ShrinkAverageFilter, ImageToImageFilter);
/** Typedef to images */
typedef TOutputImage OutputImageType;
typedef TInputImage InputImageType;
typedef typename OutputImageType::Pointer OutputImagePointer;
typedef typename InputImageType::Pointer InputImagePointer;
typedef typename InputImageType::ConstPointer InputImageConstPointer;
typedef typename TOutputImage::IndexType OutputIndexType;
typedef typename TInputImage::IndexType InputIndexType;
typedef typename TOutputImage::OffsetType OutputOffsetType;
typedef typename TInputImage::PixelType InputImagePixelType;
typedef typename TOutputImage::PixelType OutputImagePixelType;
/** Typedef to describe the output image region type. */
typedef typename TInputImage::RegionType InputImageRegionType;
/** Typedef to describe the output image region type. */
typedef typename TOutputImage::RegionType OutputImageRegionType;
/** ImageDimension enumeration. */
itkStaticConstMacro(ImageDimension, unsigned int,
TInputImage::ImageDimension);
itkStaticConstMacro(OutputImageDimension, unsigned int,
TOutputImage::ImageDimension);
typedef FixedArray< unsigned int, ImageDimension > ShrinkFactorsType;
/** Set the shrink factors. Values are clamped to
* a minimum value of 1. Default is 1 for all dimensions. */
itkSetMacro(ShrinkFactors, ShrinkFactorsType);
void SetShrinkFactors(unsigned int factor);
void SetShrinkFactor(unsigned int i, unsigned int factor)
{
m_ShrinkFactors[i] = factor;
}
/** Get the shrink factors. */
itkGetConstReferenceMacro(ShrinkFactors, ShrinkFactorsType);
/** ShrinkAverageFilter produces an image which is a different
* resolution and with a different pixel spacing than its input
* image. As such, ShrinkAverageFilter needs to provide an
* implementation for GenerateOutputInformation() in order to inform
* the pipeline execution model. The original documentation of this
* method is below.
* \sa ProcessObject::GenerateOutputInformaton() */
virtual void GenerateOutputInformation();
/** ShrinkAverageFilter needs a larger input requested region than the output
* requested region. As such, ShrinkAverageFilter needs to provide an
* implementation for GenerateInputRequestedRegion() in order to inform the
* pipeline execution model.
* \sa ProcessObject::GenerateInputRequestedRegion() */
virtual void GenerateInputRequestedRegion();
#ifdef ITK_USE_CONCEPT_CHECKING
/** Begin concept checking */
itkConceptMacro( InputConvertibleToOutputCheck,
( Concept::Convertible< typename TInputImage::PixelType, typename TOutputImage::PixelType > ) );
itkConceptMacro( SameDimensionCheck,
( Concept::SameDimension< ImageDimension, OutputImageDimension > ) );
/** End concept checking */
#endif
protected:
ShrinkAverageFilter();
~ShrinkAverageFilter() {}
void PrintSelf(std::ostream & os, Indent indent) const;
/** ShrinkAverageFilter can be implemented as a multithreaded filter.
* Therefore, this implementation provides a ThreadedGenerateData() routine
* which is called for each processing thread. The output image data is
* allocated automatically by the superclass prior to calling
* ThreadedGenerateData(). ThreadedGenerateData can only write to the
* portion of the output image specified by the parameter
* "outputRegionForThread"
*
* \sa ImageToImageFilter::ThreadedGenerateData(),
* ImageToImageFilter::GenerateData() */
void ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread,
ThreadIdType threadId);
private:
ShrinkAverageFilter(const Self &); //purposely not implemented
void operator=(const Self &); //purposely not implemented
ShrinkFactorsType m_ShrinkFactors;
};
} // end namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkShrinkAverageFilter.txx"
#endif
#endif
#ifndef __itkShrinkAverageFilter_txx
#define __itkShrinkAverageFilter_txx
#include "itkShrinkAverageFilter.h"
#include <itkImageRegionIterator.h>
#include <itkImageRegionConstIterator.h>
#include <itkContinuousIndex.h>
#include <itkObjectFactory.h>
#include <itkProgressReporter.h>
#include <itkVariableLengthVector.h>
#include <boost/type_traits.hpp>
#include <boost/mpl/if.hpp>
namespace itk
{
/**
*
*/
template< class TInputImage, class TOutputImage >
ShrinkAverageFilter< TInputImage, TOutputImage >
::ShrinkAverageFilter()
{
for ( unsigned int j = 0; j < ImageDimension; j++ )
{
m_ShrinkFactors[j] = 1;
}
}
/**
*
*/
template< class TInputImage, class TOutputImage >
void
ShrinkAverageFilter< TInputImage, TOutputImage >
::PrintSelf(std::ostream & os, Indent indent) const
{
Superclass::PrintSelf(os, indent);
os << indent << "Shrink Factor: ";
for ( unsigned int j = 0; j < ImageDimension; j++ )
{
os << m_ShrinkFactors[j] << " ";
}
os << std::endl;
}
/**
*
*/
template< class TInputImage, class TOutputImage >
void
ShrinkAverageFilter< TInputImage, TOutputImage >
::SetShrinkFactors(unsigned int factor)
{
unsigned int j;
for ( j = 0; j < ImageDimension; j++ )
{
if ( factor != m_ShrinkFactors[j] ) { break; }
}
if ( j < ImageDimension )
{
this->Modified();
for ( j = 0; j < ImageDimension; j++ )
{
m_ShrinkFactors[j] = factor;
if ( m_ShrinkFactors[j] < 1 )
{
m_ShrinkFactors[j] = 1;
}
}
}
}
/**
*
*/
template< class TInputImage, class TOutputImage >
void
ShrinkAverageFilter< TInputImage, TOutputImage >
::ThreadedGenerateData(const OutputImageRegionType &
outputRegionForThread,
ThreadIdType threadId)
{
// Get the input and output pointers
InputImageConstPointer inputPtr = this->GetInput();
OutputImagePointer outputPtr = this->GetOutput();
typename OutputImageType::SizeType outputSize =
outputRegionForThread.GetSize();
typename OutputImageType::IndexType outputIndex =
outputRegionForThread.GetIndex();
InputImageRegionType inputRegion =
inputPtr->GetLargestPossibleRegion();
typename InputImageType::SizeType inputSize;
typename InputImageType::IndexType inputIndex;
for ( unsigned int dim = 0; dim < TInputImage::ImageDimension;
dim++ )
{
inputSize[dim] = outputSize[dim] * m_ShrinkFactors[dim];
inputIndex[dim] = outputIndex[dim] *
m_ShrinkFactors[dim];
}
inputRegion.SetSize(inputSize);
inputRegion.SetIndex(inputIndex);
// number of all pixel/voxel of the output image
unsigned int outputPixelNumber = 1;
for ( unsigned int dim = 0; dim < TInputImage::ImageDimension;
dim++ )
{
outputPixelNumber *= outputSize[dim];
}
// set pixel type depending on input image
typedef typename boost::mpl::if_<
boost::is_integral<InputImagePixelType>, long int, double >::type ValType;
// sum of all pixel values of the corresponding input values
for each output pixel
typedef itk::VariableLengthVector<ValType> VectorValSumType;
VectorValSumType valSum;
valSum.SetSize(outputPixelNumber);
valSum.Fill(0);
// iterator for input image
typedef ImageRegionConstIterator< TInputImage >
InputImageIterator;
InputImageIterator inIt(inputPtr, inputRegion);
unsigned outIdx = 0;
typedef itk::VariableLengthVector<unsigned> VectorUnsignedType;
VectorUnsignedType outIdxA;
outIdxA.SetSize(OutputImageType::ImageDimension);
outIdxA.Fill(0);
InputIndexType indexInsideOutputBlock;
VectorUnsignedType outBlockSizePerDimension;
outBlockSizePerDimension.SetSize(OutputImageType::ImageDimension);
unsigned osize = 1;
for(unsigned dim = 0; dim < OutputImageType::ImageDimension;
dim++) {
outBlockSizePerDimension[dim] = osize;
osize *= outputSize[dim];
indexInsideOutputBlock[dim] = 0;
}
// walk the input image and sum all values for each
corresponding output pixel
while ( !inIt.IsAtEnd() )
{
valSum[outIdx] += inIt.Get();
++inIt;
unsigned dim = 0;
bool nextdim = false;
do {
indexInsideOutputBlock[dim]++;
nextdim = false;
if (indexInsideOutputBlock[dim] ==
m_ShrinkFactors[dim])
{
indexInsideOutputBlock[dim] = 0;
outIdx += outBlockSizePerDimension[dim];
outIdxA[dim]++;
if ( outIdxA[dim] == outputSize[dim] )
{
outIdxA[dim] = 0;
outIdx -=
outBlockSizePerDimension[dim+1];
nextdim = true;
}
}
dim++;
} while (nextdim);
}
// Define/declare an iterator that will walk the output region
for this
// thread.
typedef ImageRegionIterator< TOutputImage > OutputIterator;
OutputIterator outIt(outputPtr, outputRegionForThread);
// Number of input pixel per output pixel
int VoxelShrinkNumber = 1;
for ( unsigned int dim = 0; dim < TInputImage::ImageDimension;
dim++ )
{
VoxelShrinkNumber *= m_ShrinkFactors[dim];
}
outIdx = 0;
// walk the output pixel and set new value
while ( !outIt.IsAtEnd() )
{
outIt.Set( valSum[outIdx] / VoxelShrinkNumber );
++outIdx;
++outIt;
}
}
/**
*
*/
template< class TInputImage, class TOutputImage >
void
ShrinkAverageFilter< TInputImage, TOutputImage >
::GenerateInputRequestedRegion()
{
// Call the superclass' implementation of this method
Superclass::GenerateInputRequestedRegion();
// Get pointers to the input and output
InputImagePointer inputPtr = const_cast< TInputImage * >(
this->GetInput() );
OutputImagePointer outputPtr = this->GetOutput();
if ( !inputPtr || !outputPtr )
{
return;
}
// Compute the input requested region (size and start index)
// Use the image transformations to insure an input requested
region
// that will provide the proper range
unsigned int i;
const typename TOutputImage::SizeType &
outputRequestedRegionSize =
outputPtr->GetRequestedRegion().GetSize();
const typename TOutputImage::IndexType &
outputRequestedRegionStartIndex =
outputPtr->GetRequestedRegion().GetIndex();
// Convert the factor for convenient multiplication
typename TOutputImage::SizeType factorSize;
for ( i = 0; i < TInputImage::ImageDimension; i++ )
{
factorSize[i] = m_ShrinkFactors[i];
}
OutputIndexType outputIndex;
InputIndexType inputIndex, inputRequestedRegionIndex;
OutputOffsetType offsetIndex;
typename TInputImage::SizeType inputRequestedRegionSize;
typename TOutputImage::PointType tempPoint;
// Use this index to compute the offset everywhere in this class
outputIndex = outputPtr->GetLargestPossibleRegion().GetIndex();
// We wish to perform the following mapping of outputIndex to
// inputIndex on all points in our region
outputPtr->TransformIndexToPhysicalPoint(outputIndex,
tempPoint);
inputPtr->TransformPhysicalPointToIndex(tempPoint, inputIndex);
// Given that the size is scaled by a constant factor eq:
// inputIndex = outputIndex * factorSize
// is equivalent up to a fixed offset which we now compute
OffsetValueType zeroOffset = 0;
for ( i = 0; i < TInputImage::ImageDimension; i++ )
{
offsetIndex[i] = inputIndex[i] - outputIndex[i] *
m_ShrinkFactors[i];
// It is plausible that due to small amounts of loss of
numerical
// precision that the offset it negaive, this would
cause sampling
// out of out region, this is insurance against that
possibility
offsetIndex[i] = vnl_math_max(zeroOffset,
offsetIndex[i]);
}
inputRequestedRegionIndex = outputRequestedRegionStartIndex *
factorSize + offsetIndex;
// Originally this was
// for ( i=0; i < TInputImage::ImageDimension; ++i )
// {
// inputRequestedRegionSize[i] =
(outputRequestedRegionSize[i] - 1 ) *
// factorSize[i] + 1;
// }
// but with centered pixels we may sample edge to edge
inputRequestedRegionSize = outputRequestedRegionSize *
factorSize;
typename TInputImage::RegionType inputRequestedRegion;
inputRequestedRegion.SetIndex(inputRequestedRegionIndex);
inputRequestedRegion.SetSize(inputRequestedRegionSize);
inputRequestedRegion.Crop( inputPtr->GetLargestPossibleRegion()
);
inputPtr->SetRequestedRegion(inputRequestedRegion);
}
/**
*
*/
template< class TInputImage, class TOutputImage >
void
ShrinkAverageFilter< TInputImage, TOutputImage >
::GenerateOutputInformation()
{
// Call the superclass' implementation of this method
Superclass::GenerateOutputInformation();
// Get pointers to the input and output
InputImageConstPointer inputPtr = this->GetInput();
OutputImagePointer outputPtr = this->GetOutput();
if ( !inputPtr || !outputPtr )
{
return;
}
// Compute the output spacing, the output image size, and the
// output image start index
unsigned int i;
const typename TInputImage::SpacingType &
inputSpacing = inputPtr->GetSpacing();
const typename TInputImage::SizeType & inputSize =
inputPtr->GetLargestPossibleRegion().GetSize();
const typename TInputImage::IndexType & inputStartIndex =
inputPtr->GetLargestPossibleRegion().GetIndex();
typename TOutputImage::SpacingType outputSpacing;
typename TOutputImage::SizeType outputSize;
typename TOutputImage::IndexType outputStartIndex;
for ( i = 0; i < TOutputImage::ImageDimension; i++ )
{
outputSpacing[i] = inputSpacing[i] *
(double)m_ShrinkFactors[i];
// Round down so that all output pixels fit input input
region
outputSize[i] = static_cast<SizeValueType>(
vcl_floor( (double)inputSize[i] /
(double)m_ShrinkFactors[i] ) );
if ( outputSize[i] < 1 )
{
outputSize[i] = 1;
}
// Because of the later origin shift this starting
index is not
// critical
outputStartIndex[i] = static_cast<IndexValueType>(
vcl_ceil( (double)inputStartIndex[i] /
(double)m_ShrinkFactors[i] ) );
}
outputPtr->SetSpacing(outputSpacing);
// Compute origin offset
// The physical center's of the input and output should be the
same
ContinuousIndex< double, TOutputImage::ImageDimension >
inputCenterIndex;
ContinuousIndex< double, TOutputImage::ImageDimension >
outputCenterIndex;
for ( i = 0; i < TOutputImage::ImageDimension; i++ )
{
inputCenterIndex[i] = inputStartIndex[i] + (
inputSize[i] - 1 ) / 2.0;
outputCenterIndex[i] = outputStartIndex[i] + (
outputSize[i] - 1 ) / 2.0;
}
typename TOutputImage::PointType inputCenterPoint;
typename TOutputImage::PointType outputCenterPoint;
inputPtr->TransformContinuousIndexToPhysicalPoint(inputCenterIndex,
inputCenterPoint);
outputPtr->TransformContinuousIndexToPhysicalPoint(outputCenterIndex,
outputCenterPoint);
typename TOutputImage::PointType outputOrigin =
outputPtr->GetOrigin();
outputOrigin = outputOrigin + ( inputCenterPoint -
outputCenterPoint );
outputPtr->SetOrigin(outputOrigin);
// Set region
typename TOutputImage::RegionType outputLargestPossibleRegion;
outputLargestPossibleRegion.SetSize(outputSize);
outputLargestPossibleRegion.SetIndex(outputStartIndex);
outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
}
} // end namespace itk
#endif
_______________________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
Kitware offers ITK Training Courses, for more information visit:
http://kitware.com/products/protraining.php
Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ
Follow this link to subscribe/unsubscribe:
http://www.itk.org/mailman/listinfo/insight-developers
