Skip to content
Snippets Groups Projects
Commit 1dbcce9a authored by Guillaume Pasero's avatar Guillaume Pasero
Browse files

BUG: Mantis-875: avoid division by zero, fix unbiased estimator for ComponentCovariance

parent 052674d9
No related branches found
No related tags found
No related merge requests found
......@@ -452,13 +452,21 @@ PersistentStreamingStatisticsVectorImageFilter<TInputImage, TPrecision>
const RealPixelType& mean = this->GetMeanOutput()->Get();
double regul = 1.0;
double regulComponent = 1.0;
if( m_UseUnbiasedEstimator && nbPixels>1 )
if( m_UseUnbiasedEstimator && nbRelevantPixels>1 )
{
regul =
static_cast< double >( nbRelevantPixels ) /
( static_cast< double >( nbRelevantPixels ) - 1.0 );
}
if( m_UseUnbiasedEstimator && (nbRelevantPixels * numberOfComponent) > 1 )
{
regulComponent =
static_cast< double >(nbRelevantPixels * numberOfComponent) /
( static_cast< double >(nbRelevantPixels * numberOfComponent) - 1.0 );
}
MatrixType cov = cor;
for (unsigned int r = 0; r < numberOfComponent; ++r)
......@@ -473,8 +481,7 @@ PersistentStreamingStatisticsVectorImageFilter<TInputImage, TPrecision>
this->GetComponentMeanOutput()->Set(streamFirstOrderComponentAccumulator / (nbRelevantPixels * numberOfComponent));
this->GetComponentCorrelationOutput()->Set(streamSecondOrderComponentAccumulator / (nbRelevantPixels * numberOfComponent));
this->GetComponentCovarianceOutput()->Set(
(nbRelevantPixels * numberOfComponent) / (nbRelevantPixels * numberOfComponent - 1)
* (this->GetComponentCorrelation()
regulComponent * (this->GetComponentCorrelation()
- (this->GetComponentMean() * this->GetComponentMean())));
}
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment