From e349301a8c313653a36f50b387fcd1b657cbcd6d Mon Sep 17 00:00:00 2001
From: Manuel Grizonnet <manuel.grizonnet@orfeo-toolbox.org>
Date: Mon, 22 Apr 2013 16:09:17 +0200
Subject: [PATCH] BUG: corrections in ReflectanceToSurfaceReflectanceFilter
class
Fistly, the spectral sentivity values tabulate in metadata class must be read by this filter and not by third part application. When no value is available and when there is no user spectral sentivity file provided, default sensitivity filters to 1 are used.
Secondly, the atmospheric coefficients compute by 6S were reordered by the filter using the BandIndexToWavelenghtPosition method. It is wrong in my opinion as it shows on Pleiades image where the red band is the most impacted by the atmosphere as It must be the blue band. I change this and consider that 6S returns atmospheric parameters in the same order as the band in the image file (note that spectral sensitivity tabulates in the metadata class are provided in the same order as bands in the image file.
---
...flectanceToSurfaceReflectanceImageFilter.h | 1 +
...ectanceToSurfaceReflectanceImageFilter.txx | 46 +++++++++++++++----
2 files changed, 39 insertions(+), 8 deletions(-)
diff --git a/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.h b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.h
index 9daf0d175a..782609ab2c 100644
--- a/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.h
+++ b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.h
@@ -231,6 +231,7 @@ protected:
ReflectanceToSurfaceReflectanceImageFilter();
/** Destructor */
virtual ~ReflectanceToSurfaceReflectanceImageFilter() {}
+ void PrintSelf(std::ostream& os, itk::Indent indent) const;
/** Read the aeronet data and extract aerosol optical and water vapor amount. */
//void UpdateAeronetData( const MetaDataDictionaryType dict );
diff --git a/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.txx b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.txx
index 5f81655c1d..a30fc2057e 100644
--- a/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.txx
+++ b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.txx
@@ -99,7 +99,15 @@ ReflectanceToSurfaceReflectanceImageFilter<TInputImage, TOutputImage>
//case where filter function values are not read from an ascii file
else
{
- if (m_FilterFunctionCoef->Size() == 0)
+ try
+ {
+ m_CorrectionParameters->SetWavelengthSpectralBand(imageMetadataInterface->GetSpectralSensitivity());
+ otbMsgDevMacro(<< "use filter available in MetadataInterface " << imageMetadataInterface->GetSpectralSensitivity());
+ }
+
+ catch (itk::ExceptionObject& err)
+ {
+ if (m_FilterFunctionCoef->Size() == 0)
{
otbMsgDevMacro(<< "use dummy filter");
for (unsigned int i = 0; i < this->GetInput()->GetNumberOfComponentsPerPixel(); ++i)
@@ -108,13 +116,13 @@ ReflectanceToSurfaceReflectanceImageFilter<TInputImage, TOutputImage>
}
}
- m_CorrectionParameters->SetWavelengthSpectralBand(m_FilterFunctionCoef);
-
- Parameters2RadiativeTermsPointerType param2Terms = Parameters2RadiativeTermsType::New();
- param2Terms->SetInput(m_CorrectionParameters);
- param2Terms->Update();
- m_AtmosphericRadiativeTerms = param2Terms->GetOutput();
+ m_CorrectionParameters->SetWavelengthSpectralBand(m_FilterFunctionCoef);
+ }
}
+ Parameters2RadiativeTermsPointerType param2Terms = Parameters2RadiativeTermsType::New();
+ param2Terms->SetInput(m_CorrectionParameters);
+ param2Terms->Update();
+ m_AtmosphericRadiativeTerms = param2Terms->GetOutput();
}
template <class TInputImage, class TOutputImage>
@@ -156,7 +164,11 @@ ReflectanceToSurfaceReflectanceImageFilter<TInputImage, TOutputImage>
{
double coef;
double res;
- unsigned int wavelengthPosition = imageMetadataInterface->BandIndexToWavelengthPosition(i);
+ // Don't change the order of atmospheric parameters.
+ //FIXME Need to check in 6S that
+ //unsigned int wavelengthPosition = imageMetadataInterface->BandIndexToWavelengthPosition(i);
+ unsigned int wavelengthPosition = i;
+
coef = static_cast<double>(m_AtmosphericRadiativeTerms->GetTotalGaseousTransmission(wavelengthPosition)
* m_AtmosphericRadiativeTerms->GetDownwardTransmittance(wavelengthPosition)
* m_AtmosphericRadiativeTerms->GetUpwardTransmittance(wavelengthPosition));
@@ -168,6 +180,14 @@ ReflectanceToSurfaceReflectanceImageFilter<TInputImage, TOutputImage>
functor.SetSphericalAlbedo(static_cast<double>(m_AtmosphericRadiativeTerms->GetSphericalAlbedo(wavelengthPosition)));
this->GetFunctorVector().push_back(functor);
+
+ otbMsgDevMacro(<< "Parameters to compute surface reflectance: ");
+
+ otbMsgDevMacro(<< "Band wavelengh position: " << wavelengthPosition);
+ otbMsgDevMacro(<< "Coef (A): " << functor.GetCoefficient());
+ otbMsgDevMacro(<< "Residu: " << functor.GetResidu());
+ otbMsgDevMacro(<< "Spherical albedo: " << functor.GetSphericalAlbedo());
+
}
}
@@ -185,6 +205,16 @@ ReflectanceToSurfaceReflectanceImageFilter<TInputImage, TOutputImage>
this->UpdateFunctors();
}
+
+/* Standard "PrintSelf" method */
+template <class TInputImage, class TOutputImage>
+void
+ReflectanceToSurfaceReflectanceImageFilter<TInputImage, TOutputImage>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+ os << indent << "Correction parameters : " << m_CorrectionParameters << std::endl;
}
+} //end namespace otb
+
#endif
--
GitLab