diff --git a/include/otbStatisticsImageCustomFilter.h b/include/otbStatisticsImageCustomFilter.h deleted file mode 100644 index 8ff6fcad5493cd713c2cc8f557c710b5321c1921..0000000000000000000000000000000000000000 --- a/include/otbStatisticsImageCustomFilter.h +++ /dev/null @@ -1,197 +0,0 @@ -/*========================================================================= - - Copyright (c) Remi Cresson (IRSTEA). All rights reserved. - - - 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 __otbStatisticsImageCustomFilter_h -#define __otbStatisticsImageCustomFilter_h - -#include "itkImageToImageFilter.h" -#include "itkNumericTraits.h" -#include "itkArray.h" -#include "itkSimpleDataObjectDecorator.h" - -namespace otb -{ -/** \class StatisticsImageCustomFilter - * \brief Compute min. max, variance and mean of an Image. - * - * StatisticsImageCustomFilter 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. - * - * A value to ignore (in the stats) can be chosen - * - * \ingroup MathematicalStatisticsImageCustomFilters - * \ingroup SimpleExtractionTools - * - * \wiki - * \wikiexample{Statistics/StatisticsImageCustomFilter,Compute min\, max\, variance and mean of an Image.} - * \endwiki - */ -template< typename TInputImage > -class StatisticsImageCustomFilter: - public itk::ImageToImageFilter< TInputImage, TInputImage > -{ -public: - /** Standard Self typedef */ - typedef StatisticsImageCustomFilter Self; - typedef itk::ImageToImageFilter< TInputImage, TInputImage > Superclass; - typedef itk::SmartPointer< Self > Pointer; - typedef itk::SmartPointer< const Self > ConstPointer; - - /** Method for creation through the object factory. */ - itkNewMacro(Self); - - /** Runtime information support. */ - itkTypeMacro(StatisticsImageCustomFilter,itk::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::PixelType PixelType; - - /** Image related typedefs. */ - itkStaticConstMacro(ImageDimension, unsigned int, - TInputImage::ImageDimension); - - /** Type to use for computations. */ - typedef typename itk::NumericTraits< PixelType >::RealType RealType; - typedef unsigned long LongType; - - /** Smart Pointer type to a DataObject. */ - typedef typename itk::DataObject::Pointer DataObjectPointer; - - /** Type of DataObjects used for scalar outputs */ - typedef itk::SimpleDataObjectDecorator< RealType > RealObjectType; - typedef itk::SimpleDataObjectDecorator< LongType > LongObjectType; - typedef itk::SimpleDataObjectDecorator< PixelType > PixelObjectType; - - /** 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; - - /** Return the compute Sum of squares. */ - RealType GetSumOfSquares() const - { return this->GetSumOfSquaresOutput()->Get(); } - RealObjectType * GetSumOfSquaresOutput(); - - const RealObjectType * GetSumOfSquaresOutput() const; - - /** Return the compute count. */ - RealType GetCount() const - { return this->GetCountOutput()->Get(); } - LongObjectType * GetCountOutput(); - - const LongObjectType * GetCountOutput() const; - - /** Make a DataObject of the correct type to be used as the specified - * output. */ - typedef itk::ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType; - using Superclass::MakeOutput; - virtual DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx); - -#ifdef ITK_USE_CONCEPT_CHECKING - // Begin concept checking - itkConceptMacro( InputHasNumericTraitsCheck, - ( itk::Concept::HasNumericTraits< PixelType > ) ); - // End concept checking -#endif - - /** Accessors for ignored value */ - itkGetMacro(IgnoredInputValue, PixelType); - itkSetMacro(IgnoredInputValue, PixelType); - -protected: - StatisticsImageCustomFilter(); - ~StatisticsImageCustomFilter(){} - void PrintSelf(std::ostream & os, itk::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, - itk::ThreadIdType threadId); - -private: - StatisticsImageCustomFilter(const Self &); //purposely not implemented - void operator=(const Self &); //purposely not implemented - - itk::Array< RealType > m_ThreadSum; - itk::Array< RealType > m_SumOfSquares; - itk::Array< itk::SizeValueType > m_Count; - itk::Array< PixelType > m_ThreadMin; - itk::Array< PixelType > m_ThreadMax; - - PixelType m_IgnoredInputValue; -}; // end of class -} // end namespace itk - -#ifndef ITK_MANUAL_INSTANTIATION -#include "otbStatisticsImageCustomFilter.hxx" -#endif - -#endif diff --git a/include/otbStatisticsImageCustomFilter.hxx b/include/otbStatisticsImageCustomFilter.hxx deleted file mode 100644 index 8e801dbe3f6b0919db688ff8876d22e5762e2c3d..0000000000000000000000000000000000000000 --- a/include/otbStatisticsImageCustomFilter.hxx +++ /dev/null @@ -1,395 +0,0 @@ -/*========================================================================= - - Copyright (c) Remi Cresson (IRSTEA). All rights reserved. - - - 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 __otbStatisticsImageCustomFilter_hxx -#define __otbStatisticsImageCustomFilter_hxx -#include "otbStatisticsImageCustomFilter.h" - - -#include "itkImageScanlineIterator.h" -#include "itkProgressReporter.h" - -namespace otb -{ -template< typename TInputImage > -StatisticsImageCustomFilter< TInputImage > -::StatisticsImageCustomFilter():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->itk::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->itk::ProcessObject::SetNthOutput( i, output.GetPointer() ); - } - - // allocate the data objects for the count output (long type) - typename LongObjectType::Pointer output = - static_cast< LongObjectType * >( this->MakeOutput(7).GetPointer() ); - this->itk::ProcessObject::SetNthOutput( 7, output.GetPointer() ); - - this->GetMinimumOutput()->Set( itk::NumericTraits< PixelType >::max() ); - this->GetMaximumOutput()->Set( itk::NumericTraits< PixelType >::NonpositiveMin() ); - this->GetMeanOutput()->Set( itk::NumericTraits< RealType >::max() ); - this->GetSigmaOutput()->Set( itk::NumericTraits< RealType >::max() ); - this->GetVarianceOutput()->Set( itk::NumericTraits< RealType >::max() ); - this->GetSumOutput()->Set(itk::NumericTraits< RealType >::Zero); - this->GetCountOutput()->Set(itk::NumericTraits< LongType >::Zero); - - m_IgnoredInputValue = itk::NumericTraits<PixelType>::max(); -} - -template< typename TInputImage > -itk::DataObject::Pointer -StatisticsImageCustomFilter< TInputImage > -::MakeOutput(DataObjectPointerArraySizeType output) -{ - switch ( output ) - { - case 0: - return TInputImage::New().GetPointer(); - break; - case 1: - return PixelObjectType::New().GetPointer(); - break; - case 2: - return PixelObjectType::New().GetPointer(); - break; - case 3: - case 4: - case 5: - case 6: - return RealObjectType::New().GetPointer(); - break; - case 7: - return LongObjectType::New().GetPointer(); - break; - default: - // might as well make an image - return TInputImage::New().GetPointer(); - break; - } -} - -template< typename TInputImage > -typename StatisticsImageCustomFilter< TInputImage >::PixelObjectType * -StatisticsImageCustomFilter< TInputImage > -::GetMinimumOutput() -{ - return static_cast< PixelObjectType * >( this->itk::ProcessObject::GetOutput(1) ); -} - -template< typename TInputImage > -const typename StatisticsImageCustomFilter< TInputImage >::PixelObjectType * -StatisticsImageCustomFilter< TInputImage > -::GetMinimumOutput() const -{ - return static_cast< const PixelObjectType * >( this->itk::ProcessObject::GetOutput(1) ); -} - -template< typename TInputImage > -typename StatisticsImageCustomFilter< TInputImage >::PixelObjectType * -StatisticsImageCustomFilter< TInputImage > -::GetMaximumOutput() -{ - return static_cast< PixelObjectType * >( this->itk::ProcessObject::GetOutput(2) ); -} - -template< typename TInputImage > -const typename StatisticsImageCustomFilter< TInputImage >::PixelObjectType * -StatisticsImageCustomFilter< TInputImage > -::GetMaximumOutput() const -{ - return static_cast< const PixelObjectType * >( this->itk::ProcessObject::GetOutput(2) ); -} - -template< typename TInputImage > -typename StatisticsImageCustomFilter< TInputImage >::RealObjectType * -StatisticsImageCustomFilter< TInputImage > -::GetMeanOutput() -{ - return static_cast< RealObjectType * >( this->itk::ProcessObject::GetOutput(3) ); -} - -template< typename TInputImage > -const typename StatisticsImageCustomFilter< TInputImage >::RealObjectType * -StatisticsImageCustomFilter< TInputImage > -::GetMeanOutput() const -{ - return static_cast< const RealObjectType * >( this->itk::ProcessObject::GetOutput(3) ); -} - -template< typename TInputImage > -typename StatisticsImageCustomFilter< TInputImage >::RealObjectType * -StatisticsImageCustomFilter< TInputImage > -::GetSigmaOutput() -{ - return static_cast< RealObjectType * >( this->itk::ProcessObject::GetOutput(4) ); -} - -template< typename TInputImage > -const typename StatisticsImageCustomFilter< TInputImage >::RealObjectType * -StatisticsImageCustomFilter< TInputImage > -::GetSigmaOutput() const -{ - return static_cast< const RealObjectType * >( this->itk::ProcessObject::GetOutput(4) ); -} - -template< typename TInputImage > -typename StatisticsImageCustomFilter< TInputImage >::RealObjectType * -StatisticsImageCustomFilter< TInputImage > -::GetVarianceOutput() -{ - return static_cast< RealObjectType * >( this->itk::ProcessObject::GetOutput(5) ); -} - -template< typename TInputImage > -const typename StatisticsImageCustomFilter< TInputImage >::RealObjectType * -StatisticsImageCustomFilter< TInputImage > -::GetVarianceOutput() const -{ - return static_cast< const RealObjectType * >( this->itk::ProcessObject::GetOutput(5) ); -} - -template< typename TInputImage > -typename StatisticsImageCustomFilter< TInputImage >::RealObjectType * -StatisticsImageCustomFilter< TInputImage > -::GetSumOutput() -{ - return static_cast< RealObjectType * >( this->itk::ProcessObject::GetOutput(6) ); -} - -template< typename TInputImage > -const typename StatisticsImageCustomFilter< TInputImage >::RealObjectType * -StatisticsImageCustomFilter< TInputImage > -::GetSumOutput() const -{ - return static_cast< const RealObjectType * >( this->itk::ProcessObject::GetOutput(6) ); -} - -template< typename TInputImage > -typename StatisticsImageCustomFilter< TInputImage >::RealObjectType * -StatisticsImageCustomFilter< TInputImage > -::GetSumOfSquaresOutput() -{ - return static_cast< RealObjectType * >( this->itk::ProcessObject::GetOutput(6) ); -} - -template< typename TInputImage > -const typename StatisticsImageCustomFilter< TInputImage >::RealObjectType * -StatisticsImageCustomFilter< TInputImage > -::GetSumOfSquaresOutput() const -{ - return static_cast< const RealObjectType * >( this->itk::ProcessObject::GetOutput(6) ); -} - -template< typename TInputImage > -typename StatisticsImageCustomFilter< TInputImage >::LongObjectType * -StatisticsImageCustomFilter< TInputImage > -::GetCountOutput() -{ - return static_cast< LongObjectType * >( this->itk::ProcessObject::GetOutput(7) ); -} - -template< typename TInputImage > -const typename StatisticsImageCustomFilter< TInputImage >::LongObjectType * -StatisticsImageCustomFilter< TInputImage > -::GetCountOutput() const -{ - return static_cast< const LongObjectType * >( this->itk::ProcessObject::GetOutput(7) ); -} - -template< typename TInputImage > -void -StatisticsImageCustomFilter< 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< typename TInputImage > -void -StatisticsImageCustomFilter< TInputImage > -::BeforeThreadedGenerateData() -{ - itk::ThreadIdType 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(itk::NumericTraits< itk::SizeValueType >::Zero); - m_ThreadSum.Fill(itk::NumericTraits< RealType >::Zero); - m_SumOfSquares.Fill(itk::NumericTraits< RealType >::Zero); - m_ThreadMin.Fill( itk::NumericTraits< PixelType >::max() ); - m_ThreadMax.Fill( itk::NumericTraits< PixelType >::NonpositiveMin() ); -} - -template< typename TInputImage > -void -StatisticsImageCustomFilter< TInputImage > -::AfterThreadedGenerateData() -{ - itk::ThreadIdType i; - - itk::ThreadIdType numberOfThreads = this->GetNumberOfThreads(); - - PixelType minimum = this->GetMinimumOutput()->Get(); - PixelType maximum = this->GetMaximumOutput()->Get(); - RealType sum = this ->GetSumOutput()->Get(); - RealType sumOfSquares = this->GetSumOfSquaresOutput()->Get(); - LongType count = this->GetCountOutput()->Get(); - - RealType mean, sigma, variance; - mean = sigma = variance = itk::NumericTraits< RealType >::ZeroValue(); - - // Find the min/max over all threads and accumulate count, sum and - // sum of squares - 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 - if (count > 1) - { - mean = sum / static_cast< RealType >( count ); - - // unbiased estimate - variance = ( sumOfSquares - ( sum * sum / static_cast< RealType >( count ) ) ) - / ( static_cast< RealType >( count ) - 1 ); - sigma = std::sqrt(variance); - } - - // Update 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); - this->GetSumOfSquaresOutput()->Set(sumOfSquares); - this->GetCountOutput()->Set(count); -} - -template< typename TInputImage > -void -StatisticsImageCustomFilter< TInputImage > -::ThreadedGenerateData(const RegionType & outputRegionForThread, - itk::ThreadIdType threadId) -{ - const itk::SizeValueType size0 = outputRegionForThread.GetSize(0); - if( size0 == 0) - { - return; - } - RealType realValue; - PixelType value; - - RealType sum = itk::NumericTraits< RealType >::ZeroValue(); - RealType sumOfSquares = itk::NumericTraits< RealType >::ZeroValue(); - itk::SizeValueType count = itk::NumericTraits< itk::SizeValueType >::ZeroValue(); - PixelType min = itk::NumericTraits< PixelType >::max(); - PixelType max = itk::NumericTraits< PixelType >::NonpositiveMin(); - - itk::ImageScanlineConstIterator< TInputImage > it (this->GetInput(), outputRegionForThread); - - // support progress methods/callbacks - const size_t numberOfLinesToProcess = outputRegionForThread.GetNumberOfPixels() / size0; - itk::ProgressReporter progress( this, threadId, numberOfLinesToProcess ); - - // do the work - while ( !it.IsAtEnd() ) - { - while ( !it.IsAtEndOfLine() ) - { - value = it.Get(); - if (value != m_IgnoredInputValue) - { - realValue = static_cast< RealType >( value ); - if ( value < min ) - { - min = value; - } - if ( value > max ) - { - max = value; - } - sum += realValue; - sumOfSquares += ( realValue * realValue ); - ++count; - } - ++it; - } - it.NextLine(); - progress.CompletedPixel(); - } - - m_ThreadSum[threadId] = sum; - m_SumOfSquares[threadId] = sumOfSquares; - m_Count[threadId] = count; - m_ThreadMin[threadId] = min; - m_ThreadMax[threadId] = max; -} - -template< typename TImage > -void -StatisticsImageCustomFilter< TImage > -::PrintSelf(std::ostream & os, itk::Indent indent) const -{ - Superclass::PrintSelf(os, indent); - - os << indent << "Minimum: " - << static_cast< typename itk::NumericTraits< PixelType >::PrintType >( this->GetMinimum() ) << std::endl; - os << indent << "Maximum: " - << static_cast< typename itk::NumericTraits< PixelType >::PrintType >( this->GetMaximum() ) << std::endl; - os << indent << "Sum: " << this->GetSum() << std::endl; - os << indent << "SumOfSquares: " << this->GetSumOfSquares() << std::endl; - os << indent << "Count: " << this->GetCount() << 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