An error occurred while loading the file. Please try again.
-
coupesrases authored5ee03564
#ifndef __MultitemporalStatisticsImageFilter_hxx#define __MultitemporalStatisticsImageFilter_hxx#include "otbMultitemporalStatisticsImageFilter.h"#include "itkProgressReporter.h"namespace otb{/** * */template <class TImage, class TFunctor>PersitentMultitemporalStatisticsImageFilter<TImage, TFunctor>::PersitentMultitemporalStatisticsImageFilter() { // 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 matrices for (int i = 3; i < 7; ++i) { typename RealMatrixObjectType::Pointer output = static_cast<RealMatrixObjectType*>(this->MakeOutput(i).GetPointer()); this->itk::ProcessObject::SetNthOutput(i, output.GetPointer()); } InputImagePointer inputPtr = const_cast<TImage *>(this->GetInput(0)); unsigned int n = 1; RealMatrixType nullMatrix(n,n,itk::NumericTraits<RealType>::Zero); this->GetMeanOutput()->Set(nullMatrix); this->GetSigmaOutput()->Set(nullMatrix); }template <class TImage, class TFunctor>typename itk::DataObject::PointerPersitentMultitemporalStatisticsImageFilter<TImage, TFunctor>::MakeOutput(DataObjectPointerArraySizeType output){ switch (output) { case 0: return static_cast<itk::DataObject*>(TImage::New().GetPointer()); break; case 1: return static_cast<itk::DataObject*>(PixelObjectType::New().GetPointer()); break; case 2: return static_cast<itk::DataObject*>(PixelObjectType::New().GetPointer()); break; case 3: return static_cast<itk::DataObject*>(RealMatrixObjectType::New().GetPointer()); break; case 4: return static_cast<itk::DataObject*>(RealMatrixObjectType::New().GetPointer()); break; case 5: case 6:7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
return static_cast<itk::DataObject*>(RealObjectType::New().GetPointer());
break;
default:
// might as well make an image
return static_cast<itk::DataObject*>(TImage::New().GetPointer());
break;
}
}
template <class TImage, class TFunctor>
typename PersitentMultitemporalStatisticsImageFilter<TImage, TFunctor>::RealMatrixObjectType*
PersitentMultitemporalStatisticsImageFilter<TImage, TFunctor>
::GetMeanOutput()
{
return static_cast<RealMatrixObjectType*>(this->itk::ProcessObject::GetOutput(3));
}
template <class TImage, class TFunctor>
const typename PersitentMultitemporalStatisticsImageFilter<TImage, TFunctor>::RealMatrixObjectType*
PersitentMultitemporalStatisticsImageFilter<TImage, TFunctor>
::GetMeanOutput() const
{
return static_cast<const RealMatrixObjectType*>(this->itk::ProcessObject::GetOutput(3));
}
template <class TImage, class TFunctor>
typename PersitentMultitemporalStatisticsImageFilter<TImage, TFunctor>::RealMatrixObjectType*
PersitentMultitemporalStatisticsImageFilter<TImage, TFunctor>
::GetSigmaOutput()
{
return static_cast<RealMatrixObjectType*>(this->itk::ProcessObject::GetOutput(4));
}
template <class TImage, class TFunctor>
const typename PersitentMultitemporalStatisticsImageFilter<TImage, TFunctor>::RealMatrixObjectType*
PersitentMultitemporalStatisticsImageFilter<TImage, TFunctor>
::GetSigmaOutput() const
{
return static_cast<const RealMatrixObjectType*>(this->itk::ProcessObject::GetOutput(4));
}
template <class TImage, class TFunctor>
void
PersitentMultitemporalStatisticsImageFilter<TImage, TFunctor>
::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 TImage, class TFunctor>
void
PersitentMultitemporalStatisticsImageFilter<TImage, TFunctor>
::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
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
}
/*
* Synthetize internal counters
*/
template <class TImage, class TFunctor>
void
PersitentMultitemporalStatisticsImageFilter<TImage, TFunctor>
::Synthetize()
{
// Prepare container for final results
InputImagePointer inputPtr = const_cast<TImage *>(this->GetInput(0));
unsigned int n = inputPtr->GetNumberOfComponentsPerPixel();
m_FinalStatsContainer = ThreadResultsContainer(n);
// Gather threads containers
int numberOfThreads = this->GetNumberOfThreads();
for (unsigned int i = 0; i < numberOfThreads; ++i)
{
m_FinalStatsContainer.Update(m_StatsContainers[i]);
}
// Create final stats matrices
RealMatrixType means(n,n,itk::NumericTraits<RealType>::Zero);
RealMatrixType stds(n,n,itk::NumericTraits<RealType>::Zero);
// Update container for final results
for (unsigned int i = 0 ; i < n ; i++)
{
for (unsigned int j = 0 ; j < i ; j++)
{
long count = m_FinalStatsContainer.m_count[i][j];
RealType sum = m_FinalStatsContainer.m_sum[i][j];
RealType sumOfSquares = m_FinalStatsContainer.m_sqSum[i][j];
RealType mean = itk::NumericTraits<RealType>::Zero;
RealType sigma = itk::NumericTraits<RealType>::Zero;
if (count > 0)
{
// compute statistics
mean = sum / static_cast<RealType>(count);
if (count > 1)
{
// unbiased estimate
RealType variance = (sumOfSquares - (sum * sum / static_cast<RealType>(count)))
/ static_cast<RealType>(count - 1);
sigma = vcl_sqrt(variance);
}
}
means[i][j] = mean;
stds[i][j] = sigma;
} // next j
} // next i
// Set the outputs
this->GetMeanOutput()->Set(means);
this->GetSigmaOutput()->Set(stds);
}
/*
* Reset internal counters
*/
template <class TImage, class TFunctor>
void