diff --git a/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.h b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.h
index f6876b0f2c7dbd56d1762f0d4b4ff322603f7978..30d9c5c4d1e3e2d07645d423e4dbee02aae977ff 100644
--- a/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.h
+++ b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.h
@@ -156,7 +156,7 @@ public:
typedef AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms Parameters2RadiativeTermsType;
typedef Parameters2RadiativeTermsType::Pointer Parameters2RadiativeTermsPointerType;
typedef AtmosphericCorrectionParameters::Pointer CorrectionParametersPointerType;
- typedef AtmosphericRadiativeTerms::Pointer AtmosphericRadiativeTermsPointerType;
+ typedef AtmosphericRadiativeTerms::Pointer AtmosphericRadiativeTermsPointerType;
typedef FilterFunctionValues FilterFunctionValuesType;
@@ -173,9 +173,10 @@ public:
m_IsSetAtmosphericRadiativeTerms = true;
this->Modified();
}
+ itkGetObjectMacro(AtmosphericRadiativeTerms, AtmosphericRadiativeTerms);
/** Get/Set Atmospheric Correction Parameters. */
- itkGetObjectMacro(AtmosphericRadiativeTerms, AtmosphericRadiativeTerms);
+ itkSetObjectMacro(CorrectionParameters, AtmosphericCorrectionParameters);
itkGetObjectMacro(CorrectionParameters, AtmosphericCorrectionParameters);
/** Get/Set Aeronet file name. */
@@ -215,7 +216,6 @@ protected:
/** Initialize the functor vector */
void GenerateOutputInformation();
-
/** Fill AtmosphericRadiativeTerms using image metadata*/
void UpdateAtmosphericRadiativeTerms();
/** Update Functors parameters */
diff --git a/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.txx b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.txx
index 799ea439c06f0d5faf8db363b6091a663343a1d3..004879f09ea7f69dc9c06223a3caf27ee1a5dc35 100644
--- a/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.txx
+++ b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.txx
@@ -76,7 +76,7 @@ ReflectanceToSurfaceReflectanceImageFilter<TInputImage,TOutputImage>
{
m_CorrectionParameters->SetViewingZenithalAngle(90. - imageMetadataInterface->GetSatElevation(dict));
}
-
+
if ((m_CorrectionParameters->GetViewingAzimutalAngle() == 361.))
{
m_CorrectionParameters->SetViewingAzimutalAngle(imageMetadataInterface->GetSatAzimuth(dict));
diff --git a/Code/Radiometry/otbSurfaceAdjencyEffect6SCorrectionSchemeFilter.h b/Code/Radiometry/otbSurfaceAdjencyEffect6SCorrectionSchemeFilter.h
index 8dc7dad340413aae054f8fc1dde2e344a3248746..05e3cdeaf02b9fc978fae5a925b203e97e148077 100644
--- a/Code/Radiometry/otbSurfaceAdjencyEffect6SCorrectionSchemeFilter.h
+++ b/Code/Radiometry/otbSurfaceAdjencyEffect6SCorrectionSchemeFilter.h
@@ -28,6 +28,7 @@
#include "otbUnaryFunctorNeighborhoodImageFilter.h"
#include "itkVariableSizeMatrix.h"
#include "otbAtmosphericRadiativeTerms.h"
+#include "otbAtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms.h"
namespace otb
{
@@ -179,13 +180,22 @@ public:
typedef AtmosphericRadiativeTerms AtmosphericRadiativeTermsType;
typedef typename AtmosphericRadiativeTermsType::Pointer AtmosphericRadiativeTermsPointerType;
+ typedef AtmosphericCorrectionParameters::Pointer CorrectionParametersPointerType;
+ typedef AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms Parameters2RadiativeTermsType;
+ typedef Parameters2RadiativeTermsType::Pointer Parameters2RadiativeTermsPointerType;
+
+ typedef FilterFunctionValues FilterFunctionValuesType;
+ typedef FilterFunctionValuesType::ValuesVectorType CoefVectorType;
+ typedef std::vector<CoefVectorType> FilterFunctionCoefVectorType;
+
+ typedef itk::MetaDataDictionary MetaDataDictionaryType;
/** Storage ponderation values types*/
- typedef itk::VariableSizeMatrix<double> WeightingMatrixType;
- typedef typename std::vector<WeightingMatrixType> WeightingValuesContainerType;
+ typedef itk::VariableSizeMatrix<double> WeightingMatrixType;
+ typedef typename std::vector<WeightingMatrixType> WeightingValuesContainerType;
/** typedef for calculation*/
- typedef typename itk::ConstNeighborhoodIterator<InputImageType> NeighborIterType;
+ typedef typename itk::ConstNeighborhoodIterator<InputImageType> NeighborIterType;
/** Set/Get the Size of the neighbor window. */
void SetWindowRadius(unsigned int rad)
@@ -208,16 +218,40 @@ public:
{
m_AtmosphericRadiativeTerms = atmo;
this->SetNthInput(1, m_AtmosphericRadiativeTerms);
+ m_IsSetAtmosphericRadiativeTerms = true;
this->Modified();
}
- AtmosphericRadiativeTermsPointerType GetAtmosphericRadiativeTerms()
+ itkGetObjectMacro(AtmosphericRadiativeTerms, AtmosphericRadiativeTerms);
+
+ /** Get/Set Atmospheric Correction Parameters. */
+ itkSetObjectMacro(CorrectionParameters, AtmosphericCorrectionParameters);
+ itkGetObjectMacro(CorrectionParameters, AtmosphericCorrectionParameters);
+
+ /** Get/Set Aeronet file name. */
+ itkSetMacro(AeronetFileName, std::string);
+ itkGetMacro(AeronetFileName, std::string);
+
+ /** Get/Set Aeronet file name. */
+ itkSetMacro(FilterFunctionValuesFileName, std::string);
+ itkGetMacro(FilterFunctionValuesFileName, std::string);
+
+ /** Get/Set Filter function coef. */
+ void SetFilterFunctionCoef( FilterFunctionCoefVectorType vect )
+ {
+ m_FilterFunctionCoef = vect;
+ this->Modified();
+ }
+ FilterFunctionCoefVectorType GetFilterFunctionCoef()
{
- return m_AtmosphericRadiativeTerms;
+ return m_FilterFunctionCoef;
}
/** Compute the functor parameters */
void ComputeParameters();
-
+ /** Compute radiative terms if necessary and then updtae functors attibuts. */
+ void GenerateParameters();
+ /** Fill AtmosphericRadiativeTerms using image metadata*/
+ void UpdateAtmosphericRadiativeTerms();
protected:
SurfaceAdjencyEffect6SCorrectionSchemeFilter();
@@ -228,13 +262,12 @@ protected:
virtual void GenerateOutputInformation();
/** Initialize some accumulators before the threads run. */
- virtual void BeforeThreadedGenerateData();
+ //virtual void BeforeThreadedGenerateData();
/** If modified, we need to compute the parameters again */
virtual void Modified();
private:
-
/** Size of the window. */
unsigned int m_WindowRadius;
/** Weighting values for the neighbor pixels.*/
@@ -247,7 +280,19 @@ private:
double m_ZenithalViewingAngle;
/** Radiative terms object */
AtmosphericRadiativeTermsPointerType m_AtmosphericRadiativeTerms;
-
+ bool m_IsSetAtmosphericRadiativeTerms;
+ /** Atmospheric Correction Parameters. */
+ CorrectionParametersPointerType m_CorrectionParameters;
+ /** Path to an Aeronet data file, allows to compute aerosol optical and water vapor amounts. */
+ std::string m_AeronetFileName;
+ /** Path to an filter function values file. */
+ std::string m_FilterFunctionValuesFileName;
+ /** Contains the filter function values (each element is a vector and represnts the values for each channel) */
+ FilterFunctionCoefVectorType m_FilterFunctionCoef;
+ /** Enable/Disable GenerateParameters in GenerateOutputInformation.
+ * Usefull for image view that call GenerateOutputInformation each time you move the full area.
+ */
+ bool m_UseGenerateParameters;
};
} // end namespace otb
diff --git a/Code/Radiometry/otbSurfaceAdjencyEffect6SCorrectionSchemeFilter.txx b/Code/Radiometry/otbSurfaceAdjencyEffect6SCorrectionSchemeFilter.txx
index 3a300a463211f6e7c390bd0bac539f9fbe1c4317..1959c7c913128ce5dfb657ce659ea2ef1fb37fa3 100644
--- a/Code/Radiometry/otbSurfaceAdjencyEffect6SCorrectionSchemeFilter.txx
+++ b/Code/Radiometry/otbSurfaceAdjencyEffect6SCorrectionSchemeFilter.txx
@@ -34,8 +34,9 @@
#include "itkProgressReporter.h"
#include "otbImage.h"
#include "otbSIXSTraits.h"
-#include "otbAtmosphericRadiativeTerms.h"
#include "otbMath.h"
+#include "otbImageMetadataInterfaceFactory.h"
+#include "otbImageMetadataInterfaceBase.h"
namespace otb
{
@@ -46,8 +47,14 @@ SurfaceAdjencyEffect6SCorrectionSchemeFilter<TInputImage, TOutputImage>
{
m_WindowRadius = 1;
m_PixelSpacingInKilometers = 1.;
- m_ZenithalViewingAngle = 0;
+ m_ZenithalViewingAngle = 361.;
m_AtmosphericRadiativeTerms = AtmosphericRadiativeTermsType::New();
+ m_CorrectionParameters = AtmosphericCorrectionParameters::New();
+ m_IsSetAtmosphericRadiativeTerms = false;
+ m_AeronetFileName = "";
+ m_FilterFunctionValuesFileName = "";
+ m_FilterFunctionCoef.clear();
+ m_UseGenerateParameters = true;
}
template <class TInputImage, class TOutputImage>
@@ -62,6 +69,15 @@ SurfaceAdjencyEffect6SCorrectionSchemeFilter<TInputImage, TOutputImage>
if (!inputPtr || !outputPtr)
return;
outputPtr->SetNumberOfComponentsPerPixel(inputPtr->GetNumberOfComponentsPerPixel());
+
+ if(m_UseGenerateParameters)
+ this->GenerateParameters();
+
+ if (!m_ParametersHaveBeenComputed)
+ {
+ this->ComputeParameters();
+ m_ParametersHaveBeenComputed = true;
+ }
}
template <class TInputImage, class TOutputImage>
@@ -73,16 +89,98 @@ SurfaceAdjencyEffect6SCorrectionSchemeFilter<TInputImage, TOutputImage>
m_ParametersHaveBeenComputed = false;
}
+
template <class TInputImage, class TOutputImage>
void
-SurfaceAdjencyEffect6SCorrectionSchemeFilter<TInputImage, TOutputImage>
-::BeforeThreadedGenerateData ()
+SurfaceAdjencyEffect6SCorrectionSchemeFilter<TInputImage,TOutputImage>
+::UpdateAtmosphericRadiativeTerms()
{
- if (!m_ParametersHaveBeenComputed)
- {
- this->ComputeParameters();
- m_ParametersHaveBeenComputed = true;
- }
+ MetaDataDictionaryType dict = this->GetInput()->GetMetaDataDictionary();
+
+ ImageMetadataInterfaceBase::Pointer imageMetadataInterface = ImageMetadataInterfaceFactory::CreateIMI(dict);
+
+ if ((m_CorrectionParameters->GetDay() == 0))
+ {
+ m_CorrectionParameters->SetDay(imageMetadataInterface->GetDay(dict));
+ }
+
+ if ((m_CorrectionParameters->GetMonth() == 0))
+ {
+ m_CorrectionParameters->SetMonth(imageMetadataInterface->GetMonth(dict));
+ }
+
+ if ((m_CorrectionParameters->GetSolarZenithalAngle() == 361.))
+ {
+ m_CorrectionParameters->SetSolarZenithalAngle(90. - imageMetadataInterface->GetSunElevation(dict));
+ }
+
+ if ((m_CorrectionParameters->GetSolarAzimutalAngle() == 361.))
+ {
+ m_CorrectionParameters->SetSolarAzimutalAngle(imageMetadataInterface->GetSunAzimuth(dict));
+ }
+
+ if ((m_CorrectionParameters->GetViewingZenithalAngle() == 361.))
+ {
+ m_CorrectionParameters->SetViewingZenithalAngle(90. - imageMetadataInterface->GetSatElevation(dict));
+ }
+
+ if( m_ZenithalViewingAngle == 361.)
+ {
+ this->SetZenithalViewingAngle(90. - imageMetadataInterface->GetSatElevation(dict));
+ }
+
+ if ((m_CorrectionParameters->GetViewingAzimutalAngle() == 361.))
+ {
+ m_CorrectionParameters->SetViewingAzimutalAngle(imageMetadataInterface->GetSatAzimuth(dict));
+ }
+
+ if(m_AeronetFileName != "")
+ m_CorrectionParameters->UpdateAeronetData( m_AeronetFileName,
+ imageMetadataInterface->GetYear(dict),
+ imageMetadataInterface->GetHour(dict),
+ imageMetadataInterface->GetMinute(dict) );
+
+ // load fiter function values
+ if(m_FilterFunctionValuesFileName != "")
+ {
+ m_CorrectionParameters->LoadFilterFunctionValue( m_FilterFunctionValuesFileName );
+ }
+ // the user has set the filter function values
+ else
+ {
+ if( m_FilterFunctionCoef.size() != this->GetInput()->GetNumberOfComponentsPerPixel() )
+ {
+ itkExceptionMacro(<<"Filter Function and image channels mismatch.");
+ }
+ for(unsigned int i=0; i<this->GetInput()->GetNumberOfComponentsPerPixel(); i++)
+ {
+ FilterFunctionValuesType::Pointer functionValues = FilterFunctionValuesType::New();
+ functionValues->SetFilterFunctionValues(m_FilterFunctionCoef[i]);
+ functionValues->SetMinSpectralValue(imageMetadataInterface->GetFirstWavelengths(dict)[i]);
+ functionValues->SetMaxSpectralValue(imageMetadataInterface->GetLastWavelengths(dict)[i]);
+
+ m_CorrectionParameters->SetWavelenghtSpectralBandWithIndex(i, functionValues);
+ }
+ }
+
+
+ Parameters2RadiativeTermsPointerType param2Terms = Parameters2RadiativeTermsType::New();
+
+ param2Terms->SetInput(m_CorrectionParameters);
+ param2Terms->Update();
+ m_AtmosphericRadiativeTerms = param2Terms->GetOutput();
+}
+
+
+template <class TInputImage, class TOutputImage>
+void
+SurfaceAdjencyEffect6SCorrectionSchemeFilter<TInputImage,TOutputImage>
+::GenerateParameters()
+{
+ if(m_IsSetAtmosphericRadiativeTerms==false)
+ {
+ this->UpdateAtmosphericRadiativeTerms();
+ }
}
template <class TInputImage, class TOutputImage>
diff --git a/Testing/Code/Radiometry/otbSurfaceAdjencyEffect6SCorrectionSchemeFilter.cxx b/Testing/Code/Radiometry/otbSurfaceAdjencyEffect6SCorrectionSchemeFilter.cxx
index dab4d28f3e1f74a0d9c09b13098960d5a5832212..88e74eeca05fe289d3885e5e30fbd27f8990ea13 100644
--- a/Testing/Code/Radiometry/otbSurfaceAdjencyEffect6SCorrectionSchemeFilter.cxx
+++ b/Testing/Code/Radiometry/otbSurfaceAdjencyEffect6SCorrectionSchemeFilter.cxx
@@ -91,6 +91,9 @@ int otbSurfaceAdjencyEffect6SCorrectionSchemeFilter(int argc, char * argv[])
double aerosolOptical(0.);
+ std::vector<double> DVector(1,0);
+ std::vector< std::vector<double> > DVectorVector(1,DVector);
+
std::ifstream fin;
fin.open(paramFile);
//Read input file parameters
@@ -124,6 +127,8 @@ int otbSurfaceAdjencyEffect6SCorrectionSchemeFilter(int argc, char * argv[])
ValuesVectorType vect;
for (unsigned int j=0; j<nbChannel; j++)
{
+ DVector.clear();
+
functionValues = FilterFunctionValuesType::New();
vect.clear();
@@ -141,7 +146,7 @@ int otbSurfaceAdjencyEffect6SCorrectionSchemeFilter(int argc, char * argv[])
while (!fin.eof() && fin.good())
{
fin >> value;
- vect.push_back(value);
+ DVector/*vect*/.push_back(value);
}
fin.close();
@@ -149,17 +154,18 @@ int otbSurfaceAdjencyEffect6SCorrectionSchemeFilter(int argc, char * argv[])
functionValues->SetMinSpectralValue(minSpectralValue);
functionValues->SetMaxSpectralValue(maxSpectralValue);
functionValues->SetUserStep( val );
-
+ DVectorVector.push_back(DVector);
param->SetWavelenghtSpectralBandWithIndex(j, functionValues);
}
corrToRadia->SetInput( param );
corrToRadia->Update();
- filter->SetAtmosphericRadiativeTerms(corrToRadia->GetOutput());
+ //filter->SetAtmosphericRadiativeTerms(corrToRadia->GetOutput());
+ filter->SetFilterFunctionCoef(DVectorVector);
filter->SetWindowRadius(atoi(argv[3]));
filter->SetPixelSpacingInKilometers(static_cast<double>(atof(argv[4])));
- filter->SetZenithalViewingAngle(param->GetViewingZenithalAngle());
+ //filter->SetZenithalViewingAngle(param->GetViewingZenithalAngle());
filter->SetInput(reader->GetOutput());
writer->SetInput(filter->GetOutput());