From c254bca3d97ef9d0471691c395d789051bd78eb3 Mon Sep 17 00:00:00 2001
From: Cyrille Valladeau <cyrille.valladeau@c-s.fr>
Date: Mon, 24 Aug 2009 17:55:14 +0200
Subject: [PATCH] ENH : radio correctio in progress
---
Code/IO/otbDefaultImageMetadataInterface.h | 13 ++
Code/IO/otbIkonosImageMetadataInterface.cxx | 89 ++++++++++
Code/IO/otbIkonosImageMetadataInterface.h | 6 +
Code/IO/otbImageMetadataInterfaceBase.h | 8 +
.../IO/otbQuickBirdImageMetadataInterface.cxx | 95 +++++++++++
Code/IO/otbQuickBirdImageMetadataInterface.h | 6 +
Code/IO/otbSpotImageMetadataInterface.cxx | 102 ++++++++++++
Code/IO/otbSpotImageMetadataInterface.h | 6 +
.../otbAtmosphericCorrectionParameters.cxx | 61 +++++--
.../otbAtmosphericCorrectionParameters.h | 22 ++-
...arametersTo6SAtmosphericRadiativeTerms.cxx | 4 +-
.../otbAtmosphericRadiativeTerms.cxx | 9 ++
...flectanceToSurfaceReflectanceImageFilter.h | 25 +++
...ectanceToSurfaceReflectanceImageFilter.txx | 152 ++++++++++++++++++
14 files changed, 583 insertions(+), 15 deletions(-)
create mode 100644 Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.txx
diff --git a/Code/IO/otbDefaultImageMetadataInterface.h b/Code/IO/otbDefaultImageMetadataInterface.h
index 33d56ae6df..4c50d84104 100644
--- a/Code/IO/otbDefaultImageMetadataInterface.h
+++ b/Code/IO/otbDefaultImageMetadataInterface.h
@@ -130,6 +130,19 @@ public:
{
itkExceptionMacro("GetSatElevation not implemented in DefaultImageMetadataInterface, no captor type found");
};
+
+
+ /** Get the first wavelength for the spectral band definition */
+ VariableLengthVectorType GetFirstWavelengths( const MetaDataDictionaryType & dict ) const
+ {
+ itkExceptionMacro("GetFirstWavelengths not implemented in DefaultImageMetadataInterface, no captor type found");
+ };
+
+ /** Get the last wavelength for the spectral band definition */
+ VariableLengthVectorType GetLastWavelengths( const MetaDataDictionaryType & dict ) const
+ {
+ itkExceptionMacro("GetLastWavelengths not implemented in DefaultImageMetadataInterface, no captor type found");
+ };
bool CanRead( const MetaDataDictionaryType & dict ) const
{
diff --git a/Code/IO/otbIkonosImageMetadataInterface.cxx b/Code/IO/otbIkonosImageMetadataInterface.cxx
index 49792cd182..a11002f671 100755
--- a/Code/IO/otbIkonosImageMetadataInterface.cxx
+++ b/Code/IO/otbIkonosImageMetadataInterface.cxx
@@ -522,4 +522,93 @@ IkonosImageMetadataInterface::GetSatAzimuth( const MetaDataDictionaryType & dict
}
+IkonosImageMetadataInterface::VariableLengthVectorType
+IkonosImageMetadataInterface
+::GetFirstWavelengths( const MetaDataDictionaryType & dict ) const
+{
+ if( !this->CanRead( dict ) )
+ {
+ itkExceptionMacro(<<"Invalid Metadata, no Ikonos Image");
+ }
+
+ ImageKeywordlistType imageKeywordlist;
+
+ if (dict.HasKey(MetaDataKey::OSSIMKeywordlistKey))
+ {
+ itk::ExposeMetaData<ImageKeywordlistType>(dict, MetaDataKey::OSSIMKeywordlistKey, imageKeywordlist);
+ }
+
+ VariableLengthVectorType wavel(1);
+ wavel.Fill(0.);
+
+ ossimKeywordlist kwl;
+ imageKeywordlist.convertToOSSIMKeywordlist(kwl);
+ std::string key= "support_data.number_bands";
+ int nbBands = ossimString(kwl.find(key.c_str())).toInt();
+
+ // Panchromatic case
+ if(nbBands==1)
+ {
+ wavel.SetSize(1);
+ wavel.Fill(0.526);
+ }
+ else if(nbBands==4)
+ {
+ wavel.SetSize(4);
+ wavel[0] = 0.445;
+ wavel[1] = 0.506;
+ wavel[2] = 0.632;
+ wavel[3] = 0.757;
+ }
+ else
+ itkExceptionMacro(<<"Invalid number of bands...");
+
+ return wavel;
+}
+
+
+IkonosImageMetadataInterface::VariableLengthVectorType
+IkonosImageMetadataInterface
+::GetLastWavelengths( const MetaDataDictionaryType & dict ) const
+{
+ if( !this->CanRead( dict ) )
+ {
+ itkExceptionMacro(<<"Invalid Metadata, no Ikonos Image");
+ }
+
+ ImageKeywordlistType imageKeywordlist;
+
+ if (dict.HasKey(MetaDataKey::OSSIMKeywordlistKey))
+ {
+ itk::ExposeMetaData<ImageKeywordlistType>(dict, MetaDataKey::OSSIMKeywordlistKey, imageKeywordlist);
+ }
+
+ VariableLengthVectorType wavel(1);
+ wavel.Fill(0.);
+
+ ossimKeywordlist kwl;
+ imageKeywordlist.convertToOSSIMKeywordlist(kwl);
+ std::string key= "support_data.number_bands";
+ int nbBands = ossimString(kwl.find(key.c_str())).toInt();
+
+ // Panchromatic case
+ if(nbBands==1)
+ {
+ wavel.SetSize(1);
+ wavel.Fill(0.929);
+ }
+ else if(nbBands==4)
+ {
+ wavel.SetSize(4);
+ wavel[0] = 0.516;
+ wavel[1] = 0.595;
+ wavel[2] = 0.698;
+ wavel[3] = 0.853;
+ }
+ else
+ itkExceptionMacro(<<"Invalid number of bands...");
+
+ return wavel;
+}
+
} // end namespace otb
diff --git a/Code/IO/otbIkonosImageMetadataInterface.h b/Code/IO/otbIkonosImageMetadataInterface.h
index 161d9c7384..a8d091570d 100644
--- a/Code/IO/otbIkonosImageMetadataInterface.h
+++ b/Code/IO/otbIkonosImageMetadataInterface.h
@@ -93,6 +93,12 @@ public:
/** Get the sat azimuth from the ossim metadata */
double GetSatAzimuth( const MetaDataDictionaryType & dict ) const;
+ /** Get the first wavelength for the spectral band definition */
+ VariableLengthVectorType GetFirstWavelengths( const MetaDataDictionaryType & dict ) const;
+
+ /** Get the last wavelength for the spectral band definition */
+ VariableLengthVectorType GetLastWavelengths( const MetaDataDictionaryType & dict ) const;
+
bool CanRead( const MetaDataDictionaryType & dict) const;
diff --git a/Code/IO/otbImageMetadataInterfaceBase.h b/Code/IO/otbImageMetadataInterfaceBase.h
index 3a95fc0e97..3ddc902354 100644
--- a/Code/IO/otbImageMetadataInterfaceBase.h
+++ b/Code/IO/otbImageMetadataInterfaceBase.h
@@ -194,6 +194,14 @@ public:
virtual double GetSatAzimuth( const MetaDataDictionaryType & dict ) const =0;
otbMetadataGetMacro(SatAzimuth, double);
+ /** Get the first wavelength for the spectral band definition */
+ virtual VariableLengthVectorType GetFirstWavelengths( const MetaDataDictionaryType & dict ) const =0;
+ otbMetadataGetMacro(FirstWavelengths, VariableLengthVectorType);
+
+ /** Get the last wavelength for the spectral band definition */
+ virtual VariableLengthVectorType GetLastWavelengths( const MetaDataDictionaryType & dict ) const =0;
+ otbMetadataGetMacro(LastWavelengths, VariableLengthVectorType);
+
virtual bool CanRead( const MetaDataDictionaryType & dict ) const =0;
virtual void PrintSelf(std::ostream& os, itk::Indent indent, const MetaDataDictionaryType & dict) const;
diff --git a/Code/IO/otbQuickBirdImageMetadataInterface.cxx b/Code/IO/otbQuickBirdImageMetadataInterface.cxx
index 26d74ae29f..b5dd89f1d6 100755
--- a/Code/IO/otbQuickBirdImageMetadataInterface.cxx
+++ b/Code/IO/otbQuickBirdImageMetadataInterface.cxx
@@ -625,6 +625,101 @@ QuickBirdImageMetadataInterface::GetSatAzimuth( const MetaDataDictionaryType & d
return keywordString.toDouble();
}
+QuickBirdImageMetadataInterface::VariableLengthVectorType
+QuickBirdImageMetadataInterface
+::GetFirstWavelengths( const MetaDataDictionaryType & dict ) const
+{
+ if( !this->CanRead( dict ) )
+ {
+ itkExceptionMacro(<<"Invalid Metadata, no QuickBird Image");
+ }
+
+ ImageKeywordlistType imageKeywordlist;
+
+ if (dict.HasKey(MetaDataKey::OSSIMKeywordlistKey))
+ {
+ itk::ExposeMetaData<ImageKeywordlistType>(dict, MetaDataKey::OSSIMKeywordlistKey, imageKeywordlist);
+ }
+
+ VariableLengthVectorType wavel(1);
+ wavel.Fill(0.);
+
+ ossimKeywordlist kwl;
+ imageKeywordlist.convertToOSSIMKeywordlist(kwl);
+ std::string key = "support_data.band_id";
+ ossimString keywordStringBId = kwl.find(key.c_str());
+
+ if (keywordStringBId != ossimString("P") && keywordStringBId != ossimString("Multi"))
+ {
+ itkExceptionMacro(<<"Invalid bandID "<<keywordStringBId);
+ }
+
+ // Panchromatic case
+ if (keywordStringBId == ossimString("P") )
+ {
+ wavel.SetSize(1);
+ wavel.Fill(0.450);
+ }
+ else
+ {
+ wavel.SetSize(4);
+ wavel[0] = 0.450;
+ wavel[1] = 0.520;
+ wavel[2] = 0.630;
+ wavel[3] = 0.760;
+ }
+
+ return wavel;
+}
+
+
+QuickBirdImageMetadataInterface::VariableLengthVectorType
+QuickBirdImageMetadataInterface
+::GetLastWavelengths( const MetaDataDictionaryType & dict ) const
+{
+ if( !this->CanRead( dict ) )
+ {
+ itkExceptionMacro(<<"Invalid Metadata, no QuickBird Image");
+ }
+
+ ImageKeywordlistType imageKeywordlist;
+
+ if (dict.HasKey(MetaDataKey::OSSIMKeywordlistKey))
+ {
+ itk::ExposeMetaData<ImageKeywordlistType>(dict, MetaDataKey::OSSIMKeywordlistKey, imageKeywordlist);
+ }
+
+ VariableLengthVectorType wavel(1);
+ wavel.Fill(0.);
+
+ ossimKeywordlist kwl;
+ imageKeywordlist.convertToOSSIMKeywordlist(kwl);
+ std::string key = "support_data.band_id";
+ ossimString keywordStringBId = kwl.find(key.c_str());
+
+ if (keywordStringBId != ossimString("P") && keywordStringBId != ossimString("Multi"))
+ {
+ itkExceptionMacro(<<"Invalid bandID "<<keywordStringBId);
+ }
+
+ // Panchromatic case
+ if (keywordStringBId == ossimString("P") )
+ {
+ wavel.SetSize(1);
+ wavel.Fill(0.900);
+ }
+ else
+ {
+ wavel.SetSize(4);
+ wavel[0] = 0.520;
+ wavel[1] = 0.600;
+ wavel[2] = 0.690;
+ wavel[3] = 0.900;
+ }
+
+ return wavel;
+}
+
} // end namespace otb
diff --git a/Code/IO/otbQuickBirdImageMetadataInterface.h b/Code/IO/otbQuickBirdImageMetadataInterface.h
index 0c9000f92f..73c2c757a5 100644
--- a/Code/IO/otbQuickBirdImageMetadataInterface.h
+++ b/Code/IO/otbQuickBirdImageMetadataInterface.h
@@ -91,6 +91,12 @@ public:
/** Get the sat azimuth from the ossim metadata */
double GetSatAzimuth( const MetaDataDictionaryType & dict ) const;
+
+ /** Get the first wavelength for the spectral band definition */
+ VariableLengthVectorType GetFirstWavelengths( const MetaDataDictionaryType & dict ) const;
+
+ /** Get the last wavelength for the spectral band definition */
+ VariableLengthVectorType GetLastWavelengths( const MetaDataDictionaryType & dict ) const;
bool CanRead( const MetaDataDictionaryType & dict) const;
diff --git a/Code/IO/otbSpotImageMetadataInterface.cxx b/Code/IO/otbSpotImageMetadataInterface.cxx
index 2d03f6fa37..7b930df34e 100755
--- a/Code/IO/otbSpotImageMetadataInterface.cxx
+++ b/Code/IO/otbSpotImageMetadataInterface.cxx
@@ -532,5 +532,107 @@ SpotImageMetadataInterface::GetSatAzimuth( const MetaDataDictionaryType & dict )
return satAz;
}
+SpotImageMetadataInterface::VariableLengthVectorType
+SpotImageMetadataInterface
+::GetFirstWavelengths( const MetaDataDictionaryType & dict ) const
+{
+ if( !this->CanRead( dict ) )
+ {
+ itkExceptionMacro(<<"Invalid Metadata, no Spot Image");
+ }
+
+ ImageKeywordlistType imageKeywordlist;
+
+ if (dict.HasKey(MetaDataKey::OSSIMKeywordlistKey))
+ {
+ itk::ExposeMetaData<ImageKeywordlistType>(dict, MetaDataKey::OSSIMKeywordlistKey, imageKeywordlist);
+ }
+
+ VariableLengthVectorType wavel(1);
+ wavel.Fill(0.);
+
+
+ ossimKeywordlist kwl;
+ imageKeywordlist.convertToOSSIMKeywordlist(kwl);
+ std::string key= "support_data.number_bands";
+ int nbBands = ossimString(kwl.find(key.c_str())).toInt();
+ std::string sensorId = this->GetSensorID(dict);
+
+ // Panchromatic case
+ if(nbBands==1)
+ {
+ wavel.SetSize(1);
+ if(sensorId == "SPOT4")
+ wavel.Fill(0.610);
+ else if(sensorId == "SPOT5")
+ wavel.Fill(0.480);
+ else
+ itkExceptionMacro(<<"Invalid Spot Sensor ID");
+ }
+ else if(nbBands>1 && nbBands<5)
+ {
+ wavel.SetSize(4);
+ wavel[0] = 0.500;
+ wavel[1] = 0.610;
+ wavel[2] = 0.780;
+ wavel[3] = 1.580;
+ }
+ else
+ itkExceptionMacro(<<"Invalid number of bands...");
+
+ return wavel;
+}
+
+
+SpotImageMetadataInterface::VariableLengthVectorType
+SpotImageMetadataInterface
+::GetLastWavelengths( const MetaDataDictionaryType & dict ) const
+{
+ if( !this->CanRead( dict ) )
+ {
+ itkExceptionMacro(<<"Invalid Metadata, no Spot Image");
+ }
+
+ ImageKeywordlistType imageKeywordlist;
+
+ if (dict.HasKey(MetaDataKey::OSSIMKeywordlistKey))
+ {
+ itk::ExposeMetaData<ImageKeywordlistType>(dict, MetaDataKey::OSSIMKeywordlistKey, imageKeywordlist);
+ }
+
+ VariableLengthVectorType wavel(1);
+ wavel.Fill(0.);
+
+ ossimKeywordlist kwl;
+ imageKeywordlist.convertToOSSIMKeywordlist(kwl);
+ std::string key= "support_data.number_bands";
+ int nbBands = ossimString(kwl.find(key.c_str())).toInt();
+ std::string sensorId = this->GetSensorID(dict);
+
+ // Panchromatic case
+ if(nbBands==1)
+ {
+ wavel.SetSize(1);
+ if(sensorId == "SPOT4")
+ wavel.Fill(0.680);
+ else if(sensorId == "SPOT5")
+ wavel.Fill(0.710);
+ else
+ itkExceptionMacro(<<"Invalid Spot Sensor ID");
+ }
+ else if(nbBands>1 && nbBands<5)
+ {
+ wavel.SetSize(4);
+ wavel[0] = 0.590;
+ wavel[1] = 0.680;
+ wavel[2] = 0.890;
+ wavel[3] = 1.750;
+ }
+ else
+ itkExceptionMacro(<<"Invalid number of bands...");
+
+ return wavel;
+}
+
} // end namespace otb
diff --git a/Code/IO/otbSpotImageMetadataInterface.h b/Code/IO/otbSpotImageMetadataInterface.h
index f73e31e713..4d6b9672e1 100644
--- a/Code/IO/otbSpotImageMetadataInterface.h
+++ b/Code/IO/otbSpotImageMetadataInterface.h
@@ -93,6 +93,12 @@ public:
/** Get the sat azimuth from the ossim metadata */
double GetSatAzimuth( const MetaDataDictionaryType & dict ) const;
+ /** Get the first wavelength for the spectral band definition */
+ VariableLengthVectorType GetFirstWavelengths( const MetaDataDictionaryType & dict ) const;
+
+ /** Get the last wavelength for the spectral band definition */
+ VariableLengthVectorType GetLastWavelengths( const MetaDataDictionaryType & dict ) const;
+
bool CanRead( const MetaDataDictionaryType & dict) const;
diff --git a/Code/Radiometry/otbAtmosphericCorrectionParameters.cxx b/Code/Radiometry/otbAtmosphericCorrectionParameters.cxx
index f7539476f4..e99d56c0c3 100644
--- a/Code/Radiometry/otbAtmosphericCorrectionParameters.cxx
+++ b/Code/Radiometry/otbAtmosphericCorrectionParameters.cxx
@@ -18,6 +18,12 @@ PURPOSE. See the above copyright notices for more information.
#include "otbAtmosphericCorrectionParameters.h"
+#include "otbAeronetFileReader.h"
+#include <fstream>
+#include <iostream>
+
+
+
namespace otb
{
/*********************** FilterFunctionValues **************************/
@@ -26,6 +32,7 @@ FilterFunctionValues
{
m_MinSpectralValue = 0;
m_MaxSpectralValue = 0;
+ m_UserStep = 0.0025;
m_FilterFunctionValues.clear();
}
@@ -59,20 +66,56 @@ FilterFunctionValues
AtmosphericCorrectionParameters
::AtmosphericCorrectionParameters()
{
- m_SolarZenithalAngle = 361.;
- m_SolarAzimutalAngle = 361.;
+ m_SolarZenithalAngle = 361.;
+ m_SolarAzimutalAngle = 361.;
m_ViewingZenithalAngle = 361.;
m_ViewingAzimutalAngle = 361.;
- m_Month = 0;
- m_Day = 0;
- m_AtmosphericPressure = 1030.;
- m_WaterVaporAmount = 2.5;
- m_OzoneAmount = 0.28;
- m_AerosolModel = CONTINENTAL;
- m_AerosolOptical = 0.2;
+ m_Month = 0;
+ m_Day = 0;
+ m_AtmosphericPressure = 1030.;
+ m_WaterVaporAmount = 2.5;
+ m_OzoneAmount = 0.28;
+ m_AerosolModel = CONTINENTAL;
+ m_AerosolOptical = 0.2;
+}
+/** Get data from aeronet file*/
+void
+AtmosphericCorrectionParameters
+::UpdateAeronetData( std::string file, int year, int month, int day, int hour, int minute, double epsi )
+{
+ if(file == "")
+ itkExceptionMacro(<<"No Aeronet filename specified.");
+
+ AeronetFileReader::Pointer reader = AeronetFileReader::New();
+ reader->SetFileName(file);
+ reader->SetDay(day);
+ reader->SetMonth(month);
+ reader->SetYear(year);
+ reader->SetHour(hour);
+ reader->SetMinute(minute);
+ reader->SetEpsilon(epsi);
+ std::cout<<day<<std::endl;
+ std::cout<<month<<std::endl;
+ std::cout<<year<<std::endl;
+ std::cout<<hour<<std::endl;
+ std::cout<<minute<<std::endl;
+ std::cout<<epsi<<std::endl;
+
+ std::cout<<reader->GetDay()<<std::endl;
+ std::cout<<reader->GetMonth()<<std::endl;
+ std::cout<<reader->GetYear()<<std::endl;
+ std::cout<<reader->GetHour()<<std::endl;
+ std::cout<<reader->GetMinute()<<std::endl;
+ std::cout<<reader->GetEpsilon()<<std::endl;
+
+ reader->Update();
+
+ m_AerosolOptical = reader->GetOutput()->GetAerosolOpticalThickness();
+ m_WaterVaporAmount = reader->GetOutput()->GetWater();
}
+
/**PrintSelf method */
void
AtmosphericCorrectionParameters
diff --git a/Code/Radiometry/otbAtmosphericCorrectionParameters.h b/Code/Radiometry/otbAtmosphericCorrectionParameters.h
index 79fff9db12..353bbe6cc4 100644
--- a/Code/Radiometry/otbAtmosphericCorrectionParameters.h
+++ b/Code/Radiometry/otbAtmosphericCorrectionParameters.h
@@ -90,6 +90,7 @@ public:
/** Get user step between each wavelenght spectral band values. */
itkGetMacro(UserStep,WavelenghtSpectralBandType);
+
protected:
/** Constructor */
FilterFunctionValues();
@@ -103,8 +104,7 @@ protected:
private:
FilterFunctionValues(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
-
-
+
/** Vector that contains the filter function value. */
ValuesVectorType m_FilterFunctionValues;
/** Vector that contains the filter function value in 6S format (step of 0.0025µm).
@@ -242,6 +242,22 @@ public:
return &m_WavelenghtSpectralBand;
};
+ /** Read the aeronet data and extract aerosol optical and water vapor amount. */
+ void UpdateAeronetData( std::string file, int year, int month, int day, int hour, int minute, double epsi );
+ void UpdateAeronetData( std::string file, int year, int month, int day, int hour, int minute )
+ {
+ this->UpdateAeronetData( file, year, month, day, hour, minute, 0.4 );
+ };
+ void UpdateAeronetData( std::string file, int year, int hour, int minute, double epsi )
+ {
+ this->UpdateAeronetData( file, year, m_Month, m_Day, hour, minute, epsi );
+ };
+ void UpdateAeronetData( std::string file, int year, int hour, int minute )
+ {
+ this->UpdateAeronetData( file, year, m_Month, m_Day, hour, minute, 0.4 );
+ };
+
+
/** Constructor */
AtmosphericCorrectionParameters();
/** Destructor */
@@ -278,10 +294,8 @@ private:
AerosolModelType m_AerosolModel;
/** The Aerosol optical (radiative impact of aerosol for the reference wavelenght 550-nm) */
double m_AerosolOptical;
-
/** Wavelenght for the each spectral band definition */
WavelenghtSpectralBandVectorType m_WavelenghtSpectralBand;
-
};
} // end namespace otb
diff --git a/Code/Radiometry/otbAtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms.cxx b/Code/Radiometry/otbAtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms.cxx
index 8fc75cc5c7..05ee358b6f 100644
--- a/Code/Radiometry/otbAtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms.cxx
+++ b/Code/Radiometry/otbAtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms.cxx
@@ -155,10 +155,10 @@ AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms
double upwardDirectTransmittance(0.);
double upwardDiffuseTransmittanceForRayleigh(0.);
double upwardDiffuseTransmittanceForAerosol(0.);
-
+
for (unsigned int i=0; i<NbBand; ++i)
{
- atmosphericReflectance = 0.;
+ atmosphericReflectance = 0.;
atmosphericSphericalAlbedo = 0.;
totalGaseousTransmission = 0.;
downwardTransmittance = 0.;
diff --git a/Code/Radiometry/otbAtmosphericRadiativeTerms.cxx b/Code/Radiometry/otbAtmosphericRadiativeTerms.cxx
index e3540908cf..9f615440d7 100644
--- a/Code/Radiometry/otbAtmosphericRadiativeTerms.cxx
+++ b/Code/Radiometry/otbAtmosphericRadiativeTerms.cxx
@@ -505,6 +505,15 @@ AtmosphericRadiativeTerms
{
os << indent << "Channel "<< i << " : "<< std::endl;
//ValueType::(os,indent);
+ os << indent << "Intrinsic Atmospheric Reflectance : " << m_Values[i]->GetIntrinsicAtmosphericReflectance() << std::endl;
+ os << indent << "Shperical Albedo of the Atmosphere : " << m_Values[i]->GetSphericalAlbedo() << std::endl;
+ os << indent << "Total Gaseous Transmission : " << m_Values[i]->GetTotalGaseousTransmission() << std::endl;
+ os << indent << "Downward Transmittance of the Atmospher : " << m_Values[i]->GetDownwardTransmittance() << std::endl;
+ os << indent << "Upward Transmittance of the Atmospher : " << m_Values[i]->GetUpwardTransmittance() << std::endl;
+ os << indent << "Upward diffuse transmittance : " << m_Values[i]->GetUpwardDiffuseTransmittance() << std::endl;
+ os << indent << "Upward direct transmittance : " << m_Values[i]->GetUpwardDirectTransmittance() << std::endl;
+ os << indent << "Upward diffuse transmittance for rayleigh: " << m_Values[i]->GetUpwardDiffuseTransmittanceForRayleigh() << std::endl;
+ os << indent << "Upward diffuse transmittance for aerosols: " << m_Values[i]->GetUpwardDiffuseTransmittanceForAerosol() << std::endl;
}
}
diff --git a/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.h b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.h
index 7c2a41be29..8111427e48 100644
--- a/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.h
+++ b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.h
@@ -158,6 +158,10 @@ public:
typedef AtmosphericCorrectionParameters::Pointer CorrectionParametersPointerType;
typedef AtmosphericRadiativeTerms::Pointer AtmosphericRadiativeTermsPointerType;
+ typedef FilterFunctionValues FilterFunctionValuesType;
+ typedef FilterFunctionValuesType::ValuesVectorType CoefVectorType;
+ typedef std::vector<CoefVectorType> FilterFunctionCoefVectorType;
+
typedef itk::MetaDataDictionary MetaDataDictionaryType;
/** Get/Set Atmospheric Radiative Terms. */
@@ -177,6 +181,21 @@ public:
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_FilterFunctionCoef;
+ }
+
protected:
/** Constructor */
ReflectanceToSurfaceReflectanceImageFilter();
@@ -199,6 +218,12 @@ private:
bool m_IsSetAtmosphericRadiativeTerms;
/** 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;
+ /** BeforeThreadedGenerateData executed once or not */
+ bool m_BeforeDone;
};
} // end namespace otb
diff --git a/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.txx b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.txx
new file mode 100644
index 0000000000..dfd308f360
--- /dev/null
+++ b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter.txx
@@ -0,0 +1,152 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef __otbReflectanceToSurfaceReflectanceImageFilter_txx
+#define __otbReflectanceToSurfaceReflectanceImageFilter_txx
+
+#include "otbReflectanceToSurfaceReflectanceImageFilter.h"
+#include "otbImageMetadataInterfaceFactory.h"
+#include "otbImageMetadataInterfaceBase.h"
+#include "otbAeronetFileReader.h"
+
+namespace otb
+{
+
+/**
+ * Constructor
+ */
+template <class TInputImage, class TOutputImage>
+ReflectanceToSurfaceReflectanceImageFilter<TInputImage,TOutputImage>
+::ReflectanceToSurfaceReflectanceImageFilter()
+{
+ m_AtmosphericRadiativeTerms = AtmosphericRadiativeTerms::New();
+ m_CorrectionParameters = AtmosphericCorrectionParameters::New();
+ m_IsSetAtmosphericRadiativeTerms = false;
+ m_BeforeDone = false;
+ m_AeronetFileName = "";
+ m_FilterFunctionValuesFileName = "";
+ m_FilterFunctionCoef.clear();
+}
+
+
+template <class TInputImage, class TOutputImage>
+void
+ReflectanceToSurfaceReflectanceImageFilter<TInputImage,TOutputImage>
+::UpdateAtmosphericRadiativeTerms( const MetaDataDictionaryType dict )
+{
+ 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_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) );
+
+ Parameters2RadiativeTermsPointerType param2Terms = Parameters2RadiativeTermsType::New();
+
+ 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);
+ }
+
+ param2Terms->SetInput(m_CorrectionParameters);
+ param2Terms->Update();
+
+ m_AtmosphericRadiativeTerms = param2Terms->GetOutput();
+}
+
+
+template <class TInputImage, class TOutputImage>
+void
+ReflectanceToSurfaceReflectanceImageFilter<TInputImage,TOutputImage>
+::BeforeThreadedGenerateData()
+{
+ if(m_BeforeDone == true)
+ return;
+
+ std::cout<<"BeforeThreadedGenerateData "<<this->GetNumberOfThreads()<<std::endl;
+ if(m_IsSetAtmosphericRadiativeTerms==false)
+ this->UpdateAtmosphericRadiativeTerms(this->GetInput()->GetMetaDataDictionary());
+
+ std::cout<<"m_CorrectionParameters: "<<std::endl;
+ std::cout<<m_CorrectionParameters<<std::endl;
+ std::cout<<"m_AtmosphericRadiativeTerms: "<<std::endl;
+ std::cout<<m_AtmosphericRadiativeTerms<<std::endl;
+
+ this->GetFunctorVector().clear();
+ for (unsigned int i = 0;i<this->GetInput()->GetNumberOfComponentsPerPixel();++i)
+ {
+ double coef;
+ double res;
+ coef = static_cast<double>(m_AtmosphericRadiativeTerms->GetTotalGaseousTransmission(i)
+ * m_AtmosphericRadiativeTerms->GetDownwardTransmittance(i)
+ * m_AtmosphericRadiativeTerms->GetUpwardTransmittance(i) );
+ coef = 1. / coef;
+ res = -m_AtmosphericRadiativeTerms->GetIntrinsicAtmosphericReflectance(i) * coef;
+
+ FunctorType functor;
+ functor.SetCoefficient(coef);
+ functor.SetResidu(res);
+ functor.SetSphericalAlbedo(static_cast<double>(m_AtmosphericRadiativeTerms->GetSphericalAlbedo(i)));
+
+ this->GetFunctorVector().push_back(functor);
+ }
+ m_BeforeDone = true;
+ std::cout<<"BeforeThreadedGenerateData FIN"<<std::endl;
+}
+
+}
+
+#endif
--
GitLab