diff --git a/Modules/Filtering/Statistics/include/otbStreamingStatisticsVectorImageFilter.h b/Modules/Filtering/Statistics/include/otbStreamingStatisticsVectorImageFilter.h index d3da80dceb3d65e08965727eec485e187c8da615..fcb3d0d4f03d0042601510569e8994a449b3c993 100644 --- a/Modules/Filtering/Statistics/include/otbStreamingStatisticsVectorImageFilter.h +++ b/Modules/Filtering/Statistics/include/otbStreamingStatisticsVectorImageFilter.h @@ -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(); diff --git a/Modules/Filtering/Statistics/include/otbStreamingStatisticsVectorImageFilter.txx b/Modules/Filtering/Statistics/include/otbStreamingStatisticsVectorImageFilter.txx index e76b876e578063e4914da29921c93b591d004a35..fbd4e8004e2a6be33dea3a5315e40777b7a176ba 100644 --- a/Modules/Filtering/Statistics/include/otbStreamingStatisticsVectorImageFilter.txx +++ b/Modules/Filtering/Statistics/include/otbStreamingStatisticsVectorImageFilter.txx @@ -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;