[Insight-users] StatisticsOnImage
Stefan Klein
stefan at isi.uu.nl
Thu Jun 22 05:23:43 EDT 2006
Hi,
The itkStatisticsOnImageFilter seems to have a small bug. When I run it on
an a 2D 256x256 image of doubles containing everywhere the same value
(0.2506507647...), it returns a negative variance and, consequently, an
undefined standard deviation:
var: -3.637423687e-014
std: -1.#IND
This is solved by applying the following change at line 297 of the .txx.
// unbiased estimate
variance = (sumOfSquares - (sum*sum / static_cast<RealType>(count)))
/ (static_cast<RealType>(count) - 1);
sigma = vcl_sqrt(variance);
change to:
// unbiased estimate
variance = (sumOfSquares - (sum*sum / static_cast<RealType>(count)))
/ (static_cast<RealType>(count) - 1);
// in case of numerical errors the variance might be <0.
variance = vnl_math_max(0.0, variance);
sigma = vcl_sqrt(variance);
Furthermore, I have a feature request for this class: Mask support. I've
added to this email the code that implements this (and fixes the bug).
Regards,
Stefan
-------------- next part --------------
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: itkStatisticsImageFilter.h,v $
Language: C++
Date: $Date: 2006/06/19 12:51:56 $
Version: $Revision: 1.1 $
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 __itkStatisticsImageFilter_h
#define __itkStatisticsImageFilter_h
#include "itkImageToImageFilter.h"
#include "itkNumericTraits.h"
#include "itkArray.h"
#include "itkSimpleDataObjectDecorator.h"
#include "itkSpatialObject.h"
namespace itk {
/** \class StatisticsImageFilter
* \brief Compute min. max, variance and mean of an Image.
*
* StatisticsImageFilter computes the minimum, maximum, sum, mean, variance
* sigma of an image. The filter needs all of its input image. It
* behaves as a filter with an input and output. Thus it can be inserted
* in a pipline with other filters and the statistics will only be
* recomputed if a downstream filter changes.
*
* The filter passes its input through unmodified. The filter is
* threaded. It computes statistics in each thread then combines them in
* its AfterThreadedGenerate method.
*
* \ingroup MathematicalStatisticsImageFilters
*/
template<class TInputImage>
class ITK_EXPORT StatisticsImageFilter :
public ImageToImageFilter<TInputImage, TInputImage>
{
public:
/** Standard Self typedef */
typedef StatisticsImageFilter Self;
typedef ImageToImageFilter<TInputImage,TInputImage> Superclass;
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Runtime information support. */
itkTypeMacro(StatisticsImageFilter, ImageToImageFilter);
/** Image related typedefs. */
typedef typename TInputImage::Pointer InputImagePointer;
typedef typename TInputImage::RegionType RegionType ;
typedef typename TInputImage::SizeType SizeType ;
typedef typename TInputImage::IndexType IndexType ;
typedef typename TInputImage::PointType PointType ;
typedef typename TInputImage::PixelType PixelType ;
/** Image related typedefs. */
itkStaticConstMacro(ImageDimension, unsigned int,
TInputImage::ImageDimension ) ;
/** Type to use for computations. */
typedef typename NumericTraits<PixelType>::RealType RealType;
/** Smart Pointer type to a DataObject. */
typedef typename DataObject::Pointer DataObjectPointer;
/** Type of DataObjects used for scalar outputs */
typedef SimpleDataObjectDecorator<RealType> RealObjectType;
typedef SimpleDataObjectDecorator<PixelType> PixelObjectType;
/** Mask support.
* Type for the mask of the input image. Only pixels that are "inside"
* this mask will be considered for the computation of the statistics */
typedef SpatialObject<
itkGetStaticConstMacro(ImageDimension)> MaskType;
typedef typename MaskType::Pointer MaskPointer;
/** Return the computed Minimum. */
PixelType GetMinimum() const
{ return this->GetMinimumOutput()->Get(); }
PixelObjectType* GetMinimumOutput();
const PixelObjectType* GetMinimumOutput() const;
/** Return the computed Maximum. */
PixelType GetMaximum() const
{ return this->GetMaximumOutput()->Get(); }
PixelObjectType* GetMaximumOutput();
const PixelObjectType* GetMaximumOutput() const;
/** Return the computed Mean. */
RealType GetMean() const
{ return this->GetMeanOutput()->Get(); }
RealObjectType* GetMeanOutput();
const RealObjectType* GetMeanOutput() const;
/** Return the computed Standard Deviation. */
RealType GetSigma() const
{ return this->GetSigmaOutput()->Get(); }
RealObjectType* GetSigmaOutput();
const RealObjectType* GetSigmaOutput() const;
/** Return the computed Variance. */
RealType GetVariance() const
{ return this->GetVarianceOutput()->Get(); }
RealObjectType* GetVarianceOutput();
const RealObjectType* GetVarianceOutput() const;
/** Return the compute Sum. */
RealType GetSum() const
{ return this->GetSumOutput()->Get(); }
RealObjectType* GetSumOutput();
const RealObjectType* GetSumOutput() const;
/** Make a DataObject of the correct type to be used as the specified
* output. */
virtual DataObjectPointer MakeOutput(unsigned int idx);
/** Set/Get Mask */
itkSetObjectMacro(Mask, MaskType);
itkGetConstObjectMacro(Mask, MaskType);
protected:
StatisticsImageFilter();
~StatisticsImageFilter(){};
void PrintSelf(std::ostream& os, Indent indent) const;
/** Pass the input through unmodified. Do this by Grafting in the AllocateOutputs method. */
void AllocateOutputs();
/** Initialize some accumulators before the threads run. */
void BeforeThreadedGenerateData ();
/** Do final mean and variance computation from data accumulated in threads. */
void AfterThreadedGenerateData ();
/** Multi-thread version GenerateData. */
void ThreadedGenerateData (const RegionType&
outputRegionForThread,
int threadId) ;
// Override since the filter needs all the data for the algorithm
void GenerateInputRequestedRegion();
// Override since the filter produces all of its output
void EnlargeOutputRequestedRegion(DataObject *data);
MaskPointer m_Mask;
private:
StatisticsImageFilter(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
Array<RealType> m_ThreadSum;
Array<RealType> m_SumOfSquares;
Array<long> m_Count;
Array<PixelType> m_ThreadMin;
Array<PixelType> m_ThreadMax;
} ; // end of class
} // end namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkStatisticsImageFilter.txx"
#endif
#endif
-------------- next part --------------
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: itkStatisticsImageFilter.txx,v $
Language: C++
Date: $Date: 2006/06/19 12:51:56 $
Version: $Revision: 1.1 $
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 _itkStatisticsImageFilter_txx
#define _itkStatisticsImageFilter_txx
#include "itkStatisticsImageFilter.h"
#include "itkImageRegionIterator.h"
#include "itkImageRegionConstIterator.h"
#include "itkImageRegionConstIteratorWithIndex.h"
#include "itkNumericTraits.h"
#include "itkProgressReporter.h"
namespace itk {
template<class TInputImage>
StatisticsImageFilter<TInputImage>
::StatisticsImageFilter(): m_ThreadSum(1), m_SumOfSquares(1), m_Count(1), m_ThreadMin(1), m_ThreadMax(1)
{
// first output is a copy of the image, DataObject created by
// superclass
//
// allocate the data objects for the outputs which are
// just decorators around pixel types
for (int i=1; i < 3; ++i)
{
typename PixelObjectType::Pointer output
= static_cast<PixelObjectType*>(this->MakeOutput(i).GetPointer());
this->ProcessObject::SetNthOutput(i, output.GetPointer());
}
// allocate the data objects for the outputs which are
// just decorators around real types
for (int i=3; i < 7; ++i)
{
typename RealObjectType::Pointer output
= static_cast<RealObjectType*>(this->MakeOutput(i).GetPointer());
this->ProcessObject::SetNthOutput(i, output.GetPointer());
}
this->GetMinimumOutput()->Set( NumericTraits<PixelType>::max() );
this->GetMaximumOutput()->Set( NumericTraits<PixelType>::NonpositiveMin() );
this->GetMeanOutput()->Set( NumericTraits<RealType>::max() );
this->GetSigmaOutput()->Set( NumericTraits<RealType>::max() );
this->GetVarianceOutput()->Set( NumericTraits<RealType>::max() );
this->GetSumOutput()->Set( NumericTraits<RealType>::Zero );
this->m_Mask = 0;
}
template<class TInputImage>
DataObject::Pointer
StatisticsImageFilter<TInputImage>
::MakeOutput(unsigned int output)
{
switch (output)
{
case 0:
return static_cast<DataObject*>(TInputImage::New().GetPointer());
break;
case 1:
return static_cast<DataObject*>(PixelObjectType::New().GetPointer());
break;
case 2:
return static_cast<DataObject*>(PixelObjectType::New().GetPointer());
break;
case 3:
case 4:
case 5:
case 6:
return static_cast<DataObject*>(RealObjectType::New().GetPointer());
break;
default:
// might as well make an image
return static_cast<DataObject*>(TInputImage::New().GetPointer());
break;
}
}
template<class TInputImage>
typename StatisticsImageFilter<TInputImage>::PixelObjectType*
StatisticsImageFilter<TInputImage>
::GetMinimumOutput()
{
return static_cast<PixelObjectType*>(this->ProcessObject::GetOutput(1));
}
template<class TInputImage>
const typename StatisticsImageFilter<TInputImage>::PixelObjectType*
StatisticsImageFilter<TInputImage>
::GetMinimumOutput() const
{
return static_cast<const PixelObjectType*>(this->ProcessObject::GetOutput(1));
}
template<class TInputImage>
typename StatisticsImageFilter<TInputImage>::PixelObjectType*
StatisticsImageFilter<TInputImage>
::GetMaximumOutput()
{
return static_cast<PixelObjectType*>(this->ProcessObject::GetOutput(2));
}
template<class TInputImage>
const typename StatisticsImageFilter<TInputImage>::PixelObjectType*
StatisticsImageFilter<TInputImage>
::GetMaximumOutput() const
{
return static_cast<const PixelObjectType*>(this->ProcessObject::GetOutput(2));
}
template<class TInputImage>
typename StatisticsImageFilter<TInputImage>::RealObjectType*
StatisticsImageFilter<TInputImage>
::GetMeanOutput()
{
return static_cast<RealObjectType*>(this->ProcessObject::GetOutput(3));
}
template<class TInputImage>
const typename StatisticsImageFilter<TInputImage>::RealObjectType*
StatisticsImageFilter<TInputImage>
::GetMeanOutput() const
{
return static_cast<const RealObjectType*>(this->ProcessObject::GetOutput(3));
}
template<class TInputImage>
typename StatisticsImageFilter<TInputImage>::RealObjectType*
StatisticsImageFilter<TInputImage>
::GetSigmaOutput()
{
return static_cast<RealObjectType*>(this->ProcessObject::GetOutput(4));
}
template<class TInputImage>
const typename StatisticsImageFilter<TInputImage>::RealObjectType*
StatisticsImageFilter<TInputImage>
::GetSigmaOutput() const
{
return static_cast<const RealObjectType*>(this->ProcessObject::GetOutput(4));
}
template<class TInputImage>
typename StatisticsImageFilter<TInputImage>::RealObjectType*
StatisticsImageFilter<TInputImage>
::GetVarianceOutput()
{
return static_cast<RealObjectType*>(this->ProcessObject::GetOutput(5));
}
template<class TInputImage>
const typename StatisticsImageFilter<TInputImage>::RealObjectType*
StatisticsImageFilter<TInputImage>
::GetVarianceOutput() const
{
return static_cast<const RealObjectType*>(this->ProcessObject::GetOutput(5));
}
template<class TInputImage>
typename StatisticsImageFilter<TInputImage>::RealObjectType*
StatisticsImageFilter<TInputImage>
::GetSumOutput()
{
return static_cast<RealObjectType*>(this->ProcessObject::GetOutput(6));
}
template<class TInputImage>
const typename StatisticsImageFilter<TInputImage>::RealObjectType*
StatisticsImageFilter<TInputImage>
::GetSumOutput() const
{
return static_cast<const RealObjectType*>(this->ProcessObject::GetOutput(6));
}
template<class TInputImage>
void
StatisticsImageFilter<TInputImage>
::GenerateInputRequestedRegion()
{
Superclass::GenerateInputRequestedRegion();
if ( this->GetInput() )
{
InputImagePointer image =
const_cast< typename Superclass::InputImageType * >( this->GetInput() );
image->SetRequestedRegionToLargestPossibleRegion();
}
}
template<class TInputImage>
void
StatisticsImageFilter<TInputImage>
::EnlargeOutputRequestedRegion(DataObject *data)
{
Superclass::EnlargeOutputRequestedRegion(data);
data->SetRequestedRegionToLargestPossibleRegion();
}
template<class TInputImage>
void
StatisticsImageFilter<TInputImage>
::AllocateOutputs()
{
// Pass the input through as the output
InputImagePointer image =
const_cast< TInputImage * >( this->GetInput() );
this->GraftOutput( image );
// Nothing that needs to be allocated for the remaining outputs
}
template<class TInputImage>
void
StatisticsImageFilter<TInputImage>
::BeforeThreadedGenerateData()
{
int numberOfThreads = this->GetNumberOfThreads();
// Resize the thread temporaries
m_Count.SetSize(numberOfThreads);
m_SumOfSquares.SetSize(numberOfThreads);
m_ThreadSum.SetSize(numberOfThreads);
m_ThreadMin.SetSize(numberOfThreads);
m_ThreadMax.SetSize(numberOfThreads);
// Initialize the temporaries
m_Count.Fill(NumericTraits<long>::Zero);
m_ThreadSum.Fill(NumericTraits<RealType>::Zero);
m_SumOfSquares.Fill(NumericTraits<RealType>::Zero);
m_ThreadMin.Fill(NumericTraits<PixelType>::max());
m_ThreadMax.Fill(NumericTraits<PixelType>::NonpositiveMin());
}
template<class TInputImage>
void
StatisticsImageFilter<TInputImage>
::AfterThreadedGenerateData()
{
int i;
long count;
RealType sumOfSquares;
int numberOfThreads = this->GetNumberOfThreads();
PixelType minimum;
PixelType maximum;
RealType mean;
RealType sigma;
RealType variance;
RealType sum;
sum = sumOfSquares = NumericTraits<RealType>::Zero;
count = 0;
// Find the min/max over all threads and accumulate count, sum and
// sum of squares
minimum = NumericTraits<PixelType>::max();
maximum = NumericTraits<PixelType>::NonpositiveMin();
for( i = 0; i < numberOfThreads; i++)
{
count += m_Count[i];
sum += m_ThreadSum[i];
sumOfSquares += m_SumOfSquares[i];
if (m_ThreadMin[i] < minimum)
{
minimum = m_ThreadMin[i];
}
if (m_ThreadMax[i] > maximum)
{
maximum = m_ThreadMax[i];
}
}
// compute statistics
mean = sum / static_cast<RealType>(count);
// unbiased estimate
variance = (sumOfSquares - (sum*sum / static_cast<RealType>(count)))
/ (static_cast<RealType>(count) - 1);
// in case of numerical errors the variance might be <0.
variance = vnl_math_max(0.0, variance);
sigma = vcl_sqrt(variance);
// Set the outputs
this->GetMinimumOutput()->Set( minimum );
this->GetMaximumOutput()->Set( maximum );
this->GetMeanOutput()->Set( mean );
this->GetSigmaOutput()->Set( sigma );
this->GetVarianceOutput()->Set( variance );
this->GetSumOutput()->Set( sum );
}
template<class TInputImage>
void
StatisticsImageFilter<TInputImage>
::ThreadedGenerateData(const RegionType& outputRegionForThread,
int threadId)
{
RealType realValue;
PixelType value;
// support progress methods/callbacks
ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
if ( this->m_Mask.IsNull() )
{
ImageRegionConstIterator<TInputImage> it (this->GetInput(), outputRegionForThread);
// do the work
while (!it.IsAtEnd())
{
value = it.Get();
realValue = static_cast<RealType>( value );
if (value < m_ThreadMin[threadId])
{
m_ThreadMin[threadId] = value;
}
if (value > m_ThreadMax[threadId])
{
m_ThreadMax[threadId] = value;
}
m_ThreadSum[threadId] += realValue;
m_SumOfSquares[threadId] += (realValue * realValue);
m_Count[threadId]++;
++it;
progress.CompletedPixel();
} // end while
} //end if
else
{ //use a mask
ImageRegionConstIteratorWithIndex<TInputImage> it (this->GetInput(), outputRegionForThread);
// do the work
while (!it.IsAtEnd())
{
PointType point;
this->GetInput()->
TransformIndexToPhysicalPoint(it.GetIndex(), point);
if ( this->m_Mask->IsInside( point ) )
{
value = it.Get();
realValue = static_cast<RealType>( value );
if (value < m_ThreadMin[threadId])
{
m_ThreadMin[threadId] = value;
}
if (value > m_ThreadMax[threadId])
{
m_ThreadMax[threadId] = value;
}
m_ThreadSum[threadId] += realValue;
m_SumOfSquares[threadId] += (realValue * realValue);
m_Count[threadId]++;
}
++it;
progress.CompletedPixel();
} // end while
} // end else
}
template <class TImage>
void
StatisticsImageFilter<TImage>
::PrintSelf(std::ostream& os, Indent indent) const
{
Superclass::PrintSelf(os,indent);
os << indent << "Minimum: "
<< static_cast<typename NumericTraits<PixelType>::PrintType>(this->GetMinimum()) << std::endl;
os << indent << "Maximum: "
<< static_cast<typename NumericTraits<PixelType>::PrintType>(this->GetMaximum()) << std::endl;
os << indent << "Sum: " << this->GetSum() << std::endl;
os << indent << "Mean: " << this->GetMean() << std::endl;
os << indent << "Sigma: " << this->GetSigma() << std::endl;
os << indent << "Variance: " << this->GetVariance() << std::endl;
}
}// end namespace itk
#endif
-------------- next part --------------
More information about the Insight-users
mailing list