diff --git a/Code/Radiometry/otbAtmosphericRadiativeTerms2.h b/Code/Radiometry/otbAtmosphericRadiativeTerms2.h index 69227c57dcbb3c6f0fa51a7ab87734dcec316457..0705195d57ef8694b72f01fcb1124c3111a3d334 100644 --- a/Code/Radiometry/otbAtmosphericRadiativeTerms2.h +++ b/Code/Radiometry/otbAtmosphericRadiativeTerms2.h @@ -19,6 +19,7 @@ #define __otbAtmosphericRadiativeTerms_h #include "itkObject.h" +#include "itkDataObject.h" #include "itkObjectFactory.h" #include "itkMacro.h" #include <vector> @@ -160,7 +161,7 @@ private: * \ingroup Radiometry */ -class ITK_EXPORT AtmosphericRadiativeTerms2 : public itk::Object +class ITK_EXPORT AtmosphericRadiativeTerms2 : public itk::DataObject { public: /** Standard typedefs */ @@ -273,4 +274,4 @@ private: } // end namespace otb -#endif \ No newline at end of file +#endif diff --git a/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter2.h b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter2.h new file mode 100644 index 0000000000000000000000000000000000000000..9f30fdb6b200a87203658d4b65a47770f72943fd --- /dev/null +++ b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter2.h @@ -0,0 +1,279 @@ +/*========================================================================= + + Program: ORFEO Toolbox + Language: C++ + Date: $Date$ + Version: $Revision$ + + + Copyright (c) Centre National d'Etudes Spatiales. All rights reserved. + See OTBCopyright.txt for details. + + + Some parts of this code are derived from ITK. See ITKCopyright.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 __otbReflectanceToSurfaceReflectanceImageFilter2_h +#define __otbReflectanceToSurfaceReflectanceImageFilter2_h + +#include "otbUnaryImageFunctorWithVectorImageFilter.h" + +#include "otbRadiometryCorrectionParametersToAtmosphericRadiativeTerms.h" +#include "otbAtmosphericCorrectionParameters2.h" +#include "otbImageMetadataCorrectionParameters.h" +#include "otbAtmosphericRadiativeTerms2.h" + +#include "itkMetaDataDictionary.h" + +#include <iomanip> + +namespace otb +{ +namespace Functor +{ +/** + * \class ReflectanceToSurfaceReflectanceImageFunctor + * \brief Compute the surface reflectance pixel from a TOA reflectance. + * + * \ingroup Functor + * \ingroup Radiometry + * + */ +template <class TInput, class TOutput> +class ReflectanceToSurfaceReflectanceImageFunctor +{ +public: + ReflectanceToSurfaceReflectanceImageFunctor() + { + m_Coefficient = 1.; + m_Residu = 1.; + m_SphericalAlbedo = 1.; + } + virtual ~ReflectanceToSurfaceReflectanceImageFunctor() {} + + /** + * Set/Get the spherical albedo of the atmosphere. + */ + void SetSphericalAlbedo(double albedo) + { + m_SphericalAlbedo = albedo; + } + double GetSphericalAlbedo() + { + return m_SphericalAlbedo; + } + + /** + * Set/Get Coefficient, computed from AtmosphericRadiativeTermsPointerType datas. + */ + void SetCoefficient(double coef) + { + m_Coefficient = coef; + } + double GetCoefficient() + { + return m_Coefficient; + } + + /** + * Set/Get Residu, computed from AtmosphericRadiativeTermsPointerType datas. + */ + void SetResidu(double res) + { + m_Residu = res; + } + double GetResidu() + { + return m_Residu; + } + + inline TOutput operator ()(const TInput& inPixel) + { + TOutput outPixel; + double temp, temp2; + temp = (static_cast<double>(inPixel) + m_Residu)* m_Coefficient; + temp2 = temp / (1. + m_SphericalAlbedo * temp); + + outPixel = static_cast<TOutput>(temp2); + + return outPixel; + } + +private: + double m_SphericalAlbedo; + double m_Coefficient; + double m_Residu; +}; +} + +/** \class ReflectanceToSurfaceReflectanceImageFilter2 + * \brief Calculates the slope, the orientation incidence and exitance radius values for each pixel. + * + * + * \ingroup AtmosphericRadiativeTerms + * \ingroup AtmosphericCorrectionParameters + * \ingroup LuminanceToReflectanceImageFilter + * \ingroup ImageToReflectanceImageFilter + * \ingroup ImageMetadataCorrectionParameters + * \ingroup otbRadiometryCorrectionParametersToAtmosphericRadiativeTerms + * \ingroup Radiometry + * + * \example Radiometry/AtmosphericCorrectionSequencement.cxx + */ +template <class TInputImage, class TOutputImage> +class ITK_EXPORT ReflectanceToSurfaceReflectanceImageFilter2 : + public UnaryImageFunctorWithVectorImageFilter<TInputImage, + TOutputImage, + typename Functor::ReflectanceToSurfaceReflectanceImageFunctor< + typename TInputImage::InternalPixelType, + typename + TOutputImage::InternalPixelType> > +{ +public: + /** Extract input and output images dimensions.*/ + itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension); + itkStaticConstMacro(OutputImageDimension, unsigned int, TOutputImage::ImageDimension); + + /** "typedef" to simplify the variables definition and the declaration. */ + typedef TInputImage InputImageType; + typedef TOutputImage OutputImageType; + + typedef typename Functor::ReflectanceToSurfaceReflectanceImageFunctor<typename InputImageType::InternalPixelType, + typename OutputImageType::InternalPixelType> + FunctorType; + /** "typedef" for standard classes. */ + typedef ReflectanceToSurfaceReflectanceImageFilter2 Self; + typedef UnaryImageFunctorWithVectorImageFilter<InputImageType, OutputImageType, FunctorType> Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + + /** object factory method. */ + itkNewMacro(Self); + + /** return class name. */ + itkTypeMacro(ReflectanceToSurfaceReflectanceImageFilter2, UnaryImageFunctorWithVectorImageFilter); + + /** Supported images definition. */ + typedef typename InputImageType::PixelType InputPixelType; + typedef typename InputImageType::InternalPixelType InputInternalPixelType; + typedef typename InputImageType::RegionType InputImageRegionType; + typedef typename OutputImageType::PixelType OutputPixelType; + typedef typename OutputImageType::InternalPixelType OutputInternalPixelType; + typedef typename OutputImageType::RegionType OutputImageRegionType; + + + typedef otb::RadiometryCorrectionParametersToAtmosphericRadiativeTerms CorrectionParametersToRadiativeTermsType; + + typedef otb::AtmosphericCorrectionParameters2 AtmoCorrectionParametersType; + typedef typename AtmoCorrectionParametersType::Pointer AtmoCorrectionParametersPointerType; + + typedef otb::ImageMetadataCorrectionParameters AcquiCorrectionParametersType; + typedef typename AcquiCorrectionParametersType::Pointer AcquiCorrectionParametersPointerType; + + typedef otb::AtmosphericRadiativeTerms2 AtmosphericRadiativeTermsType; + typedef typename AtmosphericRadiativeTermsType::Pointer AtmosphericRadiativeTermsPointerType; + + + typedef otb::FilterFunctionValues FilterFunctionValuesType; + typedef FilterFunctionValuesType::WavelengthSpectralBandType ValueType; //float + typedef FilterFunctionValuesType::ValuesVectorType ValuesVectorType; //std::vector<float> + + typedef typename AcquiCorrectionParametersType::WavelengthSpectralBandVectorType WavelengthSpectralBandVectorType; + + + typedef itk::MetaDataDictionary MetaDataDictionaryType; + + /** Get/Set Atmospheric Radiative Terms. */ + void SetAtmosphericRadiativeTerms(AtmosphericRadiativeTermsPointerType atmoRadTerms) + { + m_AtmosphericRadiativeTerms = atmoRadTerms; + this->SetNthInput(1, m_AtmosphericRadiativeTerms); + m_IsSetAtmosphericRadiativeTerms = true; + this->Modified(); + } + itkGetObjectMacro(AtmosphericRadiativeTerms, AtmosphericRadiativeTermsType); + + /** Get/Set Atmospheric Correction Parameters. */ + void SetAtmoCorrectionParameters(AtmoCorrectionParametersPointerType atmoCorrTerms) + { + m_AtmoCorrectionParameters = atmoCorrTerms; + this->SetNthInput(2, m_AtmoCorrectionParameters); + m_IsSetAtmoCorrectionParameters = true; + this->Modified(); + } + itkGetObjectMacro(AtmoCorrectionParameters, AtmoCorrectionParametersType); + + /** Get/Set Acquisition Correction Parameters. */ + void SetAcquiCorrectionParameters(AcquiCorrectionParametersPointerType acquiCorrTerms) + { + m_AcquiCorrectionParameters = acquiCorrTerms; + this->SetNthInput(3, m_AcquiCorrectionParameters); + m_IsSetAcquiCorrectionParameters = true; + this->Modified(); + } + itkGetObjectMacro(AcquiCorrectionParameters, AcquiCorrectionParametersType); + + + /** Generate radiative terms from the atmospheric parameters */ + void GenerateAtmosphericRadiativeTerms(); + + /** Compute radiative terms if necessary and then update functors attributs. */ + void GenerateParameters(); + + /** Set/Get UseGenerateParameters. */ + itkSetMacro(UseGenerateParameters, bool); + itkGetMacro(UseGenerateParameters, bool); + itkBooleanMacro(UseGenerateParameters); + + /** Set/Get IsSetAtmosphericRadiativeTerms */ + itkSetMacro(IsSetAtmosphericRadiativeTerms, bool); + itkGetMacro(IsSetAtmosphericRadiativeTerms, bool); + itkBooleanMacro(IsSetAtmosphericRadiativeTerms); + +protected: + /** Constructor */ + ReflectanceToSurfaceReflectanceImageFilter2(); + /** Destructor */ + virtual ~ReflectanceToSurfaceReflectanceImageFilter2() {} + void PrintSelf(std::ostream& os, itk::Indent indent) const; + + + + /** Initialize the functor vector */ + void GenerateOutputInformation(); + /** Fill AtmosphericRadiativeTerms using image metadata*/ + void UpdateAtmosphericRadiativeTerms(); + /** Update Functors parameters */ + void UpdateFunctors(); + +private: + + bool m_IsSetAtmosphericRadiativeTerms; + bool m_IsSetAtmoCorrectionParameters; + bool m_IsSetAcquiCorrectionParameters; + + /** Enable/Disable GenerateParameters in GenerateOutputInformation. + * Useful for image view that call GenerateOutputInformation each time you move the full area. + */ + bool m_UseGenerateParameters; + + AtmosphericRadiativeTermsPointerType m_AtmosphericRadiativeTerms; + AtmoCorrectionParametersPointerType m_AtmoCorrectionParameters; + AcquiCorrectionParametersPointerType m_AcquiCorrectionParameters; + + +}; + +} // end namespace otb + +#ifndef OTB_MANUAL_INSTANTIATION +#include "otbReflectanceToSurfaceReflectanceImageFilter2.txx" +#endif + +#endif diff --git a/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter2.txx b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter2.txx new file mode 100644 index 0000000000000000000000000000000000000000..7a4b9af39b86f58592fc5812db0fa03f873bdf9c --- /dev/null +++ b/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter2.txx @@ -0,0 +1,193 @@ +/*========================================================================= + + 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 __otbReflectanceToSurfaceReflectanceImageFilter2_txx +#define __otbReflectanceToSurfaceReflectanceImageFilter2_txx + +#include "otbReflectanceToSurfaceReflectanceImageFilter2.h" +#include "otbOpticalImageMetadataInterfaceFactory.h" +#include "otbOpticalImageMetadataInterface.h" + +namespace otb +{ + +/** + * Constructor + */ +template <class TInputImage, class TOutputImage> +ReflectanceToSurfaceReflectanceImageFilter2<TInputImage, TOutputImage> +::ReflectanceToSurfaceReflectanceImageFilter2() : + m_IsSetAtmosphericRadiativeTerms(false), + m_IsSetAtmoCorrectionParameters(false), + m_IsSetAcquiCorrectionParameters(false), + m_UseGenerateParameters(true), + m_AtmosphericRadiativeTerms(NULL), + m_AtmoCorrectionParameters(NULL), + m_AcquiCorrectionParameters(NULL) + { + + } + +template <class TInputImage, class TOutputImage> +void +ReflectanceToSurfaceReflectanceImageFilter2<TInputImage, TOutputImage> +::UpdateAtmosphericRadiativeTerms() + { + if (!m_IsSetAtmosphericRadiativeTerms) + this->GenerateAtmosphericRadiativeTerms(); + + } + +template <class TInputImage, class TOutputImage> +void +ReflectanceToSurfaceReflectanceImageFilter2<TInputImage, TOutputImage> +::GenerateOutputInformation() + { + Superclass::GenerateOutputInformation(); + if (m_UseGenerateParameters) this->GenerateParameters(); + } + +template <class TInputImage, class TOutputImage> +void +ReflectanceToSurfaceReflectanceImageFilter2<TInputImage, TOutputImage> +::GenerateAtmosphericRadiativeTerms() + { + if (this->GetInput() == NULL) + { + itkExceptionMacro(<< "Input must be set before updating the atmospheric radiative terms"); + } + + if (!m_IsSetAtmoCorrectionParameters) + { + itkExceptionMacro(<< "Atmospheric correction parameters must be provided before updating the atmospheric radiative terms"); + } + + + // Acquisition parameters + if (!m_IsSetAcquiCorrectionParameters) // Get info from image metadata interface + { + MetaDataDictionaryType dict = this->GetInput()->GetMetaDataDictionary(); + OpticalImageMetadataInterface::Pointer imageMetadataInterface = OpticalImageMetadataInterfaceFactory::CreateIMI(dict); + + m_AcquiCorrectionParameters = AcquiCorrectionParametersType::New(); + + m_AcquiCorrectionParameters->SetSolarZenithalAngle(90. - imageMetadataInterface->GetSunElevation()); + m_AcquiCorrectionParameters->SetSolarAzimutalAngle(imageMetadataInterface->GetSunAzimuth()); + m_AcquiCorrectionParameters->SetViewingZenithalAngle(90. - imageMetadataInterface->GetSatElevation()); + m_AcquiCorrectionParameters->SetViewingAzimutalAngle(imageMetadataInterface->GetSatAzimuth()); + + m_AcquiCorrectionParameters->SetDay(imageMetadataInterface->GetDay()); + m_AcquiCorrectionParameters->SetMonth(imageMetadataInterface->GetMonth()); + + if (imageMetadataInterface->GetSpectralSensitivity()->Capacity() > 0) + { + m_AcquiCorrectionParameters->SetWavelengthSpectralBand(imageMetadataInterface->GetSpectralSensitivity()); + + } + else + { + otbMsgDevMacro(<< "use dummy filter"); + WavelengthSpectralBandVectorType spectralDummy; + spectralDummy->Clear(); + for (unsigned int i = 0; i < this->GetInput()->GetNumberOfComponentsPerPixel(); ++i) + { + spectralDummy->PushBack(FilterFunctionValuesType::New()); + } + m_AcquiCorrectionParameters->SetWavelengthSpectralBand(spectralDummy); + } + + } + + m_AtmosphericRadiativeTerms = CorrectionParametersToRadiativeTermsType::Compute(m_AtmoCorrectionParameters,m_AcquiCorrectionParameters); + + } + +template <class TInputImage, class TOutputImage> +void +ReflectanceToSurfaceReflectanceImageFilter2<TInputImage, TOutputImage> +::UpdateFunctors() + { + + if (this->GetInput() == NULL) + { + itkExceptionMacro(<< "Input must be set before updating the functors"); + } + + MetaDataDictionaryType dict = this->GetInput()->GetMetaDataDictionary(); + OpticalImageMetadataInterface::Pointer imageMetadataInterface = OpticalImageMetadataInterfaceFactory::CreateIMI(dict); + + this->GetFunctorVector().clear(); + + for (unsigned int i = 0; i < this->GetInput()->GetNumberOfComponentsPerPixel(); ++i) + { + double coef; + double res; + // 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)); + coef = 1. / coef; + res = -m_AtmosphericRadiativeTerms->GetIntrinsicAtmosphericReflectance(wavelengthPosition); + FunctorType functor; + functor.SetCoefficient(coef); + functor.SetResidu(res); + 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()); + } + } + +template <class TInputImage, class TOutputImage> +void +ReflectanceToSurfaceReflectanceImageFilter2<TInputImage, TOutputImage> +::GenerateParameters() + { + if (!m_IsSetAtmosphericRadiativeTerms) + { + this->UpdateAtmosphericRadiativeTerms(); + m_IsSetAtmosphericRadiativeTerms = true; + } + + this->UpdateFunctors(); + } + + +/* Standard "PrintSelf" method */ +template <class TInputImage, class TOutputImage> +void +ReflectanceToSurfaceReflectanceImageFilter2<TInputImage, TOutputImage> +::PrintSelf(std::ostream& os, itk::Indent indent) const +{ + os << indent << "Atmospheric radiative terms : " << m_AtmosphericRadiativeTerms << std::endl; + os << indent << "Atmospheric correction terms : " << m_AtmoCorrectionParameters << std::endl; + os << indent << "Acquisition correction terms : " << m_AcquiCorrectionParameters << std::endl; +} + +} //end namespace otb + +#endif diff --git a/Testing/Code/Radiometry/CMakeLists.txt b/Testing/Code/Radiometry/CMakeLists.txt index ef3fa305797d4c9f6e566ad5d112b44e96028b24..6c29e5cd94789eeda044a026674fc9016aa1ca9c 100644 --- a/Testing/Code/Radiometry/CMakeLists.txt +++ b/Testing/Code/Radiometry/CMakeLists.txt @@ -887,6 +887,38 @@ add_test(raTvReflectanceToSurfaceReflectanceImageFilter ${RADIOMETRY_TESTS3} 3 3 3 3 # upward transmittance ) +# ------- otb::ReflectanceToSurfaceReflectanceImageFilter2 ------------------------------ +add_test(raTuReflectanceToSurfaceReflectanceImageFilter2New ${RADIOMETRY_TESTS3} + otbReflectanceToSurfaceReflectanceImageFilter2New +) +add_test(raTvReflectanceToSurfaceReflectanceImageFilter2 ${RADIOMETRY_TESTS3} + --compare-image ${EPSILON_12} ${BASELINE}/raTvReflectanceToSurfaceReflectanceImageFilter.tif + ${TEMP}/raTvReflectanceToSurfaceReflectanceImageFilter.tif + otbReflectanceToSurfaceReflectanceImageFilter2Test + ${INPUTDATA}/verySmallFSATSW.tif + ${TEMP}/raTvReflectanceToSurfaceReflectanceImageFilter.tif + # SAME VALUE FOR EACH CHANNEL + 1 1 1 1 # intrinsic atmospheric reflectance + 2 2 2 2 # spherical albedo of the atmosphere + 3 3 3 3 # total transmission + 2 2 2 2 # downward transmittance + 3 3 3 3 # upward transmittance + ) + +add_test(raTvReflectanceToSurfaceReflectanceImageFilter2 ${RADIOMETRY_TESTS3} + --compare-image ${EPSILON_12} ${BASELINE}/raTvReflectanceToSurfaceReflectanceImageFilter2.tif + ${TEMP}/raTvReflectanceToSurfaceReflectanceImageFilter2.tif + otbReflectanceToSurfaceReflectanceImageFilter2Test2 + ${INPUTDATA}/ToulouseExtract_WithGeom.tif + ${TEMP}/raTvReflectanceToSurfaceReflectanceImageFilter2.tif + 1013 + 1.424 + 0.344 + 1 + 0.1 + ) + + if(OTB_DATA_USE_LARGEINPUT) #this test was comment after the refactoring of astmospheric correction classes #//TODO Need to rewrite tests with baselines generated directly with 6S. @@ -2130,6 +2162,7 @@ otbSIXSTraitsTest.cxx otbSIXSTraitsComputeAtmosphericParameters.cxx otbAtmosphericRadiativeTermsTest.cxx otbReflectanceToSurfaceReflectanceImageFilterTest.cxx +otbReflectanceToSurfaceReflectanceImageFilter2Test.cxx otbSurfaceAdjacencyEffect6SCorrectionSchemeFilterNew.cxx otbSurfaceAdjacencyEffect6SCorrectionSchemeFilter.cxx otbRomaniaReflectanceToRomaniaSurfaceReflectanceImageFilter.cxx diff --git a/Testing/Code/Radiometry/otbRadiometryTests3.cxx b/Testing/Code/Radiometry/otbRadiometryTests3.cxx index a30358b01ee4133c2118b56c26483aaef8bd45df..5787818e0f3835335c5d835603976dd27afcc295 100644 --- a/Testing/Code/Radiometry/otbRadiometryTests3.cxx +++ b/Testing/Code/Radiometry/otbRadiometryTests3.cxx @@ -39,6 +39,9 @@ void RegisterTests() REGISTER_TEST(otbReflectanceToSurfaceReflectanceImageFilterNew); REGISTER_TEST(otbReflectanceToSurfaceReflectanceImageFilterTest); REGISTER_TEST(otbReflectanceToSurfaceReflectanceImageFilterTest2); + REGISTER_TEST(otbReflectanceToSurfaceReflectanceImageFilter2New); + REGISTER_TEST(otbReflectanceToSurfaceReflectanceImageFilter2Test); + REGISTER_TEST(otbReflectanceToSurfaceReflectanceImageFilter2Test2); REGISTER_TEST(otbSurfaceAdjacencyEffect6SCorrectionSchemeFilterNew); REGISTER_TEST(otbSurfaceAdjacencyEffect6SCorrectionSchemeFilter); REGISTER_TEST(otbRomaniaReflectanceToRomaniaSurfaceReflectanceImageFilter); diff --git a/Testing/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter2Test.cxx b/Testing/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter2Test.cxx new file mode 100644 index 0000000000000000000000000000000000000000..0163cd24b6a286ec67226429d7d0ca3dfae98e63 --- /dev/null +++ b/Testing/Code/Radiometry/otbReflectanceToSurfaceReflectanceImageFilter2Test.cxx @@ -0,0 +1,171 @@ +/*========================================================================= + + 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. + +=========================================================================*/ +#include "itkMacro.h" + +#include <fstream> +#include <iostream> + +#include "otbReflectanceToSurfaceReflectanceImageFilter2.h" +#include "otbVectorImage.h" +#include "otbImageFileReader.h" +#include "otbImageFileWriter.h" +#include "otbAtmosphericRadiativeTerms2.h" +#include "otbAtmosphericCorrectionParameters2.h" + +#include "itkMetaDataDictionary.h" +#include "otbOpticalImageMetadataInterfaceFactory.h" +#include "otbOpticalImageMetadataInterface.h" + +int otbReflectanceToSurfaceReflectanceImageFilter2New(int argc, char * argv[]) +{ + const unsigned int Dimension = 2; + typedef double PixelType; + typedef otb::VectorImage<PixelType, Dimension> InputImageType; + + typedef otb::ReflectanceToSurfaceReflectanceImageFilter2<InputImageType, + InputImageType> + ReflectanceToSurfaceReflectanceImageFilter2Type; + + // Instantiating object + ReflectanceToSurfaceReflectanceImageFilter2Type::Pointer filter + = ReflectanceToSurfaceReflectanceImageFilter2Type::New(); + + return EXIT_SUCCESS; +} + +int otbReflectanceToSurfaceReflectanceImageFilter2Test(int argc, char * argv[]) +{ + const char * inputFileName = argv[1]; + const char * outputFileName = argv[2]; + + const unsigned int Dimension = 2; + typedef double PixelType; + typedef otb::VectorImage<PixelType, Dimension> InputImageType; + typedef otb::VectorImage<PixelType, Dimension> OutputImageType; + typedef otb::ImageFileReader<InputImageType> ReaderType; + typedef otb::ImageFileWriter<OutputImageType> WriterType; + + typedef otb::ReflectanceToSurfaceReflectanceImageFilter2<InputImageType, + OutputImageType> + ReflectanceToSurfaceReflectanceImageFilter2Type; + + typedef otb::AtmosphericRadiativeTerms2::DataVectorType DataVectorType; + otb::AtmosphericRadiativeTerms2::Pointer atmoRadTerms = otb::AtmosphericRadiativeTerms2::New(); + + + ReaderType::Pointer reader = ReaderType::New(); + WriterType::Pointer writer = WriterType::New(); + reader->SetFileName(inputFileName); + writer->SetFileName(outputFileName); + + reader->UpdateOutputInformation(); + unsigned int nbChannel = reader->GetOutput()->GetNumberOfComponentsPerPixel(); + + DataVectorType intrinsic; + DataVectorType albedo; + DataVectorType gaseous; + DataVectorType downTrans; + DataVectorType upTrans; + + for (unsigned int j = 0; j < nbChannel; ++j) + { + intrinsic.push_back(static_cast<double>(atof(argv[3 + j]))); + albedo.push_back(static_cast<double>(atof(argv[3 + j + nbChannel]))); + gaseous.push_back(static_cast<double>(atof(argv[3 + j + 2 * nbChannel]))); + downTrans.push_back(static_cast<double>(atof(argv[3 + j + 3 * nbChannel]))); + upTrans.push_back(static_cast<double>(atof(argv[3 + j + 4 * nbChannel]))); + } + + atmoRadTerms->SetIntrinsicAtmosphericReflectances(intrinsic); + atmoRadTerms->SetSphericalAlbedos(albedo); + atmoRadTerms->SetTotalGaseousTransmissions(gaseous); + atmoRadTerms->SetDownwardTransmittances(downTrans); + atmoRadTerms->SetUpwardTransmittances(upTrans); + + // Instantiating object + ReflectanceToSurfaceReflectanceImageFilter2Type::Pointer filter + = ReflectanceToSurfaceReflectanceImageFilter2Type::New(); + filter->SetAtmosphericRadiativeTerms(atmoRadTerms); + filter->SetInput(reader->GetOutput()); + writer->SetInput(filter->GetOutput()); + + writer->Update(); + + return EXIT_SUCCESS; +} + +//Check the correct generation of the atmospheric parameters +int otbReflectanceToSurfaceReflectanceImageFilter2Test2(int argc, char * argv[]) +{ + const char * inputFileName = argv[1]; + const char * outputFileName = argv[2]; + + const unsigned int Dimension = 2; + typedef double PixelType; + typedef otb::VectorImage<PixelType, Dimension> InputImageType; + typedef otb::VectorImage<PixelType, Dimension> OutputImageType; + typedef otb::ImageFileReader<InputImageType> ReaderType; + typedef otb::ImageFileWriter<OutputImageType> WriterType; + typedef otb::ReflectanceToSurfaceReflectanceImageFilter2<InputImageType, + OutputImageType> + ReflectanceToSurfaceReflectanceImageFilter2Type; + + typedef otb::AtmosphericCorrectionParameters2 AtmoCorrectionParametersType; + typedef AtmoCorrectionParametersType::AerosolModelType AerosolModelType; + AtmoCorrectionParametersType::Pointer paramAtmo = AtmoCorrectionParametersType::New(); + + + paramAtmo->SetAtmosphericPressure(static_cast<double>(atof(argv[3]))); + paramAtmo->SetWaterVaporAmount(static_cast<double>(atof(argv[4]))); + paramAtmo->SetOzoneAmount(static_cast<double>(atof(argv[5]))); + paramAtmo->SetAerosolModel(static_cast<AerosolModelType>(atoi(argv[6]))); + paramAtmo->SetAerosolOptical(static_cast<double>(atof(argv[7]))); + + + ReaderType::Pointer reader = ReaderType::New(); + WriterType::Pointer writer = WriterType::New(); + reader->SetFileName(inputFileName); + writer->SetFileName(outputFileName); + + reader->UpdateOutputInformation(); + + /*itk::MetaDataDictionary dict = reader->GetMetaDataDictionary(); + otb::OpticalImageMetadataInterface::Pointer imageMetadataInterface = otb::OpticalImageMetadataInterfaceFactory::CreateIMI(dict); + + std::cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< " << imageMetadataInterface->GetDay() << " <<<<<<<<<<<<<<<<<<<<<" << std::endl; + std::cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< " << imageMetadataInterface->GetMonth() << " <<<<<<<<<<<<<<<<<<<<<" << std::endl; + std::cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< " << imageMetadataInterface->GetSpectralSensitivity() << " <<<<<<<<<<<<<<<<<<<<<" << std::endl;*/ + + // Instantiating object + ReflectanceToSurfaceReflectanceImageFilter2Type::Pointer filter + = ReflectanceToSurfaceReflectanceImageFilter2Type::New(); + filter->SetInput(reader->GetOutput()); + filter->SetAtmoCorrectionParameters(paramAtmo); + writer->SetInput(filter->GetOutput()); + + writer->Update(); + + /*filter->GenerateParameters(); + otb::AtmosphericRadiativeTerms2::Pointer terms = filter->GetAtmosphericRadiativeTerms2(); + + std::ofstream fout (outputFileName); + fout << terms; + fout.close();*/ + + return EXIT_SUCCESS; +}