diff --git a/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.h b/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.h
index 3caeb00321d67454dabcf32a03f495a429cc2ac7..371ff0dffa6750142a066ddca9761f918fec0858 100644
--- a/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.h
+++ b/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.h
@@ -33,8 +33,8 @@
namespace otb
{
-/** \class PersistentStatisticsVectorImageFilter
- * \brief Compute min. max, covariance of a large image using streaming
+/** \class PersistentStreamingStatisticsVectorImageFilter
+ * \brief Compute covariance & correlation of a large image using streaming
*
* This filter persists its temporary data. It means that if you Update it n times on n different
* requested regions, the output statistics will be the statitics of the whole set of n regions.
@@ -49,13 +49,13 @@ namespace otb
* \ingroup MathematicalStatisticsImageFilters
*
*/
-template<class TInputImage>
-class ITK_EXPORT PersistentStatisticsVectorImageFilter :
+template<class TInputImage, class TPrecision >
+class ITK_EXPORT PersistentStreamingStatisticsVectorImageFilter :
public PersistentImageFilter<TInputImage, TInputImage>
{
public:
/** Standard Self typedef */
- typedef PersistentStatisticsVectorImageFilter Self;
+ typedef PersistentStreamingStatisticsVectorImageFilter Self;
typedef PersistentImageFilter<TInputImage, TInputImage> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
@@ -64,56 +64,53 @@ public:
itkNewMacro(Self);
/** Runtime information support. */
- itkTypeMacro(PersistentStatisticsVectorImageFilter, PersistentImageFilter);
+ itkTypeMacro(PersistentStreamingStatisticsVectorImageFilter, PersistentImageFilter);
/** Image related typedefs. */
- typedef TInputImage ImageType;
- 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;
- typedef typename TInputImage::InternalPixelType InternalPixelType;
+ typedef TInputImage ImageType;
+ typedef typename ImageType::Pointer InputImagePointer;
+ typedef typename ImageType::RegionType RegionType;
+ typedef typename ImageType::SizeType SizeType;
+ typedef typename ImageType::IndexType IndexType;
+ typedef typename ImageType::PixelType PixelType;
+ typedef typename ImageType::InternalPixelType InternalPixelType;
- itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension);
+ typedef TPrecision PrecisionType;
+ typedef PrecisionType RealType;
/** Image related typedefs. */
itkStaticConstMacro(ImageDimension, unsigned int, TInputImage::ImageDimension);
- /** Type to use for computations. */
- typedef typename itk::NumericTraits<InternalPixelType>::RealType RealType;
- typedef itk::VariableLengthVector<RealType> RealPixelType;
-
/** Smart Pointer type to a DataObject. */
typedef typename itk::DataObject::Pointer DataObjectPointer;
- /** Type of DataObjects used for scalar outputs */
- typedef typename itk::VariableSizeMatrix<RealType> MatrixType;
- typedef typename std::vector<MatrixType> ArrayMatrixType;
- typedef typename itk::Array<long> ArrayLongPixelType;
- typedef typename std::vector<RealPixelType> ArrayRealPixelType;
- typedef typename std::vector<PixelType> ArrayPixelType;
- typedef itk::SimpleDataObjectDecorator<RealPixelType> RealPixelObjectType;
+ /** Type to use for computations. */
+ typedef itk::VariableSizeMatrix<PrecisionType> MatrixType;
+ typedef itk::VariableLengthVector<PrecisionType> RealPixelType;
+
+ /** Type of DataObjects used for outputs */
+ typedef itk::SimpleDataObjectDecorator<IndexType> IndexObjectType;
typedef itk::SimpleDataObjectDecorator<PixelType> PixelObjectType;
+ typedef itk::SimpleDataObjectDecorator<RealPixelType> RealPixelObjectType;
typedef itk::SimpleDataObjectDecorator<MatrixType> MatrixObjectType;
- /** Return the computed Minimum. */
+ /** Return the computed Min */
PixelType GetMinimum() const
{
- return this->GetMinimumOutput()->Get();
+ return this->GetMinOutput()->Get();
}
PixelObjectType* GetMinimumOutput();
const PixelObjectType* GetMinimumOutput() const;
- /** Return the computed Maximum. */
+ /** Return the computed Max index */
PixelType GetMaximum() const
{
- return this->GetMaximumOutput()->Get();
+ return this->GetMaxOutput()->Get();
}
PixelObjectType* GetMaximumOutput();
const PixelObjectType* GetMaximumOutput() const;
- /** Return the computed Mean. */
+ /** Return the computed Mean. */
RealPixelType GetMean() const
{
return this->GetMeanOutput()->Get();
@@ -121,7 +118,7 @@ public:
RealPixelObjectType* GetMeanOutput();
const RealPixelObjectType* GetMeanOutput() const;
- /** Return the compute Sum. */
+ /** Return the computed Sum. */
RealPixelType GetSum() const
{
return this->GetSumOutput()->Get();
@@ -129,6 +126,14 @@ public:
RealPixelObjectType* GetSumOutput();
const RealPixelObjectType* GetSumOutput() const;
+ /** Return the computed Covariance. */
+ MatrixType GetCorrelation() const
+ {
+ return this->GetCorrelation()->Get();
+ }
+ MatrixObjectType* GetCorrelationOutput();
+ const MatrixObjectType* GetCorrelationOutput() const;
+
/** Return the computed Covariance. */
MatrixType GetCovariance() const
{
@@ -142,45 +147,64 @@ public:
*/
virtual DataObjectPointer MakeOutput(unsigned int idx);
+ virtual void Reset(void);
+
+ virtual void Synthetize(void);
+
+ itkSetMacro(EnableMinMax, bool);
+ itkGetMacro(EnableMinMax, bool);
+
+ itkSetMacro(EnableFirstOrderStats, bool);
+ itkGetMacro(EnableFirstOrderStats, bool);
+
+ itkSetMacro(EnableSecondOrderStats, bool);
+ itkGetMacro(EnableSecondOrderStats, bool);
+
+protected:
+ PersistentStreamingStatisticsVectorImageFilter();
+
+ virtual ~PersistentStreamingStatisticsVectorImageFilter() {}
+
/** Pass the input through unmodified. Do this by Grafting in the
* AllocateOutputs method.
*/
virtual void AllocateOutputs();
+
virtual void GenerateOutputInformation();
- virtual void Synthetize(void);
- virtual void Reset(void);
-protected:
- PersistentStatisticsVectorImageFilter();
- virtual ~PersistentStatisticsVectorImageFilter() {}
+
virtual void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
/** Multi-thread version GenerateData. */
void ThreadedGenerateData(const RegionType& outputRegionForThread, int threadId);
private:
- PersistentStatisticsVectorImageFilter(const Self &); //purposely not implemented
+ PersistentStreamingStatisticsVectorImageFilter(const Self &); //purposely not implemented
void operator =(const Self&); //purposely not implemented
- ArrayRealPixelType m_ThreadSum;
- ArrayLongPixelType m_Count;
- ArrayPixelType m_ThreadMin;
- ArrayPixelType m_ThreadMax;
- ArrayMatrixType m_XX;
+ bool m_EnableMinMax;
+ bool m_EnableFirstOrderStats;
+ bool m_EnableSecondOrderStats;
+
+ std::vector<PixelType> m_ThreadMin;
+ std::vector<PixelType> m_ThreadMax;
+ std::vector<RealPixelType> m_ThreadFirstOrderAccumulators;
+ std::vector<MatrixType> m_ThreadSecondOrderAccumulators;
-}; // end of class PersistentStatisticsVectorImageFilter
+}; // end of class PersistentStreamingStatisticsVectorImageFilter
/**===========================================================================*/
/** \class StreamingStatisticsVectorImageFilter
* \brief This class streams the whole input image through the PersistentStatisticsImageFilter.
*
- * This way, it allows to compute the first order global statistics of this image. It calls the
- * Reset() method of the PersistentStatisticsImageFilter before streaming the image and the
- * Synthetize() method of the PersistentStatisticsImageFilter after having streamed the image
+ * This way, it allows to compute the first and second order global statistics of this image. It calls the
+ * Reset() method of the PersistentStreamingStatisticsVectorImageFilter before streaming the image and the
+ * Synthetize() method of the PersistentStreamingStatisticsVectorImageFilter after having streamed the image
* to compute the statistics. The accessor on the results are wrapping the accessors of the
- * internal PersistentStatisticsImageFilter.
+ * internal PersistentStreamingStatisticsVectorImageFilter.
*
- * \sa PersistentStatisticsVectorImageFilter
+ * \sa PersistentStreamingStatisticsVectorImageFilter
* \sa PersistentImageFilter
* \sa PersistentFilterStreamingDecorator
* \sa StreamingImageVirtualWriter
@@ -189,15 +213,15 @@ private:
* \ingroup MathematicalStatisticsImageFilters
*/
-template<class TInputImage>
+template<class TInputImage, class TPrecision = typename itk::NumericTraits<typename TInputImage::InternalPixelType>::RealType>
class ITK_EXPORT StreamingStatisticsVectorImageFilter :
- public PersistentFilterStreamingDecorator<PersistentStatisticsVectorImageFilter<TInputImage> >
+ public PersistentFilterStreamingDecorator<PersistentStreamingStatisticsVectorImageFilter<TInputImage, TPrecision> >
{
public:
/** Standard Self typedef */
typedef StreamingStatisticsVectorImageFilter Self;
typedef PersistentFilterStreamingDecorator
- <PersistentStatisticsVectorImageFilter<TInputImage> > Superclass;
+ <PersistentStreamingStatisticsVectorImageFilter<TInputImage, TPrecision> > Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
@@ -209,18 +233,13 @@ public:
typedef TInputImage InputImageType;
typedef typename Superclass::FilterType StatFilterType;
- typedef typename StatFilterType::PixelType PixelType;
- typedef typename StatFilterType::RealType RealType;
- typedef typename StatFilterType::RealPixelType RealPixelType;
- typedef typename StatFilterType::MatrixType MatrixType;
- typedef typename StatFilterType::ArrayMatrixType ArrayMatrixType;
- typedef typename StatFilterType::ArrayLongPixelType ArrayLongPixelType;
- typedef typename StatFilterType::ArrayRealPixelType ArrayRealPixelType;
- typedef typename StatFilterType::ArrayPixelType ArrayPixelType;
-
- /** Type of DataObjects used for scalar outputs */
+
+ /** Type of DataObjects used for outputs */
+ typedef typename StatFilterType::PixelType PixelType;
+ typedef typename StatFilterType::RealType RealType;
+ typedef typename StatFilterType::RealPixelType RealPixelType;
typedef typename StatFilterType::RealPixelObjectType RealPixelObjectType;
- typedef typename StatFilterType::PixelObjectType PixelObjectType;
+ typedef typename StatFilterType::MatrixType MatrixType;
typedef typename StatFilterType::MatrixObjectType MatrixObjectType;
void SetInput(InputImageType * input)
@@ -232,30 +251,30 @@ public:
return this->GetFilter()->GetInput();
}
- /** Return the computed Minimum. */
- PixelType GetMinimum() const
+ /** Return the computed Mean. */
+ RealPixelType GetMinimum() const
{
return this->GetFilter()->GetMinimumOutput()->Get();
}
- PixelObjectType* GetMinimumOutput()
+ RealPixelObjectType* GetMinimumOutput()
{
return this->GetFilter()->GetMinimumOutput();
}
- const PixelObjectType* GetMinimumOutput() const
+ const RealPixelObjectType* GetMinimumOutput() const
{
return this->GetFilter()->GetMinimumOutput();
}
- /** Return the computed Maximum. */
- PixelType GetMaximum() const
+ /** Return the computed Mean. */
+ RealPixelType GetMaximum() const
{
return this->GetFilter()->GetMaximumOutput()->Get();
}
- PixelObjectType* GetMaximumOutput()
+ RealPixelObjectType* GetMaximumOutput()
{
return this->GetFilter()->GetMaximumOutput();
}
- const PixelObjectType* GetMaximumOutput() const
+ const RealPixelObjectType* GetMaximumOutput() const
{
return this->GetFilter()->GetMaximumOutput();
}
@@ -274,7 +293,7 @@ public:
return this->GetFilter()->GetMeanOutput();
}
- /** Return the compute Sum. */
+ /** Return the computed Mean. */
RealPixelType GetSum() const
{
return this->GetFilter()->GetSumOutput()->Get();
@@ -302,9 +321,33 @@ public:
return this->GetFilter()->GetCovarianceOutput();
}
+ /** Return the computed Covariance. */
+ MatrixType GetCorrelation() const
+ {
+ return this->GetFilter()->GetCorrelationOutput()->Get();
+ }
+ MatrixObjectType* GetCorrelationOutput()
+ {
+ return this->GetFilter()->GetCorrelationOutput();
+ }
+ const MatrixObjectType* GetCorrelationOutput() const
+ {
+ return this->GetFilter()->GetCorrelationOutput();
+ }
+
+ otbSetObjectMemberMacro(Filter, EnableMinMax, bool);
+ otbGetObjectMemberMacro(Filter, EnableMinMax, bool);
+
+ otbSetObjectMemberMacro(Filter, EnableFirstOrderStats, bool);
+ otbGetObjectMemberMacro(Filter, EnableFirstOrderStats, bool);
+
+ otbSetObjectMemberMacro(Filter, EnableSecondOrderStats, bool);
+ otbGetObjectMemberMacro(Filter, EnableSecondOrderStats, bool);
+
protected:
/** Constructor */
- StreamingStatisticsVectorImageFilter() {};
+ StreamingStatisticsVectorImageFilter() {}
+
/** Destructor */
virtual ~StreamingStatisticsVectorImageFilter() {}
diff --git a/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.txx b/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.txx
index 6dbd12acabe9001ca11e25c036c7befcca9e0a5c..c576cef3fd2a8f540a383758a2adac4360fa1ade 100644
--- a/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.txx
+++ b/Code/BasicFilters/otbStreamingStatisticsVectorImageFilter.txx
@@ -31,40 +31,27 @@ for details.
namespace otb
{
-template<class TInputImage>
-PersistentStatisticsVectorImageFilter<TInputImage>
-::PersistentStatisticsVectorImageFilter()
+template<class TInputImage, class TPrecision>
+PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>
+::PersistentStreamingStatisticsVectorImageFilter()
+ : m_EnableMinMax(true),
+ m_EnableFirstOrderStats(true),
+ m_EnableSecondOrderStats(true)
{
// 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 < 5; ++i)
+ // just decorators around vector/matrix types
+ for (unsigned int i = 1; i < 7; ++i)
{
- typename RealPixelObjectType::Pointer output = static_cast<RealPixelObjectType*>(this->MakeOutput(i).GetPointer());
- this->itk::ProcessObject::SetNthOutput(i, output.GetPointer());
+ this->itk::ProcessObject::SetNthOutput(i, this->MakeOutput(i).GetPointer());
}
-
- // allocate the data objects for the outputs which are
- // just decorators around matrix types
- typename MatrixObjectType::Pointer output = static_cast<MatrixObjectType*>(this->MakeOutput(5).GetPointer());
- this->itk::ProcessObject::SetNthOutput(5, output.GetPointer());
}
-template<class TInputImage>
+template<class TInputImage, class TPrecision>
itk::DataObject::Pointer
-PersistentStatisticsVectorImageFilter<TInputImage>
+PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>
::MakeOutput(unsigned int output)
{
switch (output)
@@ -74,13 +61,17 @@ PersistentStatisticsVectorImageFilter<TInputImage>
break;
case 1:
case 2:
+ // min/max
return static_cast<itk::DataObject*>(PixelObjectType::New().GetPointer());
break;
case 3:
case 4:
+ // mean / sum
return static_cast<itk::DataObject*>(RealPixelObjectType::New().GetPointer());
break;
case 5:
+ case 6:
+ // covariance / correlation
return static_cast<itk::DataObject*>(MatrixObjectType::New().GetPointer());
break;
default:
@@ -88,92 +79,107 @@ PersistentStatisticsVectorImageFilter<TInputImage>
return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
break;
}
-
}
-template<class TInputImage>
-typename PersistentStatisticsVectorImageFilter<TInputImage>::PixelObjectType*
-PersistentStatisticsVectorImageFilter<TInputImage>
+template<class TInputImage, class TPrecision>
+typename PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>::PixelObjectType*
+PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>
::GetMinimumOutput()
{
return static_cast<PixelObjectType*>(this->itk::ProcessObject::GetOutput(1));
}
-template<class TInputImage>
-const typename PersistentStatisticsVectorImageFilter<TInputImage>::PixelObjectType*
-PersistentStatisticsVectorImageFilter<TInputImage>
+template<class TInputImage, class TPrecision>
+const typename PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>::PixelObjectType*
+PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>
::GetMinimumOutput() const
{
return static_cast<const PixelObjectType*>(this->itk::ProcessObject::GetOutput(1));
}
-template<class TInputImage>
-typename PersistentStatisticsVectorImageFilter<TInputImage>::PixelObjectType*
-PersistentStatisticsVectorImageFilter<TInputImage>
+template<class TInputImage, class TPrecision>
+typename PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>::PixelObjectType*
+PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>
::GetMaximumOutput()
{
return static_cast<PixelObjectType*>(this->itk::ProcessObject::GetOutput(2));
}
-template<class TInputImage>
-const typename PersistentStatisticsVectorImageFilter<TInputImage>::PixelObjectType*
-PersistentStatisticsVectorImageFilter<TInputImage>
+template<class TInputImage, class TPrecision>
+const typename PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>::PixelObjectType*
+PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>
::GetMaximumOutput() const
{
return static_cast<const PixelObjectType*>(this->itk::ProcessObject::GetOutput(2));
}
-template<class TInputImage>
-typename PersistentStatisticsVectorImageFilter<TInputImage>::RealPixelObjectType*
-PersistentStatisticsVectorImageFilter<TInputImage>
+template<class TInputImage, class TPrecision>
+typename PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>::RealPixelObjectType*
+PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>
::GetMeanOutput()
{
return static_cast<RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(3));
}
-template<class TInputImage>
-const typename PersistentStatisticsVectorImageFilter<TInputImage>::RealPixelObjectType*
-PersistentStatisticsVectorImageFilter<TInputImage>
+template<class TInputImage, class TPrecision>
+const typename PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>::RealPixelObjectType*
+PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>
::GetMeanOutput() const
{
return static_cast<const RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(3));
}
-template<class TInputImage>
-typename PersistentStatisticsVectorImageFilter<TInputImage>::RealPixelObjectType*
-PersistentStatisticsVectorImageFilter<TInputImage>
+template<class TInputImage, class TPrecision>
+typename PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>::RealPixelObjectType*
+PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>
::GetSumOutput()
{
return static_cast<RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(4));
}
-template<class TInputImage>
-const typename PersistentStatisticsVectorImageFilter<TInputImage>::RealPixelObjectType*
-PersistentStatisticsVectorImageFilter<TInputImage>
+template<class TInputImage, class TPrecision>
+const typename PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>::RealPixelObjectType*
+PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>
::GetSumOutput() const
{
return static_cast<const RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(4));
}
-template<class TInputImage>
-typename PersistentStatisticsVectorImageFilter<TInputImage>::MatrixObjectType*
-PersistentStatisticsVectorImageFilter<TInputImage>
-::GetCovarianceOutput()
+template<class TInputImage, class TPrecision>
+typename PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>::MatrixObjectType*
+PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>
+::GetCorrelationOutput()
{
return static_cast<MatrixObjectType*>(this->itk::ProcessObject::GetOutput(5));
}
-template<class TInputImage>
-const typename PersistentStatisticsVectorImageFilter<TInputImage>::MatrixObjectType*
-PersistentStatisticsVectorImageFilter<TInputImage>
-::GetCovarianceOutput() const
+template<class TInputImage, class TPrecision>
+const typename PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>::MatrixObjectType*
+PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>
+::GetCorrelationOutput() const
{
return static_cast<const MatrixObjectType*>(this->itk::ProcessObject::GetOutput(5));
}
-template<class TInputImage>
+template<class TInputImage, class TPrecision>
+typename PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>::MatrixObjectType*
+PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>
+::GetCovarianceOutput()
+{
+ return static_cast<MatrixObjectType*>(this->itk::ProcessObject::GetOutput(6));
+}
+
+template<class TInputImage, class TPrecision>
+const typename PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>::MatrixObjectType*
+PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>
+::GetCovarianceOutput() const
+{
+ return static_cast<const MatrixObjectType*>(this->itk::ProcessObject::GetOutput(6));
+}
+
+template<class TInputImage, class TPrecision>
void
-PersistentStatisticsVectorImageFilter<TInputImage>
+PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>
::GenerateOutputInformation()
{
Superclass::GenerateOutputInformation();
@@ -189,9 +195,9 @@ PersistentStatisticsVectorImageFilter<TInputImage>
}
}
-template<class TInputImage>
+template<class TInputImage, class TPrecision>
void
-PersistentStatisticsVectorImageFilter<TInputImage>
+PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>
::AllocateOutputs()
{
// This is commented to prevent the streaming of the whole image for the first stream strip
@@ -201,9 +207,9 @@ PersistentStatisticsVectorImageFilter<TInputImage>
// Nothing that needs to be allocated for the remaining outputs
}
-template<class TInputImage>
+template<class TInputImage, class TPrecision>
void
-PersistentStatisticsVectorImageFilter<TInputImage>
+PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>
::Reset()
{
TInputImage * inputPtr = const_cast<TInputImage *>(this->GetInput());
@@ -212,255 +218,206 @@ PersistentStatisticsVectorImageFilter<TInputImage>
unsigned int numberOfThreads = this->GetNumberOfThreads();
unsigned int numberOfComponent = inputPtr->GetNumberOfComponentsPerPixel();
- // Variable Initialization
- PixelType tempPixel;
- tempPixel.SetSize(numberOfComponent);
- tempPixel.Fill(itk::NumericTraits<InternalPixelType>::NonpositiveMin());
- this->GetMaximumOutput()->Set(tempPixel);
-
- tempPixel.Fill(itk::NumericTraits<InternalPixelType>::max());
- this->GetMinimumOutput()->Set(tempPixel);
+ if (m_EnableMinMax)
+ {
+ PixelType tempPixel;
+ tempPixel.SetSize(numberOfComponent);
- RealPixelType tempRealPixel;
- tempRealPixel.SetSize(numberOfComponent);
- tempRealPixel.Fill(itk::NumericTraits<RealType>::max());
- this->GetMeanOutput()->Set(tempRealPixel);
+ tempPixel.Fill(itk::NumericTraits<InternalPixelType>::max());
+ this->GetMinimumOutput()->Set(tempPixel);
- tempRealPixel.Fill(itk::NumericTraits<RealType>::Zero);
- this->GetMeanOutput()->Set(tempRealPixel);
+ tempPixel.Fill(itk::NumericTraits<InternalPixelType>::NonpositiveMin());
+ this->GetMaximumOutput()->Set(tempPixel);
- MatrixType tempMatrix;
- tempMatrix.SetSize(numberOfComponent, numberOfComponent);
- tempMatrix.Fill(itk::NumericTraits<RealType>::Zero);
+ PixelType tempTemporiesPixel;
+ tempTemporiesPixel.SetSize(numberOfComponent);
+ tempTemporiesPixel.Fill(itk::NumericTraits<InternalPixelType>::max());
+ m_ThreadMin = std::vector<PixelType>(numberOfThreads, tempTemporiesPixel);
- // Initialize the tempories
- m_Count.SetSize(numberOfThreads);
- m_Count.Fill(itk::NumericTraits<long>::Zero);
+ tempTemporiesPixel.Fill(itk::NumericTraits<InternalPixelType>::NonpositiveMin());
+ m_ThreadMax = std::vector<PixelType>(numberOfThreads, tempTemporiesPixel);
+ }
- PixelType tempTemporiesPixel;
- tempTemporiesPixel.SetSize(numberOfComponent);
- tempTemporiesPixel.Fill(itk::NumericTraits<InternalPixelType>::max());
- m_ThreadMin = ArrayPixelType(numberOfThreads, tempTemporiesPixel);
+ if (m_EnableSecondOrderStats)
+ {
+ m_EnableFirstOrderStats = true;
+ }
- tempTemporiesPixel.Fill(itk::NumericTraits<InternalPixelType>::NonpositiveMin());
- m_ThreadMax = ArrayPixelType(numberOfThreads, tempTemporiesPixel);
+ if (m_EnableFirstOrderStats)
+ {
+ RealPixelType zeroRealPixel;
+ zeroRealPixel.SetSize(numberOfComponent);
+ zeroRealPixel.Fill(itk::NumericTraits<PrecisionType>::Zero);
+ this->GetMeanOutput()->Set(zeroRealPixel);
+ this->GetSumOutput()->Set(zeroRealPixel);
+
+ m_ThreadFirstOrderAccumulators.resize(numberOfThreads);
+ std::fill(m_ThreadFirstOrderAccumulators.begin(), m_ThreadFirstOrderAccumulators.end(), zeroRealPixel);
+ }
- RealPixelType tempTemporiesRealPixel;
- tempTemporiesRealPixel.SetSize(numberOfComponent);
- tempTemporiesRealPixel.Fill(itk::NumericTraits<RealType>::Zero);
- m_ThreadSum = ArrayRealPixelType(numberOfThreads, tempTemporiesRealPixel);
+ if (m_EnableSecondOrderStats)
+ {
+ MatrixType zeroMatrix;
+ zeroMatrix.SetSize(numberOfComponent, numberOfComponent);
+ zeroMatrix.Fill(itk::NumericTraits<PrecisionType>::Zero);
+ this->GetCovarianceOutput()->Set(zeroMatrix);
+ this->GetCorrelationOutput()->Set(zeroMatrix);
+
+ m_ThreadSecondOrderAccumulators.resize(numberOfThreads);
+ std::fill(m_ThreadSecondOrderAccumulators.begin(), m_ThreadSecondOrderAccumulators.end(), zeroMatrix);
+ }
- m_XX = ArrayMatrixType(numberOfThreads, tempMatrix);
}
-template<class TInputImage>
+template<class TInputImage, class TPrecision>
void
-PersistentStatisticsVectorImageFilter<TInputImage>
+PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>
::Synthetize()
{
- int i;
- long count;
-
- int numberOfThreads = this->GetNumberOfThreads();
- unsigned int numberOfComponent = this->GetInput()->GetNumberOfComponentsPerPixel();
-
- PixelType minimumVector;
- minimumVector.SetSize(numberOfComponent);
- minimumVector.Fill(itk::NumericTraits<InternalPixelType>::max());
-
- PixelType maximumVector;
- maximumVector.SetSize(numberOfComponent);
- maximumVector.Fill(itk::NumericTraits<InternalPixelType>::NonpositiveMin());
-
- RealPixelType sumVector;
- sumVector.SetSize(numberOfComponent);
- sumVector.Fill(itk::NumericTraits<RealType>::Zero);
-
- RealPixelType meanVector;
- meanVector.SetSize(numberOfComponent);
- MatrixType crossedMatrix;
- crossedMatrix.SetSize(numberOfComponent, numberOfComponent);
- crossedMatrix.Fill(itk::NumericTraits<RealType>::Zero);
- count = 0;
-
- // Find the min/max over all threads and accumulate count, sum and
- // sum of squares
- for (i = 0; i < numberOfThreads; ++i)
+ TInputImage * inputPtr = const_cast<TInputImage *>(this->GetInput());
+ const unsigned int nbPixels = inputPtr->GetLargestPossibleRegion().GetNumberOfPixels();
+ const unsigned int numberOfComponent = inputPtr->GetNumberOfComponentsPerPixel();
+
+ PixelType minimum;
+ minimum.SetSize(numberOfComponent);
+ minimum.Fill(itk::NumericTraits<InternalPixelType>::max());
+ PixelType maximum;
+ maximum.SetSize(numberOfComponent);
+ maximum.Fill(itk::NumericTraits<InternalPixelType>::NonpositiveMin());
+
+ RealPixelType streamFirstOrderAccumulator(numberOfComponent);
+ streamFirstOrderAccumulator.Fill(itk::NumericTraits<PrecisionType>::Zero);
+ MatrixType streamSecondOrderAccumulator(numberOfComponent, numberOfComponent);
+ streamSecondOrderAccumulator.Fill(itk::NumericTraits<PrecisionType>::Zero);
+
+ // Accumulate results from all threads
+ const unsigned int numberOfThreads = this->GetNumberOfThreads();
+ for (unsigned int threadId = 0; threadId < numberOfThreads; ++threadId)
{
- count += m_Count[i];
- /** TODO
- * To modify using + method operator. If we use it now -> exceptionmacro (no GetClassName...)
- * crossedMatrix += m_XX[i];
- **/
- if ((m_XX[i].Rows() != crossedMatrix.Rows()) || (m_XX[i].Cols() != crossedMatrix.Cols()))
+ if (m_EnableMinMax)
{
- itkExceptionMacro(<< "Matrix with size (" << m_XX[i].Rows() << "," <<
- m_XX[i].Cols() << ") cannot be subtracted from matrix with size (" <<
- crossedMatrix.Rows() << "," << crossedMatrix.Cols());
- }
+ const PixelType& threadMin = m_ThreadMin [threadId];
+ const PixelType& threadMax = m_ThreadMax [threadId];
- for (unsigned int r = 0; r < m_XX[i].Rows(); ++r)
- {
- for (unsigned int c = 0; c < m_XX[i].Cols(); ++c)
+ for (unsigned int j = 0; j < numberOfComponent; ++j)
{
- crossedMatrix(r, c) += m_XX[i](r, c);
+ if (threadMin[j] < minimum[j])
+ {
+ minimum[j] = threadMin[j];
+ }
+ if (threadMax[j] > maximum[j])
+ {
+ maximum[j] = threadMax[j];
+ }
}
}
- //**** END TODO ****
- for (unsigned int j = 0; j < numberOfComponent; ++j)
- {
- sumVector[j] += m_ThreadSum[i][j];
- if (m_ThreadMin[i][j] < minimumVector[j])
- {
- minimumVector[j] = m_ThreadMin[i][j];
- }
- if (m_ThreadMax[i][j] > maximumVector[j])
- {
- maximumVector[j] = m_ThreadMax[i][j];
- }
- }
- } // end for( i = 0; i < numberOfThreads; ++i)
+ if (m_EnableFirstOrderStats)
+ streamFirstOrderAccumulator += m_ThreadFirstOrderAccumulators [threadId];
- for (unsigned int j = 0; j < numberOfComponent; ++j)
- {
- // compute statistics
- meanVector[j] = sumVector[j] / static_cast<RealType>(count);
+ if (m_EnableSecondOrderStats)
+ streamSecondOrderAccumulator += m_ThreadSecondOrderAccumulators[threadId];
}
- // Compute Matrix Covariance
- MatrixType pixelSumMatrix;
- pixelSumMatrix.SetSize(static_cast<unsigned int>(numberOfComponent), 1);
- pixelSumMatrix.Fill(itk::NumericTraits<RealType>::Zero);
- for (unsigned int j = 0; j < numberOfComponent; ++j)
+ // Final calculations
+ if (m_EnableMinMax)
{
- pixelSumMatrix(j, 0) = sumVector[j];
+ this->GetMinimumOutput()->Set(minimum);
+ this->GetMaximumOutput()->Set(maximum);
}
- MatrixType covMatrix, covMatrixTemp, tempTranspose;
- covMatrix.SetSize(static_cast<unsigned int>(numberOfComponent), static_cast<unsigned int>(numberOfComponent));
- covMatrixTemp.SetSize(static_cast<unsigned int>(numberOfComponent), static_cast<unsigned int>(numberOfComponent));
- tempTranspose.SetSize(1, static_cast<unsigned int>(numberOfComponent));
-
- covMatrix = crossedMatrix / static_cast<RealType>(count);
- pixelSumMatrix /= static_cast<RealType>(count);
- tempTranspose = pixelSumMatrix.GetTranspose();
- covMatrixTemp = pixelSumMatrix * tempTranspose;
- /** TODO
- * To modify using - method operator. If we use it now -> exceptionmacro (no GetClassName...)
- *covMatrix -= covMatrixTemp;
- **/
- if ((covMatrix.Rows() != covMatrixTemp.Rows()) || (covMatrix.Cols() != covMatrixTemp.Cols()))
+ if (m_EnableFirstOrderStats)
{
- itkExceptionMacro(<< "Matrix with size (" << covMatrix.Rows() << "," <<
- covMatrix.Cols() << ") cannot be subtracted from matrix with size (" <<
- covMatrixTemp.Rows() << "," << covMatrixTemp.Cols());
+ this->GetMeanOutput()->Set(streamFirstOrderAccumulator / nbPixels);
+ this->GetSumOutput()->Set(streamFirstOrderAccumulator);
}
- for (unsigned int r = 0; r < covMatrix.Rows(); ++r)
+ if (m_EnableSecondOrderStats)
{
- for (unsigned int c = 0; c < covMatrix.Cols(); ++c)
+ MatrixType cor = streamSecondOrderAccumulator / nbPixels;
+ this->GetCorrelationOutput()->Set(cor);
+
+ const RealPixelType& mean = this->GetMeanOutput()->Get();
+ const double regul = static_cast<double>(nbPixels) / (nbPixels - 1);
+ MatrixType cov = cor;
+ for (unsigned int r = 0; r < numberOfComponent; ++r)
{
- covMatrix(r, c) -= covMatrixTemp(r, c);
+ for (unsigned int c = 0; c < numberOfComponent; ++c)
+ {
+ cov(r, c) = regul * (cov(r, c) - mean[r] * mean[c]);
+ }
}
+ this->GetCovarianceOutput()->Set(cov);
}
- //**** END TODO ****/
-
- // Set the outputs
- this->GetMinimumOutput()->Set(minimumVector);
- this->GetMaximumOutput()->Set(maximumVector);
- this->GetMeanOutput()->Set(meanVector);
- this->GetSumOutput()->Set(sumVector);
- this->GetCovarianceOutput()->Set(covMatrix);
}
-template<class TInputImage>
+template<class TInputImage, class TPrecision>
void
-PersistentStatisticsVectorImageFilter<TInputImage>
+PersistentStreamingStatisticsVectorImageFilter<TInputImage,TPrecision>
::ThreadedGenerateData(const RegionType& outputRegionForThread, int threadId)
{
- /**
- * Grab the input
- */
- InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput());
- // support progress methods/callbacks
+ // Support progress methods/callbacks
itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
- MatrixType pixelVector, pixelTransposeVector, pixelSumVector, tempMatrix;
-
- pixelVector.SetSize(this->GetInput()->GetNumberOfComponentsPerPixel(), 1);
- pixelVector.Fill(itk::NumericTraits<RealType>::Zero);
- pixelTransposeVector.SetSize(1, this->GetInput()->GetNumberOfComponentsPerPixel());
- pixelTransposeVector.Fill(itk::NumericTraits<RealType>::Zero);
- pixelSumVector.SetSize(this->GetInput()->GetNumberOfComponentsPerPixel(),
- this->GetInput()->GetNumberOfComponentsPerPixel());
- pixelSumVector.Fill(itk::NumericTraits<RealType>::Zero);
- tempMatrix.SetSize(this->GetInput()->GetNumberOfComponentsPerPixel(), this->GetInput()->GetNumberOfComponentsPerPixel());
+ // Grab the input
+ InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput());
+ PixelType& threadMin = m_ThreadMin [threadId];
+ PixelType& threadMax = m_ThreadMax [threadId];
+ RealPixelType& threadFirstOrder = m_ThreadFirstOrderAccumulators [threadId];
+ MatrixType& threadSecondOrder = m_ThreadSecondOrderAccumulators[threadId];
itk::ImageRegionConstIteratorWithIndex<TInputImage> it(inputPtr, outputRegionForThread);
- it.GoToBegin();
- // do the work
- while (!it.IsAtEnd())
+
+ for (it.GoToBegin(); !it.IsAtEnd(); ++it, progress.CompletedPixel())
{
- IndexType index = it.GetIndex();
- PixelType vectorValue = it.Get();
- for (unsigned int j = 0; j < vectorValue.GetSize(); ++j)
- {
- InternalPixelType value = vectorValue[j];
+ const PixelType& vectorValue = it.Get();
- RealType realValue = static_cast<RealType>(value);
- if (value < m_ThreadMin[threadId][j])
- {
- m_ThreadMin[threadId][j] = value;
- }
- if (value > m_ThreadMax[threadId][j])
+ if (m_EnableMinMax)
+ {
+ for (unsigned int j = 0; j < vectorValue.GetSize(); ++j)
{
- m_ThreadMax[threadId][j] = value;
+ if (vectorValue[j] < threadMin[j])
+ {
+ threadMin[j] = vectorValue[j];
+ }
+ if (vectorValue[j] > threadMax[j])
+ {
+ threadMax[j] = vectorValue[j];
+ }
}
- m_ThreadSum[threadId][j] += realValue;
- pixelVector(j, 0) = realValue;
}
- ++it;
- progress.CompletedPixel();
- pixelTransposeVector = pixelVector.GetTranspose();
- /** TODO
- * To modify using + method operator. If we use it now -> exceptionmacro (no GetClassName...)
- * m_XX[threadId]+=pixelVector*pixelTransposeVector;
- **/
- tempMatrix = pixelVector * pixelTransposeVector;
- if ((m_XX[threadId].Rows() != tempMatrix.Rows()) || (m_XX[threadId].Cols() != tempMatrix.Cols()))
+ if (m_EnableFirstOrderStats)
{
- itkExceptionMacro(<< "Matrix with size (" << m_XX[threadId].Rows() << "," <<
- m_XX[threadId].Cols() << ") cannot be subtracted from matrix with size (" <<
- tempMatrix.Rows() << "," << tempMatrix.Cols());
+ threadFirstOrder += vectorValue;
}
- for (unsigned int r = 0; r < m_XX[threadId].Rows(); ++r)
+ if (m_EnableSecondOrderStats)
{
- for (unsigned int c = 0; c < m_XX[threadId].Cols(); ++c)
+ for (unsigned int r = 0; r < threadSecondOrder.Rows(); ++r)
{
- m_XX[threadId](r, c) += tempMatrix(r, c);
+ for (unsigned int c = 0; c < threadSecondOrder.Cols(); ++c)
+ {
+ threadSecondOrder(r,c) += vectorValue[r] * vectorValue[c];
+ }
}
}
- //**** END TODO ****
- m_Count[threadId]++;
+
}
}
-template <class TImage>
+template <class TImage, class TPrecision>
void
-PersistentStatisticsVectorImageFilter<TImage>
+PersistentStreamingStatisticsVectorImageFilter<TImage,TPrecision>
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
Superclass::PrintSelf(os, indent);
- os << indent << "Minimum: " << this->GetMinimumOutput()->Get() << std::endl;
- os << indent << "Maximum: " << this->GetMaximumOutput()->Get() << std::endl;
- os << indent << "Sum: " << this->GetSumOutput()->Get() << std::endl;
- os << indent << "Mean: " << this->GetMeanOutput()->Get() << std::endl;
- os << indent << "Covariance: " << this->GetCovarianceOutput()->Get() << std::endl;
-
+ os << indent << "Min: " << this->GetMinimumOutput()->Get() << std::endl;
+ os << indent << "Max: " << this->GetMaximumOutput()->Get() << std::endl;
+ os << indent << "Mean: " << this->GetMeanOutput()->Get() << std::endl;
+ os << indent << "Covariance: " << this->GetCovarianceOutput()->Get() << std::endl;
+ os << indent << "Correlation: " << this->GetCorrelationOutput()->Get() << std::endl;
}
} // end namespace otb