Commit c60ea4f3 authored by Stéphane Albert's avatar Stéphane Albert

ENH: Use of Scott's formula to compute number of bins of histogram; Cloned...

ENH: Use of Scott's formula to compute number of bins of histogram; Cloned otbStreamingHistogramVectorImageFilter in order to upgrade interface (number of bins per band).
parent 9fad82fa
......@@ -61,8 +61,6 @@ namespace
/*****************************************************************************/
/* STATIC IMPLEMENTATION SECTION */
const int HistogramModel::BINS_OVERSAMPLING_RATE = 3;
/*****************************************************************************/
/* CLASS IMPLEMENTATION SECTION */
......
......@@ -104,8 +104,11 @@ public:
//
// This typedef is used for compatibility with
// itk::Histogram<>::MeasurementType.
itk::NumericTraits< DefaultImageType::InternalPixelType >::RealType
MeasurementType;
itk::NumericTraits< DefaultImageType::InternalPixelType >::RealType RealType;
/**
*/
typedef RealType MeasurementType;
/**
* \class BuildContext
......@@ -209,9 +212,7 @@ protected:
//
// Private constants.
private:
/**
*/
static const int BINS_OVERSAMPLING_RATE;
//
// Private types.
......@@ -228,10 +229,10 @@ private:
private:
/**
*/
#if 0
/*
template< typename TImage >
void template_BuildModel_I();
#endif
*/
/**
*/
......@@ -274,12 +275,13 @@ private slots:
//
// OTB includes (sorted by alphabetic order)
#include "otbStreamingHistogramVectorImageFilter.h"
#include "otbStreamingMinMaxVectorImageFilter.h"
#include "otbStreamingStatisticsVectorImageFilter.h"
//
// Monteverdi includes (sorted by alphabetic order)
#include "Core/mvdAbstractImageModel.h"
#include "Core/mvdStreamingHistogramVectorImageFilter.h"
namespace mvd
{
......@@ -363,7 +365,7 @@ HistogramModel
}
/*******************************************************************************/
#if 0
/*
template< typename TImage >
void
HistogramModel
......@@ -406,13 +408,6 @@ HistogramModel
filterMinMax->Update();
/*
// Extract min/MAX intensities for each band.
// itk::VariableLengthVector< FLOAT_TYPE >
typename MinMaxFilter::PixelType lSrcMin( f );
typename MinMaxFilter::PixelType lSrcMax( filterMinMax->GetMaximum() );
*/
// Extract-convert-remember min/MAX intensities for each band.
m_MinPixel = filterMinMax->GetMinimum();
m_MaxPixel = filterMinMax->GetMaximum();
......@@ -464,7 +459,7 @@ HistogramModel
.arg( QDateTime::currentDateTime().toString( Qt::ISODate ) )
.arg( lMain.elapsed() );
}
#endif
*/
/*******************************************************************************/
template< typename TImageModel >
......@@ -492,6 +487,9 @@ HistogramModel
TImageModel* imageModel = parentImageModel->GetQuicklookModel();
assert( imageModel!=NULL );
CountType components = imageModel->ToImage()->GetNumberOfComponentsPerPixel();
assert( components>0 );
//
// 1st pass: process min/MAX for each band.
......@@ -500,6 +498,7 @@ HistogramModel
lPass1.start();
/*
// Connect min/MAX pipe-section.
typedef
otb::StreamingMinMaxVectorImageFilter<
......@@ -512,18 +511,70 @@ HistogramModel
filterMinMax->GetFilter()->SetNoDataFlag( imageProperties->IsNoDataEnabled() );
filterMinMax->GetFilter()->SetNoDataValue( imageProperties->GetNoData() );
// Connect statistics pipe-section.
filterMinMax->Update();
/*
// Extract min/MAX intensities for each bands.
typename MinMaxFilter::PixelType lSrcMin( filterMinMax->GetMinimum() );
typename MinMaxFilter::PixelType lSrcMax( filterMinMax->GetMaximum() );
*/
// Extract-convert-remember min/MAX intensities for each band.
m_MinPixel = filterMinMax->GetMinimum();
m_MaxPixel = filterMinMax->GetMaximum();
// qDebug() << "min:" << m_MinPixel << "max:" << m_MaxPixel;
// std::cout << "min:" << m_MinPixel << "\tmax:" << m_MaxPixel << std::endl;
*/
// Define histogram-filter type.
typedef
otb::StreamingHistogramVectorImageFilter<
typename TImageModel::SourceImageType >
HistogramFilter;
// Connect statistics pipe-section.
typedef
otb::StreamingStatisticsVectorImageFilter<
typename TImageModel::SourceImageType >
StatisticsFilter;
typename StatisticsFilter::Pointer filterStats( StatisticsFilter::New() );
filterStats->SetInput( imageModel->ToImage() );
filterStats->SetEnableMinMax( true );
filterStats->SetIgnoreUserDefinedValue( imageProperties->IsNoDataEnabled() );
filterStats->SetUserIgnoredValue( imageProperties->GetNoData() );
// Connect statistics pipe-section.
filterStats->Update();
// qDebug() << "min:" << m_MinPixel << "max:" << m_MaxPixel;
// std::cout << "min:" << m_MinPixel << "\tmax:" << m_MaxPixel << std::endl;
// Extract-convert-remember min/MAX intensities for each band.
m_MinPixel = filterStats->GetMinimum();
m_MaxPixel = filterStats->GetMaximum();
// Extract sigmas for each band from covariance matrix.
typename StatisticsFilter::MatrixType covariance(
filterStats->GetFilter()->GetCovariance()
);
typename HistogramFilter::FilterType::CountVectorType bins( components );
typename StatisticsFilter::RealPixelType sums( filterStats->GetSum() );
typename StatisticsFilter::RealPixelType means( filterStats->GetMean() );
for( CountType i=0; i<components; ++i )
{
RealType sigma = sqrt( covariance( i, i ) );
RealType n = sums[ i ] / means[ i ];
// Scott's formula
// See http://en.wikipedia.org/wiki/Histogram#Number_of_bins_and_width
RealType h = 3.5 * sigma / pow( n, 1.0 / 3.0 );
bins[ i ] = ceil( m_MaxPixel[ i ] - m_MinPixel[ i ] ) / h;
}
std::cout << bins;
qDebug() << tr( "%1: Pass #1 - done (%2 ms)." )
.arg( QDateTime::currentDateTime().toString( Qt::ISODate ) )
.arg( lPass1.elapsed() );
......@@ -537,10 +588,6 @@ HistogramModel
lPass2.start();
// Connect histogram-generator pipe-section.
typedef
otb::StreamingHistogramVectorImageFilter<
typename TImageModel::SourceImageType >
HistogramFilter;
typename HistogramFilter::Pointer histogramFilter( HistogramFilter::New() );
......@@ -549,9 +596,7 @@ HistogramModel
// Setup histogram filter.
histogramFilter->GetFilter()->SetHistogramMin( m_MinPixel );
histogramFilter->GetFilter()->SetHistogramMax( m_MaxPixel );
histogramFilter->GetFilter()->SetNumberOfBins(
HistogramModel::BINS_OVERSAMPLING_RATE * 256
);
histogramFilter->GetFilter()->SetNumberOfBins( bins );
histogramFilter->GetFilter()->SetSubSamplingRate( 1 );
histogramFilter->GetFilter()->SetNoDataFlag(
imageProperties->IsNoDataEnabled() );
......
/*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.txt for details.
Some parts of this code are derived from ITK. See ITKCopyright.txt
for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANT2ABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef __otbStreamingHistogramVectorImageFilter_h
#define __otbStreamingHistogramVectorImageFilter_h
#include "otbPersistentImageFilter.h"
#include "otbPersistentFilterStreamingDecorator.h"
#include "otbObjectList.h"
#include "itkStatisticsAlgorithm.h"
#include "itkDenseFrequencyContainer.h"
#include "itkNumericTraits.h"
#include "itkHistogram.h"
namespace otb
{
/** \class PersistentHistogramVectorImageFilter
* \brief Compute the histogram 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.
*
* To reset the temporary data, one should call the Reset() function.
*
* To get the statistics once the regions have been processed via the pipeline, use the Synthetize() method.
*
* \sa PersistentImageFilter
* \ingroup Streamed
* \ingroup Multithreaded
* \ingroup MathematicalStatisticsImageFilters
*
*/
template<class TInputImage>
class ITK_EXPORT PersistentHistogramVectorImageFilter :
public PersistentImageFilter<TInputImage, TInputImage>
{
public:
/** Standard Self typedef */
typedef PersistentHistogramVectorImageFilter Self;
typedef PersistentImageFilter<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(PersistentHistogramVectorImageFilter, 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;
itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension);
/** 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;
/** Types for histogram */
typedef itk::Statistics::DenseFrequencyContainer DFContainerType;
typedef
typename itk::NumericTraits<InternalPixelType>::RealType
HistogramMeasurementRealType;
typedef
itk::Statistics::Histogram<HistogramMeasurementRealType, 1, DFContainerType>
HistogramType;
typedef itk::VariableLengthVector< unsigned int > CountVectorType;
typedef PixelType MeasurementVectorType;
typedef ObjectList<HistogramType> HistogramListType;
typedef typename HistogramListType::Pointer HistogramListPointerType;
typedef typename std::vector<HistogramListPointerType> ArrayHistogramListType;
/** Set the no data value. These value are ignored in histogram
* computation if NoDataFlag is On
*/
itkSetMacro(NoDataValue, InternalPixelType);
/** Set the no data value. These value are ignored in histogram
* computation if NoDataFlag is On
*/
itkGetConstReferenceMacro(NoDataValue, InternalPixelType);
/** Set the NoDataFlag. If set to true, samples with values equal to
* m_NoDataValue are ignored.
*/
itkSetMacro(NoDataFlag, bool);
/** Get the NoDataFlag. If set to true, samples with values equal to
* m_NoDataValue are ignored.
*/
itkGetMacro(NoDataFlag, bool);
/** Toggle the NoDataFlag. If set to true, samples with values equal to
* m_NoDataValue are ignored.
*/
itkBooleanMacro(NoDataFlag);
inline void SetNumberOfBins( unsigned int i, CountVectorType::ValueType size )
{
m_Size[ i ] = size;
}
inline void SetNumberOfBins( const CountVectorType& size )
{
m_Size = size;
}
/** Return the computed histogram list */
HistogramListType* GetHistogramListOutput();
const HistogramListType* GetHistogramListOutput() const;
/** Get the minimum values for histograms */
itkGetConstReferenceMacro(HistogramMin,MeasurementVectorType);
/** Set the minimum values for histograms */
itkSetMacro(HistogramMin,MeasurementVectorType);
/** Get the maximum values for histograms */
itkGetConstReferenceMacro(HistogramMax,MeasurementVectorType);
/** Set the maximum values for histograms */
itkSetMacro(HistogramMax,MeasurementVectorType);
/** Set the subsampling rate */
itkSetMacro(SubSamplingRate, unsigned int);
/** Get the subsampling rate */
itkGetMacro(SubSamplingRate, unsigned int);
/** Make a DataObject of the correct type to be used as the specified
* output.
*/
virtual DataObjectPointer MakeOutput(unsigned int idx);
/** 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:
PersistentHistogramVectorImageFilter();
virtual ~PersistentHistogramVectorImageFilter() {}
virtual void PrintSelf(std::ostream& os, itk::Indent indent) const;
/** Multi-thread version GenerateData. */
void ThreadedGenerateData(const RegionType& outputRegionForThread, int threadId);
private:
PersistentHistogramVectorImageFilter(const Self &); //purposely not implemented
void operator =(const Self&); //purposely not implemented
ArrayHistogramListType m_ThreadHistogramList;
CountVectorType m_Size;
MeasurementVectorType m_HistogramMin;
MeasurementVectorType m_HistogramMax;
bool m_NoDataFlag;
InternalPixelType m_NoDataValue;
/** Set the subsampling along each direction */
unsigned int m_SubSamplingRate;
}; // end of class PersistentStatisticsVectorImageFilter
/**===========================================================================*/
/** \class StreamingHistogramVectorImageFilter
* \brief This class streams the whole input image through the PersistentHistogramVectorImageFilter.
*
* This way, it allows to compute the min/max of this image. It calls the
* Reset() method of the PersistentHistogramVectorImageFilter before streaming the image and the
* Synthetize() method of the PersistentHistogramVectorImageFilter after having streamed the image
* to compute the statistics. The accessor on the results are wrapping the accessors of the
* internal PersistentMinMaxImageFilter.
*
* \sa PersistentStatisticsVectorImageFilter
* \sa PersistentImageFilter
* \sa PersistentFilterStreamingDecorator
* \sa StreamingImageVirtualWriter
* \ingroup Streamed
* \ingroup Multithreaded
* \ingroup MathematicalStatisticsImageFilters
*/
template<class TInputImage>
class ITK_EXPORT StreamingHistogramVectorImageFilter :
public PersistentFilterStreamingDecorator<PersistentHistogramVectorImageFilter<TInputImage> >
{
public:
/** Standard Self typedef */
typedef StreamingHistogramVectorImageFilter Self;
typedef PersistentFilterStreamingDecorator
<PersistentHistogramVectorImageFilter<TInputImage> > Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Type macro */
itkNewMacro(Self);
/** Creation through object factory macro */
itkTypeMacro(StreamingHistogramVectorImageFilter, PersistentFilterStreamingDecorator);
typedef TInputImage InputImageType;
typedef typename Superclass::FilterType InternalFilterType;
/** Types needed for histograms */
typedef typename InternalFilterType::HistogramType HistogramType;
typedef typename InternalFilterType::HistogramListType HistogramListType;
void SetInput(InputImageType * input)
{
this->GetFilter()->SetInput(input);
}
const InputImageType * GetInput()
{
return this->GetFilter()->GetInput();
}
/** Return the computed histogram */
HistogramListType* GetHistogramList()
{
return this->GetFilter()->GetHistogramListOutput();
}
protected:
/** Constructor */
StreamingHistogramVectorImageFilter() {};
/** Destructor */
virtual ~StreamingHistogramVectorImageFilter() {}
private:
StreamingHistogramVectorImageFilter(const Self &); //purposely not implemented
void operator =(const Self&); //purposely not implemented
};
} // end namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
#include "mvdStreamingHistogramVectorImageFilter.txx"
#endif
#endif
/*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See OTBCopyright.txt for details.
Some parts of this code are derived from ITK. See ITKCopyright.txt
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 __otbStreamingHistogramVectorImageFilter_txx
#define __otbStreamingHistogramVectorImageFilter_txx
#include "Core/mvdStreamingHistogramVectorImageFilter.h"
#include "itkImageRegionIterator.h"
#include "itkImageRegionConstIteratorWithIndex.h"
#include "itkNumericTraits.h"
#include "itkProgressReporter.h"
#include "otbMacro.h"
namespace otb
{
template<class TInputImage>
PersistentHistogramVectorImageFilter<TInputImage>
::PersistentHistogramVectorImageFilter() :
m_ThreadHistogramList(),
m_Size(),
m_HistogramMin(),
m_HistogramMax(),
m_NoDataFlag(false),
m_NoDataValue(itk::NumericTraits<InternalPixelType>::Zero),
m_SubSamplingRate(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 and histogram list
m_Size.Fill(255);
HistogramListPointerType output = static_cast<HistogramListType*>(this->MakeOutput(1).GetPointer());
this->itk::ProcessObject::SetNthOutput(1, output.GetPointer());
}
template<class TInputImage>
itk::DataObject::Pointer
PersistentHistogramVectorImageFilter<TInputImage>
::MakeOutput(unsigned int output)
{
itk::DataObject::Pointer ret;
switch (output)
{
case 0:
ret = static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
break;
case 1:
ret = static_cast<itk::DataObject*>(HistogramListType::New().GetPointer());
break;
}
return ret;
}
template<class TInputImage>
typename PersistentHistogramVectorImageFilter<TInputImage>::HistogramListType*
PersistentHistogramVectorImageFilter<TInputImage>
::GetHistogramListOutput()
{
return static_cast<HistogramListType*>(this->itk::ProcessObject::GetOutput(1));
}
template<class TInputImage>
const typename PersistentHistogramVectorImageFilter<TInputImage>::HistogramListType*
PersistentHistogramVectorImageFilter<TInputImage>
::GetHistogramListOutput() const
{
return static_cast<const HistogramListType*>(this->itk::ProcessObject::GetOutput(1));
}
template<class TInputImage>
void
PersistentHistogramVectorImageFilter<TInputImage>
::GenerateOutputInformation()
{
Superclass::GenerateOutputInformation();
if (this->GetInput())
{
this->GetOutput()->CopyInformation(this->GetInput());
this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
{
this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
}
}
}
template<class TInputImage>
void
PersistentHistogramVectorImageFilter<TInputImage>
::AllocateOutputs()
{
// This is commented to prevent the streaming of the whole image for the first stream strip
// It shall not cause any problem because the output image of this filter is not intended to be used.
//InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
//this->GraftOutput( image );
// Nothing that needs to be allocated for the remaining outputs
}
template<class TInputImage>
void
PersistentHistogramVectorImageFilter<TInputImage>
::Reset()
{
TInputImage * inputPtr = const_cast<TInputImage *>(this->GetInput());
inputPtr->UpdateOutputInformation();
unsigned int numberOfThreads = this->GetNumberOfThreads();
unsigned int numberOfComponent = inputPtr->GetNumberOfComponentsPerPixel();
bool clipBins = true;
// if histogram Min and Max have the wrong size : set to default [0, 255]
if (m_HistogramMin.Size() != numberOfComponent ||
m_HistogramMax.Size() != numberOfComponent)
{
m_HistogramMin.SetSize(numberOfComponent);
m_HistogramMax.SetSize(numberOfComponent);
m_HistogramMin.Fill(itk::NumericTraits<InternalPixelType>::Zero);
m_HistogramMax.Fill(255);
}
// Setup output histogram
HistogramListType* outputHisto = this->GetHistogramListOutput();
outputHisto->Clear();
for (unsigned int k=0; k<numberOfComponent; ++k)
{
typename HistogramType::MeasurementVectorType bandMin, bandMax;
bandMin[0] = m_HistogramMin[k];
bandMax[0] = m_HistogramMax[k];
typename HistogramType::Pointer histogram = HistogramType::New();
histogram->SetClipBinsAtEnds(clipBins);
typename HistogramType::SizeType size;
size.Fill( m_Size[ k ] );
histogram->Initialize( size, bandMin, bandMax );
outputHisto->PushBack(histogram);
}
// Setup HistogramLists for each thread
m_ThreadHistogramList.clear();
for (unsigned int i=0; i<numberOfThreads; ++i)
{
HistogramListPointerType histoList = HistogramListType::New();
histoList->Clear();
for (unsigned int k=0; k<numberOfComponent; ++k)
{
typename HistogramType::MeasurementVectorType bandMin, bandMax;
bandMin[0] = m_HistogramMin[k];
bandMax[0] = m_HistogramMax[k];
typename HistogramType::Pointer histogram = HistogramType::New();
histogram->SetClipBinsAtEnds(clipBins);
typename HistogramType::SizeType size;
size.Fill( m_Size[ k ] );
histogram->Initialize(size, bandMin, bandMax );
histoList->PushBack(histogram);
}
m_ThreadHistogramList.push_back(histoList);
}
}
template<class TInputImage>
void
PersistentHistogramVectorImageFilter<TInputImage>
::Synthetize()
{
HistogramListType* outputHisto = this->GetHistogramListOutput();
int numberOfThreads = this->GetNumberOfThreads();
unsigned int numberOfComponent = this->GetInput()->GetNumberOfComponentsPerPixel();
// copy histograms to output
for (int i = 0; i < numberOfThreads; ++i)
{
for (unsigned int j = 0; j < numberOfComponent; ++j)
{
HistogramType* outHisto = outputHisto->GetNthElement(j);
HistogramType* threadHisto = m_ThreadHistogramList[i]->GetNthElement(j);
typename HistogramType::Iterator iterOutput = outHisto->Begin();
typename HistogramType::Iterator iterThread = threadHisto->Begin();
while (iterOutput != outHisto->End() && iterThread != threadHisto->End())
{
iterOutput.SetFrequency(iterOutput.GetFrequency()+iterThread.GetFrequency());
++iterOutput;
++iterThread;
}
}
}
}
template<class TInputImage>
void
PersistentHistogramVectorImageFilter<TInputImage>
::ThreadedGenerateData(const RegionType& outputRegionForThread, int threadId)
{
/**
* Grab the input
*/
InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput());
// support progress methods/callbacks
itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
typename HistogramType::IndexType index;
itk::ImageRegionConstIteratorWithIndex<TInputImage> it(inputPtr, outputRegionForThread);
it.GoToBegin();
bool skipSample = false;
// do the work
while (!it.IsAtEnd())
{
if (m_SubSamplingRate > 1)
{
skipSample = false;
for (unsigned int i=0; i<InputImageDimension; ++i)
{