Commit 6daa2d60 authored by Ludovic Hussonnois's avatar Ludovic Hussonnois

Merge branch 'bugfix-jira-1226' into develop

parents 6950330c dfb05aa0
......@@ -86,44 +86,39 @@ private:
// Samples
typedef double ValueType;
typedef itk::VariableLengthVector<ValueType> MeasurementType;
typedef itk::VariableSizeMatrix<ValueType> MatrixValueType;
unsigned int nbSamples = 0;
unsigned int nbBands = 0;
FloatVectorImageListType* imageList = GetParameterImageList("il");
FloatVectorImageListType::InternalContainerSizeType nbImages = imageList->Size();
// Initialization, all image have same size and number of band/component
FloatVectorImageType* firstImage = imageList->GetNthElement(0);
nbBands = firstImage->GetNumberOfComponentsPerPixel();
// Build a Measurement Vector of mean
MeasurementType mean;
MatrixValueType mean(nbBands, static_cast<unsigned int>(nbImages));
mean.Fill(itk::NumericTraits<MatrixValueType::ValueType>::Zero);
// Build a MeasurementVector of variance
MeasurementType variance;
// Build a Measurement Matrix of variance
MatrixValueType variance(nbBands, static_cast<unsigned int>(nbImages));
variance.Fill(itk::NumericTraits<MatrixValueType::ValueType>::Zero);
FloatVectorImageListType* imageList = GetParameterImageList("il");
// Build a Measurement Matrix of nbSamples
MatrixValueType nbSamples(nbBands, static_cast<unsigned int>(nbImages));
nbSamples.Fill(itk::NumericTraits<MatrixValueType::ValueType>::Zero);
//Iterate over all input images
for (unsigned int imageId = 0; imageId < imageList->Size(); ++imageId)
for (unsigned int imageId = 0; imageId < nbImages; ++imageId)
{
FloatVectorImageType* image = imageList->GetNthElement(imageId);
if (nbBands == 0)
{
nbBands = image->GetNumberOfComponentsPerPixel();
}
else if (nbBands != image->GetNumberOfComponentsPerPixel())
if (nbBands != image->GetNumberOfComponentsPerPixel())
{
itkExceptionMacro(<< "The image #" << imageId + 1 << " has " << image->GetNumberOfComponentsPerPixel()
<< " bands, while the image #1 has " << nbBands );
}
FloatVectorImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();
//Set the measurement vectors size if it's the first iteration
if (imageId == 0)
{
mean.SetSize(nbBands);
mean.Fill(0.);
variance.SetSize(nbBands);
variance.Fill(0.);
}
// Compute Statistics of each VectorImage
StreamingStatisticsVImageFilterType::Pointer statsEstimator = StreamingStatisticsVImageFilterType::New();
std::ostringstream processName;
......@@ -131,41 +126,88 @@ private:
AddProcess(statsEstimator->GetStreamer(), processName.str().c_str());
statsEstimator->SetInput(image);
statsEstimator->GetStreamer()->SetAutomaticAdaptativeStreaming(GetParameterInt("ram"));
if( HasValue( "bv" )==true )
if( HasValue( "bv" ) )
{
statsEstimator->SetIgnoreUserDefinedValue(true);
statsEstimator->SetUserIgnoredValue(GetParameterFloat("bv"));
}
statsEstimator->Update();
mean += statsEstimator->GetMean();
for (unsigned int itBand = 0; itBand < nbBands; itBand++)
MeasurementType nbRelevantPixels = statsEstimator->GetNbRelevantPixels();
MeasurementType meanPerBand = statsEstimator->GetMean();
for(unsigned int itBand = 0; itBand < nbBands; itBand++)
{
mean(itBand, imageId) = meanPerBand[itBand];
variance(itBand, imageId) = (statsEstimator->GetCovariance())( itBand, itBand );
nbSamples(itBand, imageId) = nbRelevantPixels[itBand];
}
}
// Compute total mean and pooled variation for each band of the image list
MeasurementType totalSamplesPerBand;
totalSamplesPerBand.SetSize(nbBands);
totalSamplesPerBand.Fill(itk::NumericTraits<MeasurementType::ValueType>::Zero);
MeasurementType totalMeanPerBand;
totalMeanPerBand.SetSize(nbBands);
totalMeanPerBand.Fill(itk::NumericTraits<MeasurementType::ValueType>::Zero);
MeasurementType totalVariancePerBand;
totalVariancePerBand.SetSize(nbBands);
totalVariancePerBand.Fill(itk::NumericTraits<MeasurementType::ValueType>::Zero);
for (unsigned int imageId = 0; imageId < nbImages; ++imageId)
{
for(unsigned int itBand = 0; itBand < nbBands; itBand++)
{
variance[itBand] += (size[0] * size[1] - 1) * (statsEstimator->GetCovariance())(itBand, itBand);
MeasurementType::ValueType nbSample = nbSamples(itBand, imageId);
totalSamplesPerBand[itBand] += nbSample;
totalMeanPerBand[itBand] += mean(itBand, imageId) * nbSample;
totalVariancePerBand[itBand] += variance(itBand, imageId) * (nbSample - 1);
}
//Increment nbSamples
nbSamples += size[0] * size[1];
}
//Divide by the number of input images to get the mean over all layers
mean /= imageList->Size();
//Apply the pooled variance formula
variance /= (nbSamples - imageList->Size());
// Check 0 division
for(unsigned int itBand = 0; itBand < nbBands; itBand++)
{
MeasurementType::ValueType nbSample = totalSamplesPerBand[itBand];
if ( nbSample > nbImages )
{
totalVariancePerBand[itBand] /= (nbSample - nbImages);
}
else
{
totalVariancePerBand[itBand] = itk::NumericTraits<ValueType>::Zero;
}
if ( nbSample != 0 )
{
totalMeanPerBand[itBand] /= nbSample;
}
else
{
totalMeanPerBand[itBand] = itk::NumericTraits<ValueType>::Zero;
}
}
MeasurementType stddev;
stddev.SetSize(nbBands);
stddev.Fill(0.);
for (unsigned int i = 0; i < variance.GetSize(); ++i)
stddev.Fill(itk::NumericTraits<MeasurementType::ValueType>::Zero);
for (unsigned int i = 0; i < totalVariancePerBand.GetSize(); ++i)
{
stddev[i] = vcl_sqrt(variance[i]);
stddev[i] = vcl_sqrt(totalVariancePerBand[i]);
}
if( HasValue( "out" )==true )
if( HasValue( "out" ) )
{
// Write the Statistics via the statistic writer
typedef otb::StatisticsXMLFileWriter<MeasurementType> StatisticsWriter;
StatisticsWriter::Pointer writer = StatisticsWriter::New();
writer->SetFileName(GetParameterString("out"));
writer->AddInput("mean", mean);
writer->AddInput("mean", totalMeanPerBand);
writer->AddInput("stddev", stddev);
writer->Update();
}
......
......@@ -88,6 +88,7 @@ public:
/** Type to use for computations. */
typedef itk::VariableSizeMatrix<PrecisionType> MatrixType;
typedef itk::VariableLengthVector<PrecisionType> RealPixelType;
typedef itk::VariableLengthVector<unsigned long> CountType;
/** Type of DataObjects used for outputs */
typedef itk::SimpleDataObjectDecorator<RealType> RealObjectType;
......@@ -95,6 +96,15 @@ public:
typedef itk::SimpleDataObjectDecorator<PixelType> PixelObjectType;
typedef itk::SimpleDataObjectDecorator<RealPixelType> RealPixelObjectType;
typedef itk::SimpleDataObjectDecorator<MatrixType> MatrixObjectType;
typedef itk::SimpleDataObjectDecorator<CountType> CountObjectType;
/** Return the number of relevant pixels **/
CountType GetNbRelevantPixels() const
{
return this->GetNbRelevantPixelsOutput()->Get();
}
CountObjectType* GetNbRelevantPixelsOutput();
const CountObjectType* GetNbRelevantPixelsOutput() const;
/** Return the computed Min */
PixelType GetMinimum() const
......@@ -300,6 +310,8 @@ public:
typedef typename StatFilterType::RealPixelObjectType RealPixelObjectType;
typedef typename StatFilterType::MatrixType MatrixType;
typedef typename StatFilterType::MatrixObjectType MatrixObjectType;
typedef typename StatFilterType::CountType CountType;
typedef typename StatFilterType::CountObjectType CountObjectType;
typedef typename StatFilterType::InternalPixelType InternalPixelType;
......@@ -313,6 +325,21 @@ public:
return this->GetFilter()->GetInput();
}
/** Return the number of relevant pixels **/
CountType GetNbRelevantPixels() const
{
return this->GetFilter()->GetNbRelevantPixelsOutput()->Get();
}
CountObjectType* GetNbRelevantPixelsOutput()
{
return this->GetFilter()->GetNbRelevantPixelsOutput();
}
const CountObjectType* GetNbRelevantPixelsOutput() const
{
return this->GetFilter()->GetNbRelevantPixelsOutput();
}
/** Return the computed Minimum. */
RealPixelType GetMinimum() const
{
......@@ -383,7 +410,7 @@ public:
return this->GetFilter()->GetCovarianceOutput();
}
/** Return the computed Covariance. */
/** Return the computed Correlation. */
MatrixType GetCorrelation() const
{
return this->GetFilter()->GetCorrelationOutput()->Get();
......@@ -425,7 +452,7 @@ public:
return this->GetFilter()->GetComponentCovarianceOutput();
}
/** Return the computed Covariance. */
/** Return the computed Correlation. */
RealType GetComponentCorrelation() const
{
return this->GetFilter()->GetComponentCorrelationOutput()->Get();
......
......@@ -43,7 +43,7 @@ PersistentStreamingStatisticsVectorImageFilter<TInputImage, TPrecision>
// allocate the data objects for the outputs which are
// just decorators around vector/matrix types
for (unsigned int i = 1; i < 10; ++i)
for (unsigned int i = 1; i < 11; ++i)
{
this->itk::ProcessObject::SetNthOutput(i, this->MakeOutput(i).GetPointer());
}
......@@ -83,6 +83,9 @@ PersistentStreamingStatisticsVectorImageFilter<TInputImage, TPrecision>
// component mean, component covariance, component correlation
return static_cast<itk::DataObject*>(RealObjectType::New().GetPointer());
break;
case 10:
// relevant pixel
return static_cast<itk::DataObject*>(CountObjectType::New().GetPointer());
default:
// might as well make an image
return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
......@@ -90,6 +93,23 @@ PersistentStreamingStatisticsVectorImageFilter<TInputImage, TPrecision>
}
}
template<class TInputImage, class TPrecision>
typename PersistentStreamingStatisticsVectorImageFilter<TInputImage, TPrecision>::CountObjectType*
PersistentStreamingStatisticsVectorImageFilter<TInputImage, TPrecision>
::GetNbRelevantPixelsOutput()
{
return static_cast<CountObjectType*>(this->itk::ProcessObject::GetOutput(10));
}
template<class TInputImage, class TPrecision>
const typename PersistentStreamingStatisticsVectorImageFilter<TInputImage, TPrecision>::CountObjectType*
PersistentStreamingStatisticsVectorImageFilter<TInputImage, TPrecision>
::GetNbRelevantPixelsOutput() const
{
return static_cast<const CountObjectType*>(this->itk::ProcessObject::GetOutput(10));
}
template<class TInputImage, class TPrecision>
typename PersistentStreamingStatisticsVectorImageFilter<TInputImage, TPrecision>::PixelObjectType*
PersistentStreamingStatisticsVectorImageFilter<TInputImage, TPrecision>
......@@ -419,10 +439,15 @@ PersistentStreamingStatisticsVectorImageFilter<TInputImage, TPrecision>
);
}
unsigned int nbRelevantPixels =
unsigned int nbRelevantPixel =
nbPixels - (ignoredInfinitePixelCount + ignoredUserPixelCount);
if( nbRelevantPixels==0 )
CountType nbRelevantPixels(numberOfComponent);
nbRelevantPixels.Fill(nbRelevantPixel);
this->GetNbRelevantPixelsOutput()->Set(nbRelevantPixels);
if( nbRelevantPixel==0 )
{
itkExceptionMacro(
"Statistics cannot be calculated with zero relevant pixels."
......@@ -438,15 +463,15 @@ PersistentStreamingStatisticsVectorImageFilter<TInputImage, TPrecision>
if (m_EnableFirstOrderStats)
{
this->GetComponentMeanOutput()->Set(streamFirstOrderComponentAccumulator / (nbRelevantPixels * numberOfComponent));
this->GetComponentMeanOutput()->Set(streamFirstOrderComponentAccumulator / (nbRelevantPixel * numberOfComponent));
this->GetMeanOutput()->Set(streamFirstOrderAccumulator / nbRelevantPixels);
this->GetMeanOutput()->Set(streamFirstOrderAccumulator / nbRelevantPixel);
this->GetSumOutput()->Set(streamFirstOrderAccumulator);
}
if (m_EnableSecondOrderStats)
{
MatrixType cor = streamSecondOrderAccumulator / nbRelevantPixels;
MatrixType cor = streamSecondOrderAccumulator / nbRelevantPixel;
this->GetCorrelationOutput()->Set(cor);
const RealPixelType& mean = this->GetMeanOutput()->Get();
......@@ -454,18 +479,18 @@ PersistentStreamingStatisticsVectorImageFilter<TInputImage, TPrecision>
double regul = 1.0;
double regulComponent = 1.0;
if( m_UseUnbiasedEstimator && nbRelevantPixels>1 )
if( m_UseUnbiasedEstimator && nbRelevantPixel>1 )
{
regul =
static_cast< double >( nbRelevantPixels ) /
( static_cast< double >( nbRelevantPixels ) - 1.0 );
static_cast< double >( nbRelevantPixel ) /
( static_cast< double >( nbRelevantPixel ) - 1.0 );
}
if( m_UseUnbiasedEstimator && (nbRelevantPixels * numberOfComponent) > 1 )
if( m_UseUnbiasedEstimator && (nbRelevantPixel * numberOfComponent) > 1 )
{
regulComponent =
static_cast< double >(nbRelevantPixels * numberOfComponent) /
( static_cast< double >(nbRelevantPixels * numberOfComponent) - 1.0 );
static_cast< double >(nbRelevantPixel * numberOfComponent) /
( static_cast< double >(nbRelevantPixel * numberOfComponent) - 1.0 );
}
MatrixType cov = cor;
......@@ -478,8 +503,8 @@ PersistentStreamingStatisticsVectorImageFilter<TInputImage, TPrecision>
}
this->GetCovarianceOutput()->Set(cov);
this->GetComponentMeanOutput()->Set(streamFirstOrderComponentAccumulator / (nbRelevantPixels * numberOfComponent));
this->GetComponentCorrelationOutput()->Set(streamSecondOrderComponentAccumulator / (nbRelevantPixels * numberOfComponent));
this->GetComponentMeanOutput()->Set(streamFirstOrderComponentAccumulator / (nbRelevantPixel * numberOfComponent));
this->GetComponentCorrelationOutput()->Set(streamSecondOrderComponentAccumulator / (nbRelevantPixel * numberOfComponent));
this->GetComponentCovarianceOutput()->Set(
regulComponent * (this->GetComponentCorrelation()
- (this->GetComponentMean() * this->GetComponentMean())));
......@@ -586,6 +611,7 @@ PersistentStreamingStatisticsVectorImageFilter<TImage, TPrecision>
os << indent << "Mean: " << this->GetMeanOutput()->Get() << std::endl;
os << indent << "Covariance: " << this->GetCovarianceOutput()->Get() << std::endl;
os << indent << "Correlation: " << this->GetCorrelationOutput()->Get() << std::endl;
os << indent << "Relevant pixel: " << this->GetNbRelevantPixelsOutput()->Get() << std::endl;
os << indent << "Component Mean: " << this->GetComponentMeanOutput()->Get() << std::endl;
os << indent << "Component Covariance: " << this->GetComponentCovarianceOutput()->Get() << std::endl;
os << indent << "Component Correlation: " << this->GetComponentCorrelationOutput()->Get() << std::endl;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment