From 22280c7df731b77d6e4c0470b619395e69896667 Mon Sep 17 00:00:00 2001 From: Christophe Palmann <christophe.palmann@c-s.fr> Date: Wed, 19 Mar 2014 17:22:06 +0100 Subject: [PATCH] ENH:Improvement of optical calibration process support of user's parameters in optical calibration application add filters for luminance to image, reflectance to luminance, reflectance to image add related tests --- .../Radiometry/otbOpticalCalibration.cxx | 700 +++++++++++++----- .../otbLuminanceToImageImageFilter.h | 225 ++++++ .../otbReflectanceToImageImageFilter.h | 351 +++++++++ .../otbReflectanceToLuminanceImageFilter.h | 329 ++++++++ .../Applications/Radiometry/CMakeLists.txt | 34 + Testing/Code/Radiometry/CMakeLists.txt | 412 ++++++++++- .../otbLuminanceToImageImageFilter.cxx | 68 ++ .../otbLuminanceToImageImageFilterAuto.cxx | 74 ++ .../otbLuminanceToImageImageFilterNew.cxx | 37 + .../Code/Radiometry/otbRadiometryTests2.cxx | 9 + .../otbReflectanceToImageImageFilter.cxx | 183 +++++ .../otbReflectanceToImageImageFilterAuto.cxx | 72 ++ .../otbReflectanceToImageImageFilterNew.cxx | 37 + .../otbReflectanceToLuminanceImageFilter.cxx | 90 +++ ...bReflectanceToLuminanceImageFilterAuto.cxx | 58 ++ ...tbReflectanceToLuminanceImageFilterNew.cxx | 37 + 16 files changed, 2539 insertions(+), 177 deletions(-) create mode 100644 Code/Radiometry/otbLuminanceToImageImageFilter.h create mode 100644 Code/Radiometry/otbReflectanceToImageImageFilter.h create mode 100644 Code/Radiometry/otbReflectanceToLuminanceImageFilter.h create mode 100644 Testing/Code/Radiometry/otbLuminanceToImageImageFilter.cxx create mode 100644 Testing/Code/Radiometry/otbLuminanceToImageImageFilterAuto.cxx create mode 100644 Testing/Code/Radiometry/otbLuminanceToImageImageFilterNew.cxx create mode 100644 Testing/Code/Radiometry/otbReflectanceToImageImageFilter.cxx create mode 100644 Testing/Code/Radiometry/otbReflectanceToImageImageFilterAuto.cxx create mode 100644 Testing/Code/Radiometry/otbReflectanceToImageImageFilterNew.cxx create mode 100644 Testing/Code/Radiometry/otbReflectanceToLuminanceImageFilter.cxx create mode 100644 Testing/Code/Radiometry/otbReflectanceToLuminanceImageFilterAuto.cxx create mode 100644 Testing/Code/Radiometry/otbReflectanceToLuminanceImageFilterNew.cxx diff --git a/Applications/Radiometry/otbOpticalCalibration.cxx b/Applications/Radiometry/otbOpticalCalibration.cxx index 527040e11f..07bc1e45d9 100644 --- a/Applications/Radiometry/otbOpticalCalibration.cxx +++ b/Applications/Radiometry/otbOpticalCalibration.cxx @@ -20,6 +20,8 @@ #include "otbImageToLuminanceImageFilter.h" #include "otbLuminanceToReflectanceImageFilter.h" +#include "otbLuminanceToImageImageFilter.h" +#include "otbReflectanceToLuminanceImageFilter.h" #include "otbReflectanceToSurfaceReflectanceImageFilter.h" #include "otbMultiplyByScalarImageFilter.h" #include "otbClampVectorImageFilter.h" @@ -27,12 +29,20 @@ #include "otbGroundSpacingImageFunction.h" #include "vnl/vnl_random.h" +#include <fstream> +#include <sstream> +#include <vector> +#include <itkVariableLengthVector.h> + + + namespace otb { enum { - Level_TOA, + Level_IM_TOA, + Level_TOA_IM, Level_TOC }; @@ -63,14 +73,20 @@ public: itkTypeMacro(OpticalCalibration, Application); - typedef ImageToLuminanceImageFilter<UInt16VectorImageType, + typedef ImageToLuminanceImageFilter<DoubleVectorImageType, DoubleVectorImageType> ImageToLuminanceImageFilterType; typedef LuminanceToReflectanceImageFilter<DoubleVectorImageType, DoubleVectorImageType> LuminanceToReflectanceImageFilterType; + typedef LuminanceToImageImageFilter<DoubleVectorImageType, + DoubleVectorImageType> LuminanceToImageImageFilterType; + + typedef ReflectanceToLuminanceImageFilter<DoubleVectorImageType, + DoubleVectorImageType> ReflectanceToLuminanceImageFilterType; + typedef otb::MultiplyByScalarImageFilter<DoubleVectorImageType, - DoubleVectorImageType> ScaleFilterType; + DoubleVectorImageType> ScaleFilterOutDoubleType; typedef otb::ClampVectorImageFilter<DoubleVectorImageType, DoubleVectorImageType> ClampFilterType; @@ -85,23 +101,75 @@ public: typedef otb::SurfaceAdjacencyEffect6SCorrectionSchemeFilter<DoubleVectorImageType,DoubleVectorImageType> SurfaceAdjacencyEffect6SCorrectionSchemeFilterType; - typedef otb::GroundSpacingImageFunction<UInt16VectorImageType> GroundSpacingImageType; + typedef otb::GroundSpacingImageFunction<DoubleVectorImageType> GroundSpacingImageType; - typedef UInt16VectorImageType::IndexType IndexType; + typedef DoubleVectorImageType::IndexType IndexType; typedef GroundSpacingImageType::FloatType FloatType; typedef GroundSpacingImageType::ValueType ValueType; - typedef IndexType::IndexValueType IndexValueType; + typedef IndexType::IndexValueType IndexValueType; private: + + bool m_update1stTime; + string m_inImageName; + void DoInit() { SetName("OpticalCalibration"); - SetDescription("Perform optical calibration TOA/TOC (Top Of Atmosphere/Top Of Canopy). Supported sensors: QuickBird, Ikonos, WorldView2, Formosat, Spot5, Pleiades"); + SetDescription("Perform optical calibration TOA/TOC (Top Of Atmosphere/Top Of Canopy). Supported sensors: QuickBird, Ikonos, WorldView2, Formosat, Spot5, Pleiades, Spot6"); // Documentation SetDocName("Optical calibration"); - SetDocLongDescription("The application allows to convert pixel values from DN (for Digital Numbers) to physically interpretable and comparable values. Calibrated values are called surface reflectivity and its values lie in the range [0, 1].\nThe first level is called Top Of Atmosphere (TOA) reflectivity. It takes into account the sensor gain, sensor spectral response and the solar illumination.\nThe second level is called Top Of Canopy (TOC) reflectivity. In addition to sensor gain and solar illumination, it takes into account the optical thickness of the atmosphere, the atmospheric pressure, the water vapor amount, the ozone amount, as well as the composition and amount of aerosol gasses.\nIt is also possible to indicate an AERONET file which contains atmospheric parameters (version 1 and version 2 of Aeronet file are supported."); + SetDocLongDescription("The application allows to convert pixel values from DN (for Digital Numbers) to physically interpretable and comparable values. Calibrated values are called surface reflectivity and its values lie in the range [0, 1].\nThe first level is called Top Of Atmosphere (TOA) reflectivity. It takes into account the sensor gain, sensor spectral response and the solar illuminations.\nThe second level is called Top Of Canopy (TOC) reflectivity. In addition to sensor gain and solar illuminations, it takes into account the optical thickness of the atmosphere, the atmospheric pressure, the water vapor amount, the ozone amount, as well as the composition and amount of aerosol gasses.\nIt is also possible to indicate an AERONET file which contains atmospheric parameters (version 1 and version 2 of Aeronet file are supported.\n" +"\n--------------------------\n\n" +"If the sensor is not supported by the metadata interface factory of OTB, users still have the possibility to give the needed parameters to the application.\n" +"For TOA conversion, these parameters are : \n" +"- day and month of acquisition, or flux normalization coefficient;\n" +"- sun elevation angle;\n" +"- gains and biases, one pair of values for each band (passed by a file);\n" +"- solar illuminations, one value for each band (passed by a file).\n\n" +"For the conversion from DN (for Digital Numbers) to spectral radiance (or 'TOA radiance') L, the following formula is used :\n\n" + +"(1) L(b) = DN(b)/gain(b)+bias(b) (in W/m²/steradians/micrometers) with b being a band ID.\n\n" + +"These values are provided by the user thanks to a simple txt file with two lines, one for the gains and one for the biases.\n" +"Each value must be separated with colons (:), with eventual spaces. Blank lines are not allowed. If a line begins with the '#' symbol, then it is considered as comments.\n" +"Note that sometimes, the values provided by certain metadata files assume the formula L(b) = gain(b)*DC(b)+bias(b).\n" +"In this case, be sure to provide the inverse gain values so that the application can correctly interpret them.\n\n" + +"In order to convert TOA radiance to TOA reflectance, the following formula is used :\n\n" + +"(2) R(b) = (pi*L(b)*d²) / (ESUN(b)*cos(θ)) (no dimension) where : \n\n" + +"- L(b) is the spectral radiance for band b \n" +"- pi is the famous mathematical constant (3.14159...) \n" +"- d is the earth-sun distance (in astronomical units) and depends on the acquisition's day and month \n" +"- ESUN(b) is the mean TOA solar irradiance (or solar illumination) in W/m²/micrometers\n" +"- θ is the solar zenith angle in degrees. \n" +"Note that the application asks for the solar elevation angle, and will perfom the conversion to the zenith angle itself (ze. angle = 90° - el. angle).\n" +"Note also that ESUN(b) not only depends on the band b, but also on the spectral sensitivity of the sensor in this particular band. " +"In other words, the influence of spectral sensitivities is included within the ESUN different values.\n" +"These values are provided by the user thanks to a txt file following the same convention as before.\n" +"Instead of providing the date of acquisition, the user can also provide a flux normalization coefficient 'fn'. " +"The formula used instead will be the following : \n\n" +"(3) R(b) = (pi*L(b)) / (ESUN(b)*fn²*cos(θ)) \n\n" +"Whatever the formula used (2 or 3), the user should pay attention to the interpretation of the parameters he will provide to the application, " +"by taking into account the original formula that the metadata files assum.\n\n" + +"Below, we give two examples of txt files containing information about gains/biases and solar illuminations :\n\n" +"- gainbias.txt :\n" +"# Gain values for each band. Each value must be separated with colons (:), with eventual spaces. Blank lines not allowed.\n" +"10.4416 : 9.529 : 8.5175 : 14.0063\n" +"# Bias values for each band.\n" +"0.0 : 0.0 : 0.0 : 0.0\n\n" +"- solarillumination.txt : \n" +"# Solar illumination values in watt/m2/micron ('micron' means actually 'for each band').\n" +"# Each value must be separated with colons (:), with eventual spaces. Blank lines not allowed.\n" +"1540.494123 : 1826.087443 : 1982.671954 : 1094.747446\n\n" + +"Finally, the 'Logs' tab provides usefull messages that can help the user in knowing the process different status." +); SetDocLimitations("None"); SetDocAuthors("OTB-Team"); SetDocSeeAlso("The OTB CookBook"); @@ -111,19 +179,24 @@ private: AddParameter(ParameterType_InputImage, "in", "Input"); SetParameterDescription("in", "Input image filename (values in DN)"); + AddParameter(ParameterType_String, "sensor", "Sensor ID"); + SetParameterDescription("sensor", "Sensor ID"); + MandatoryOff("sensor"); + AddParameter(ParameterType_OutputImage, "out", "Output"); SetParameterDescription("out","Output calibrated image filename"); AddRAMParameter(); AddParameter(ParameterType_Choice, "level", "Calibration Level"); - AddChoice("level.toa", "TOA : Top Of Atmosphere"); + AddChoice("level.toa", "Image to TOA reflectance"); + AddChoice("level.toatoim", "TOA reflectance to Image"); AddChoice("level.toc", "TOC : Top Of Canopy (EXPERIMENTAL)"); SetParameterString("level", "toa"); AddParameter(ParameterType_Empty, "milli", "Convert to milli reflectance"); SetParameterDescription("milli", "Flag to use milli-reflectance instead of reflectance.\n" - "This allows to save the image with integer pixel type (in the range [0, 1000] instead of floating point in the range [0, 1]. In order to do that, use this option and set the output pixel type (-out filename uint16 for example)"); + "This allows to save the image with integer pixel type (in the range [0, 1000] instead of floating point in the range [0, 1]. In order to do that, use this option and set the output pixel type (-out filename double for example)"); DisableParameter("milli"); MandatoryOff("milli"); @@ -132,14 +205,43 @@ private: EnableParameter("clamp"); MandatoryOff("clamp"); - AddParameter(ParameterType_InputFilename, "rsr", "Relative Spectral Response File"); - std::ostringstream oss; - oss << "Sensor relative spectral response file"<<std::endl; - oss << "By default the application gets these informations in the metadata"; - SetParameterDescription("rsr", oss.str()); - MandatoryOff("rsr"); - - AddParameter(ParameterType_Group,"atmo","Atmospheric parameters"); + //Acquisition parameters (TOA) + AddParameter(ParameterType_Group,"acquisition","Acquisition parameters (TOA)"); + SetParameterDescription("acquisition","This group allows to set the parameters related to the acquisition conditions."); + //Day + AddParameter(ParameterType_Int, "acquisition.day", "Day"); + SetParameterDescription("acquisition.day", "Day (1-31)"); + SetMinimumParameterIntValue("acquisition.day", 1); + SetMaximumParameterIntValue("acquisition.day", 31); + MandatoryOff("acquisition.day"); + //Month + AddParameter(ParameterType_Int, "acquisition.month", "Month"); + SetParameterDescription("acquisition.month", "Month (1-12)"); + SetMinimumParameterIntValue("acquisition.month", 1); + SetMaximumParameterIntValue("acquisition.month", 12); + MandatoryOff("acquisition.month"); + //Flux normalization coefficient + AddParameter(ParameterType_Float, "acquisition.fluxnormalizationcoefficient", "Flux Normalization"); + SetParameterDescription("acquisition.fluxnormalizationcoefficient", "Flux Normalization Coefficient"); + SetMinimumParameterFloatValue("acquisition.fluxnormalizationcoefficient", 0.); + MandatoryOff("acquisition.fluxnormalizationcoefficient"); + //Sun elevation angle + AddParameter(ParameterType_Float, "acquisition.sunelevationangle", "Sun elevation angle (°)"); + SetParameterDescription("acquisition.sunelevationangle", "Sun elevation angle"); + SetMinimumParameterFloatValue("acquisition.sunelevationangle", 0.); + SetMaximumParameterFloatValue("acquisition.sunelevationangle", 120.); + MandatoryOff("acquisition.sunelevationangle"); + //Gain & bias + AddParameter(ParameterType_InputFilename, "acquisition.gainbias", "Gains | biases"); + SetParameterDescription("acquisition.gainbias", "Gains | biases"); + MandatoryOff("acquisition.gainbias"); + //Solar illuminations + AddParameter(ParameterType_InputFilename, "acquisition.solarilluminations", "Solar illuminations"); + SetParameterDescription("acquisition.solarilluminations", "Solar illuminations (one value per band)"); + MandatoryOff("acquisition.solarilluminations"); + + //Atmospheric parameters (TOC) + AddParameter(ParameterType_Group,"atmo","Atmospheric parameters (TOC)"); SetParameterDescription("atmo","This group allows to set the atmospheric parameters."); AddParameter(ParameterType_Choice, "atmo.aerosol", "Aerosol Model"); AddChoice("atmo.aerosol.noaersol", "No Aerosol Model"); @@ -173,208 +275,454 @@ private: AddParameter(ParameterType_InputFilename, "atmo.aeronet", "Aeronet File"); SetParameterDescription("atmo.aeronet","Aeronet file containing atmospheric parameters"); MandatoryOff("atmo.aeronet"); + + AddParameter(ParameterType_InputFilename, "atmo.rsr", "Relative Spectral Response File"); + std::ostringstream oss; + oss << "Sensor relative spectral response file"<<std::endl; + oss << "By default the application gets these informations in the metadata"; + SetParameterDescription("atmo.rsr", oss.str()); + MandatoryOff("atmo.rsr"); // Window radius for adjacency effects correction - AddParameter(ParameterType_Int, "radius", "Window radius"); - SetParameterDescription("radius","Window radius for adjacency effects corrections"); - MandatoryOff("radius"); - SetDefaultParameterInt("radius", 2); + AddParameter(ParameterType_Int, "atmo.radius", "Window radius (adjacency effects)"); + SetParameterDescription("atmo.radius","Window radius for adjacency effects corrections"); + MandatoryOff("atmo.radius"); + SetDefaultParameterInt("atmo.radius", 2); // Doc example parameter settings SetDocExampleParameterValue("in", "QB_1_ortho.tif"); SetDocExampleParameterValue("level", "toa"); SetDocExampleParameterValue("out", "OpticalCalibration.tif"); + + m_update1stTime = true; + m_inImageName = ""; } void DoUpdateParameters() { - // Nothing to update - } + + string tempName = GetParameterString("in"); + + if (!tempName.empty()) + { + + if (tempName != m_inImageName) + { + m_inImageName = tempName; + m_update1stTime = true; + } + + GetLogger()->Info("\nFile : " + m_inImageName + "\n"); + + //Check if valid metadata informations are available to compute ImageToLuminance and LuminanceToReflectance + DoubleVectorImageType::Pointer inImage = GetParameterDoubleVectorImage("in"); + itk::MetaDataDictionary dict = inImage->GetMetaDataDictionary(); + OpticalImageMetadataInterface::Pointer lImageMetadataInterface = OpticalImageMetadataInterfaceFactory::CreateIMI(dict); + + string IMIName( lImageMetadataInterface->GetNameOfClass() ) , IMIOptDfltName("OpticalDefaultImageMetadataInterface"); + if ( (IMIName != IMIOptDfltName) && (m_update1stTime) ) + { + itk::VariableLengthVector<double> vlvector; + std::stringstream ss; + + vlvector = lImageMetadataInterface->GetPhysicalGain(); + for(int k=0; k<vlvector.Size(); k++) + ss << vlvector[k] << " "; + ss << " | "; + vlvector = lImageMetadataInterface->GetPhysicalBias(); + for(int k=0; k<vlvector.Size(); k++) + ss << vlvector[k] << " "; + ss << " (custom values must be passed by a file)"; + SetParameterString("acquisition.gainbias",ss.str()); + + ss.str(std::string()); + vlvector = lImageMetadataInterface->GetSolarIrradiance(); + for(int k=0; k<vlvector.Size(); k++) + ss << vlvector[k] << " "; + ss << " (custom values must be passed by a file)"; + SetParameterString("acquisition.solarilluminations",ss.str()); + + + SetParameterInt("acquisition.day", lImageMetadataInterface->GetDay()); + SetParameterInt("acquisition.month", lImageMetadataInterface->GetMonth()); + SetParameterFloat("acquisition.sunelevationangle", lImageMetadataInterface->GetSunElevation()); + + SetParameterString("sensor", lImageMetadataInterface->GetSensorID()); + GetLogger()->Info( "\n-------------------------------------------------------------\n" + "Sensor ID : " + lImageMetadataInterface->GetSensorID() + "\n" + "\n-------------------------------------------------------------\n"); + + m_update1stTime=false; + + } + if ( (IMIName == IMIOptDfltName) && (m_update1stTime) ) + { + GetLogger()->Info("\n-------------------------------------------------------------\n" + "Sensor ID : unknown...\n" + "The application didn't manage to find an appropriate metadata interface; " + "custom values must be provided in order to perform TOA conversion.\nPlease, set the following fields :\n" + "- day and month of acquisition, or flux normalization coefficient;\n" + "- sun elevation angle;\n" + "- gains and biases for each band (passed by a file, see documentation);\n" + "- solar illuminationss for each band (passed by a file, see documentation).\n" + "-------------------------------------------------------------\n"); + + SetParameterString("sensor", "Unknown (see 'Logs' tabs)"); + SetParameterString("acquisition.gainbias",""); + SetParameterString("acquisition.solarilluminations",""); + SetParameterInt("acquisition.day", 1); + SetParameterInt("acquisition.month", 1); + SetParameterFloat("acquisition.fluxnormalizationcoefficient", 0); + SetParameterFloat("acquisition.sunelevationangle", 0); + + EnableParameter("acquisition.gainbias"); + EnableParameter("acquisition.solarilluminations"); + EnableParameter("acquisition.day"); + EnableParameter("acquisition.month"); + DisableParameter("acquisition.fluxnormalizationcoefficient"); + EnableParameter("acquisition.sunelevationangle"); + + + m_update1stTime=false; + } + } + + + } void DoExecute() { - UInt16VectorImageType::Pointer inImage = GetParameterUInt16VectorImage("in"); - - //Check if valid metadata informations are available to compute ImageToLuminance and LuminanceToReflectance - itk::MetaDataDictionary dict = inImage->GetMetaDataDictionary(); - OpticalImageMetadataInterface::Pointer lImageMetadataInterface = OpticalImageMetadataInterfaceFactory::CreateIMI(dict); - - // Test if needed data are available : an exception will be thrown - // if one the following Get* return failure. the exception is then - // caught in the Wrapper::Application class which redirect it to - // the logger - // ImageToLuminance - lImageMetadataInterface->GetPhysicalGain(); - lImageMetadataInterface->GetPhysicalBias(); - - // LuminanceToReflectance - lImageMetadataInterface->GetDay(); - lImageMetadataInterface->GetMonth(); - - lImageMetadataInterface->GetSolarIrradiance(); - lImageMetadataInterface->GetSunElevation(); - + //Main filters instanciations m_ImageToLuminanceFilter = ImageToLuminanceImageFilterType::New(); m_LuminanceToReflectanceFilter = LuminanceToReflectanceImageFilterType::New(); m_ReflectanceToSurfaceReflectanceFilter = ReflectanceToSurfaceReflectanceImageFilterType::New(); + m_ReflectanceToLuminanceFilter = ReflectanceToLuminanceImageFilterType::New(); + m_LuminanceToImageFilter = LuminanceToImageImageFilterType::New(); - m_ImageToLuminanceFilter->SetInput(inImage); - m_LuminanceToReflectanceFilter->SetInput(m_ImageToLuminanceFilter->GetOutput()); - m_ReflectanceToSurfaceReflectanceFilter->SetInput(m_LuminanceToReflectanceFilter->GetOutput()); - - m_ScaleFilter = ScaleFilterType::New(); + //Other instanciations + m_ScaleFilter = ScaleFilterOutDoubleType::New(); m_ScaleFilter->InPlaceOn(); - m_ClampFilter = ClampFilterType::New(); + //Check if valid metadata informations are available to compute ImageToLuminance and LuminanceToReflectance + DoubleVectorImageType::Pointer inImage = GetParameterDoubleVectorImage("in"); + itk::MetaDataDictionary dict = inImage->GetMetaDataDictionary(); + OpticalImageMetadataInterface::Pointer lImageMetadataInterface = OpticalImageMetadataInterfaceFactory::CreateIMI(dict); + + string IMIName( lImageMetadataInterface->GetNameOfClass() ) , IMIOptDfltName("OpticalDefaultImageMetadataInterface"); + if (IMIName != IMIOptDfltName) + { + // Test if needed data are available : an exception will be thrown + // if one the following Get* return failure. the exception is then + // caught in the Wrapper::Application class which redirect it to + // the logger + + // ImageToLuminance + lImageMetadataInterface->GetPhysicalGain(); + lImageMetadataInterface->GetPhysicalBias(); + + // LuminanceToReflectance + lImageMetadataInterface->GetDay(); + lImageMetadataInterface->GetMonth(); + + lImageMetadataInterface->GetSolarIrradiance(); + lImageMetadataInterface->GetSunElevation(); + lImageMetadataInterface->GetSolarIrradiance(); + } + else //No metadata interface + { + // Image to luminance / luminance to image + GetLogger()->Info("Trying to get gains/biases information...\n"); + string filename(GetParameterString("acquisition.gainbias")); + if (!filename.empty()) + { + std::ifstream file(filename.c_str(), std::ios::in); + + if(file) + { + + string line; unsigned int numLine = 0; + while (getline(file, line)) + { + if (line[0]!='#') + { + numLine++; + + std::vector<double> values; + std::istringstream iss(line); + string value; double dvalue; + while ( getline( iss, value, ':' ) ) + { + std::istringstream iss2(value); + iss2 >> dvalue; + values.push_back(dvalue); + } + + itk::VariableLengthVector<double> vlvector; + + vlvector.SetData(values.data(),values.size(),false); + + switch (numLine) + { + case 1 : m_ImageToLuminanceFilter->SetAlpha(vlvector); + m_LuminanceToImageFilter->SetAlpha(vlvector); + GetLogger()->Info("Trying to get gains/biases information... OK (1/2)\n"); + break; + case 2 : m_ImageToLuminanceFilter->SetBeta(vlvector); + m_LuminanceToImageFilter->SetBeta(vlvector); + GetLogger()->Info("Trying to get gains/biases information... OK (2/2)\n"); + break; + default : itkExceptionMacro(<< "File : " << filename << " contains wrong number of lines (needs two, one for gains and one for biases)"); + } + } + } + file.close(); + + } + else + itkExceptionMacro(<< "File : " << filename << " couldn't be opened"); + } + else + itkExceptionMacro(<< "Please, select a file containing gain/bias values for each band"); + + // Luminance to reflectance / reflectance to Luminance + m_LuminanceToReflectanceFilter->SetElevationSolarAngle(GetParameterFloat("acquisition.sunelevationangle")); + m_ReflectanceToLuminanceFilter->SetElevationSolarAngle(GetParameterFloat("acquisition.sunelevationangle")); + + if ( (IsParameterEnabled("acquisition.day")) && (IsParameterEnabled("acquisition.month")) ) + { + m_LuminanceToReflectanceFilter->SetDay(GetParameterInt("acquisition.day")); + m_LuminanceToReflectanceFilter->SetMonth(GetParameterInt("acquisition.month")); + + m_ReflectanceToLuminanceFilter->SetDay(GetParameterInt("acquisition.day")); + m_ReflectanceToLuminanceFilter->SetMonth(GetParameterInt("acquisition.month")); + } + else if (IsParameterEnabled("acquisition.fluxnormalizationcoefficient")) + { + m_LuminanceToReflectanceFilter->SetFluxNormalizationCoefficient(GetParameterFloat("acquisition.fluxnormalizationcoefficient")); + + m_ReflectanceToLuminanceFilter->SetFluxNormalizationCoefficient(GetParameterFloat("acquisition.fluxnormalizationcoefficient")); + } + else + itkExceptionMacro(<< "Please, set the day and month fields, OR set the flux normalization coefficient field"); + + GetLogger()->Info("Trying to get solar illuminations information...\n"); + string filename2(GetParameterString("acquisition.solarilluminations")); + if(!filename2.empty()) + { + std::ifstream file2(filename2.c_str(), std::ios::in); + if(file2) + { + string line; + while (getline(file2, line)) + { + if (line[0]!='#') + { + std::vector<double> values; + std::istringstream iss(line); + string value; double dvalue; + while ( getline( iss, value, ':' ) ) + { + std::istringstream iss2(value); + iss2 >> dvalue; + + values.push_back(dvalue); + } + + itk::VariableLengthVector<double> vlvector; + vlvector.SetData(values.data(),values.size(),false); + + m_LuminanceToReflectanceFilter->SetSolarIllumination(vlvector); + m_ReflectanceToLuminanceFilter->SetSolarIllumination(vlvector); + + GetLogger()->Info("Trying to get solar illuminations information... OK\n"); + } + } + file2.close(); + } + else + itkExceptionMacro(<< "File : " << filename2 << " couldn't be opened"); + } + else + itkExceptionMacro(<< "Please, select a file containing solar illumination values for each band"); + + } + switch ( GetParameterInt("level") ) { - case Level_TOA: + case Level_IM_TOA: { - GetLogger()->Info("Compute Top of Atmosphere reflectance\n"); - m_LuminanceToReflectanceFilter->SetUseClamp(IsParameterEnabled("clamp")); + GetLogger()->Info("Compute Top of Atmosphere reflectance\n"); + + //Pipeline + m_ImageToLuminanceFilter->SetInput(inImage); + m_LuminanceToReflectanceFilter->SetInput(m_ImageToLuminanceFilter->GetOutput()); + m_LuminanceToReflectanceFilter->SetUseClamp(IsParameterEnabled("clamp")); + m_LuminanceToReflectanceFilter->UpdateOutputInformation(); + m_ScaleFilter->SetInput(m_LuminanceToReflectanceFilter->GetOutput()); + } + break; + case Level_TOA_IM: + { + + GetLogger()->Info("Convert Top of Atmosphere reflectance to image DN\n"); + + //Pipeline + m_ReflectanceToLuminanceFilter->SetInput(inImage); + m_LuminanceToImageFilter->SetInput(m_ReflectanceToLuminanceFilter->GetOutput()); + m_LuminanceToImageFilter->UpdateOutputInformation(); + m_ScaleFilter->SetInput(m_LuminanceToImageFilter->GetOutput()); - m_LuminanceToReflectanceFilter->UpdateOutputInformation(); - m_ScaleFilter->SetInput(m_LuminanceToReflectanceFilter->GetOutput()); } break; case Level_TOC: { - GetLogger()->Info("Compute Top of Canopy reflectance\n"); - m_ReflectanceToSurfaceReflectanceFilter->SetIsSetAtmosphericRadiativeTerms(false); - m_ReflectanceToSurfaceReflectanceFilter->SetUseGenerateParameters(true); - m_ReflectanceToSurfaceReflectanceFilter->UpdateOutputInformation(); - m_ReflectanceToSurfaceReflectanceFilter->SetUseGenerateParameters(false); - - m_AtmosphericParam = m_ReflectanceToSurfaceReflectanceFilter->GetCorrectionParameters(); - //AerosolModelType aeroMod = AtmosphericCorrectionParametersType::NO_AEROSOL; - - switch ( GetParameterInt("atmo.aerosol") ) - { - case Aerosol_Desertic: - { - // Aerosol_Desertic correspond to 4 in the enum but actually in - // the class atmosphericParam it is known as parameter 5 - m_AtmosphericParam->SetAerosolModel(static_cast<AerosolModelType>(5)); - } - break; - default: - { - m_AtmosphericParam->SetAerosolModel(static_cast<AerosolModelType>(GetParameterInt("atmo.aerosol"))); - } - break; - } - // Set the atmospheric param - m_AtmosphericParam->SetOzoneAmount(GetParameterFloat("atmo.oz")); - m_AtmosphericParam->SetWaterVaporAmount(GetParameterFloat("atmo.wa")); - m_AtmosphericParam->SetAtmosphericPressure(GetParameterFloat("atmo.pressure")); - m_AtmosphericParam->SetAerosolOptical(GetParameterFloat("atmo.opt")); - - // Relative Spectral Response File - if (IsParameterEnabled("rsr")) - { - m_ReflectanceToSurfaceReflectanceFilter->SetFilterFunctionValuesFileName(GetParameterString("rsr")); - } - - // Aeronet file - if (IsParameterEnabled("atmo.aeronet")) - { - GetLogger()->Info("Use aeronet file to retrieve atmospheric parameters"); - m_ReflectanceToSurfaceReflectanceFilter->SetAeronetFileName(GetParameterString("atmo.aeronet")); - } - - m_ReflectanceToSurfaceReflectanceFilter->SetIsSetAtmosphericRadiativeTerms(false); - m_ReflectanceToSurfaceReflectanceFilter->SetUseGenerateParameters(true); - m_ReflectanceToSurfaceReflectanceFilter->GenerateParameters(); - m_ReflectanceToSurfaceReflectanceFilter->UpdateOutputInformation(); - m_ReflectanceToSurfaceReflectanceFilter->SetUseGenerateParameters(false); - - // std::ostringstream oss_atmo; - // oss_atmo << "Atmospheric parameters: " << std::endl; - // oss_atmo << m_AtmosphericParam; - - // GetLogger()->Info(oss_atmo.str()); - - std::ostringstream oss; - oss.str(""); - oss << std::endl << m_AtmosphericParam; - - AtmosphericRadiativeTerms::Pointer atmoTerms = m_ReflectanceToSurfaceReflectanceFilter->GetAtmosphericRadiativeTerms(); - oss << std::endl << std::endl << atmoTerms; - - GetLogger()->Info("Atmospheric correction parameters compute by 6S : " + oss.str()); - - //Compute adjacency effect - // m_SurfaceAdjacencyEffect6SCorrectionSchemeFilter - // = SurfaceAdjacencyEffect6SCorrectionSchemeFilterType::New(); - - // m_SurfaceAdjacencyEffect6SCorrectionSchemeFilter->SetInput(m_ReflectanceToSurfaceReflectanceFilter->GetOutput()); - // m_SurfaceAdjacencyEffect6SCorrectionSchemeFilter-> - // SetAtmosphericRadiativeTerms( - // m_ReflectanceToSurfaceReflectanceFilter->GetAtmosphericRadiativeTerms()); - // m_SurfaceAdjacencyEffect6SCorrectionSchemeFilter->SetZenithalViewingAngle( - // m_AtmosphericParam->GetViewingZenithalAngle()); - // m_SurfaceAdjacencyEffect6SCorrectionSchemeFilter->SetWindowRadius(GetParameterInt("radius")); - - // //estimate ground spacing in kilometers - // GroundSpacingImageType::Pointer groundSpacing = GroundSpacingImageType::New(); - - // groundSpacing->SetInputImage(inImage); - // IndexType index; - - // vnl_random rand; - - // index[0] = static_cast<IndexValueType>(rand.lrand32(0, inImage->GetLargestPossibleRegion().GetSize()[0])); - // index[1] = static_cast<IndexValueType>(rand.lrand32(0, inImage->GetLargestPossibleRegion().GetSize()[1])); - // FloatType tmpSpacing = groundSpacing->EvaluateAtIndex(index); - - // const float spacingInKilometers = (std::max(tmpSpacing[0], tmpSpacing[1])) / 1000.; - - // // std::ostringstream oss2; - // // oss2.str(""); - // // oss2 << spacingInKilometers; - - // // GetLogger()->Info("Spacing in kilometers " + oss2.str()); - // m_SurfaceAdjacencyEffect6SCorrectionSchemeFilter-> - // SetPixelSpacingInKilometers(spacingInKilometers); - - // //rescale the surface reflectance in milli-reflectance - // m_SurfaceAdjacencyEffect6SCorrectionSchemeFilter->UpdateOutputInformation(); - // //m_SurfaceAdjacencyEffect6SCorrectionSchemeFilter->Update(); - // m_ScaleFilter->SetInput(m_SurfaceAdjacencyEffect6SCorrectionSchemeFilter->GetOutput()); - if (!IsParameterEnabled("clamp")) - { - m_ScaleFilter->SetInput(m_ReflectanceToSurfaceReflectanceFilter->GetOutput()); - } - else - { - GetLogger()->Info("Clamp values between [0, 100]"); - m_ClampFilter->SetInput(m_ReflectanceToSurfaceReflectanceFilter->GetOutput()); - m_ClampFilter->ClampOutside(0.0, 1.0); - m_ScaleFilter->SetInput(m_ClampFilter->GetOutput()); - } + GetLogger()->Info("Convert Top of Canopy reflectance\n"); + + //Pipeline + m_ImageToLuminanceFilter->SetInput(inImage); + m_LuminanceToReflectanceFilter->SetInput(m_ImageToLuminanceFilter->GetOutput()); + m_ReflectanceToSurfaceReflectanceFilter->SetInput(m_LuminanceToReflectanceFilter->GetOutput()); + + m_ReflectanceToSurfaceReflectanceFilter->SetIsSetAtmosphericRadiativeTerms(false); + m_ReflectanceToSurfaceReflectanceFilter->SetUseGenerateParameters(true); + m_ReflectanceToSurfaceReflectanceFilter->UpdateOutputInformation(); + m_ReflectanceToSurfaceReflectanceFilter->SetUseGenerateParameters(false); + + m_AtmosphericParam = m_ReflectanceToSurfaceReflectanceFilter->GetCorrectionParameters(); + //AerosolModelType aeroMod = AtmosphericCorrectionParametersType::NO_AEROSOL; + + switch ( GetParameterInt("atmo.aerosol") ) + { + case Aerosol_Desertic: + { + // Aerosol_Desertic correspond to 4 in the enum but actually in + // the class atmosphericParam it is known as parameter 5 + m_AtmosphericParam->SetAerosolModel(static_cast<AerosolModelType>(5)); + } + break; + default: + { + m_AtmosphericParam->SetAerosolModel(static_cast<AerosolModelType>(GetParameterInt("atmo.aerosol"))); + } + break; + } + // Set the atmospheric param + m_AtmosphericParam->SetOzoneAmount(GetParameterFloat("atmo.oz")); + m_AtmosphericParam->SetWaterVaporAmount(GetParameterFloat("atmo.wa")); + m_AtmosphericParam->SetAtmosphericPressure(GetParameterFloat("atmo.pressure")); + m_AtmosphericParam->SetAerosolOptical(GetParameterFloat("atmo.opt")); + + // Relative Spectral Response File + if (IsParameterEnabled("rsr")) + { + m_ReflectanceToSurfaceReflectanceFilter->SetFilterFunctionValuesFileName(GetParameterString("rsr")); + } + + // Aeronet file + if (IsParameterEnabled("atmo.aeronet")) + { + GetLogger()->Info("Use aeronet file to retrieve atmospheric parameters"); + m_ReflectanceToSurfaceReflectanceFilter->SetAeronetFileName(GetParameterString("atmo.aeronet")); + } + + m_ReflectanceToSurfaceReflectanceFilter->SetIsSetAtmosphericRadiativeTerms(false); + m_ReflectanceToSurfaceReflectanceFilter->SetUseGenerateParameters(true); + m_ReflectanceToSurfaceReflectanceFilter->GenerateParameters(); + m_ReflectanceToSurfaceReflectanceFilter->UpdateOutputInformation(); + m_ReflectanceToSurfaceReflectanceFilter->SetUseGenerateParameters(false); + + // std::ostringstream oss_atmo; + // oss_atmo << "Atmospheric parameters: " << std::endl; + // oss_atmo << m_AtmosphericParam; + + // GetLogger()->Info(oss_atmo.str()); + + std::ostringstream oss; + oss.str(""); + oss << std::endl << m_AtmosphericParam; + + AtmosphericRadiativeTerms::Pointer atmoTerms = m_ReflectanceToSurfaceReflectanceFilter->GetAtmosphericRadiativeTerms(); + oss << std::endl << std::endl << atmoTerms; + + GetLogger()->Info("Atmospheric correction parameters compute by 6S : " + oss.str()); + + //Compute adjacency effect + // m_SurfaceAdjacencyEffect6SCorrectionSchemeFilter + // = SurfaceAdjacencyEffect6SCorrectionSchemeFilterType::New(); + + // m_SurfaceAdjacencyEffect6SCorrectionSchemeFilter->SetInput(m_ReflectanceToSurfaceReflectanceFilter->GetOutput()); + // m_SurfaceAdjacencyEffect6SCorrectionSchemeFilter-> + // SetAtmosphericRadiativeTerms( + // m_ReflectanceToSurfaceReflectanceFilter->GetAtmosphericRadiativeTerms()); + // m_SurfaceAdjacencyEffect6SCorrectionSchemeFilter->SetZenithalViewingAngle( + // m_AtmosphericParam->GetViewingZenithalAngle()); + // m_SurfaceAdjacencyEffect6SCorrectionSchemeFilter->SetWindowRadius(GetParameterInt("radius")); + + // //estimate ground spacing in kilometers + // GroundSpacingImageType::Pointer groundSpacing = GroundSpacingImageType::New(); + + // groundSpacing->SetInputImage(inImage); + // IndexType index; + + // vnl_random rand; + + // index[0] = static_cast<IndexValueType>(rand.lrand32(0, inImage->GetLargestPossibleRegion().GetSize()[0])); + // index[1] = static_cast<IndexValueType>(rand.lrand32(0, inImage->GetLargestPossibleRegion().GetSize()[1])); + // FloatType tmpSpacing = groundSpacing->EvaluateAtIndex(index); + + // const float spacingInKilometers = (std::max(tmpSpacing[0], tmpSpacing[1])) / 1000.; + + // // std::ostringstream oss2; + // // oss2.str(""); + // // oss2 << spacingInKilometers; + + // // GetLogger()->Info("Spacing in kilometers " + oss2.str()); + // m_SurfaceAdjacencyEffect6SCorrectionSchemeFilter-> + // SetPixelSpacingInKilometers(spacingInKilometers); + + // //rescale the surface reflectance in milli-reflectance + // m_SurfaceAdjacencyEffect6SCorrectionSchemeFilter->UpdateOutputInformation(); + // //m_SurfaceAdjacencyEffect6SCorrectionSchemeFilter->Update(); + // m_ScaleFilter->SetInput(m_SurfaceAdjacencyEffect6SCorrectionSchemeFilter->GetOutput()); + if (!IsParameterEnabled("clamp")) + { + m_ScaleFilter->SetInput(m_ReflectanceToSurfaceReflectanceFilter->GetOutput()); + } + else + { + GetLogger()->Info("Clamp values between [0, 100]"); + m_ClampFilter->SetInput(m_ReflectanceToSurfaceReflectanceFilter->GetOutput()); + m_ClampFilter->ClampOutside(0.0, 1.0); + m_ScaleFilter->SetInput(m_ClampFilter->GetOutput()); + } } break; - } + } // Output Image double scale = 1.0; - if (IsParameterEnabled("milli")) { - GetLogger()->Info("Use milli-reflectance"); - scale = 1000.; + GetLogger()->Info("Use milli-reflectance"); + if ( (GetParameterInt("level") == Level_IM_TOA) || (GetParameterInt("level") == Level_TOC) ) + scale = 1000.; + if (GetParameterInt("level") == Level_TOA_IM) + scale = 1. / 1000.; } + m_ScaleFilter->SetCoef(scale); - m_ScaleFilter->SetCoef(scale); + SetParameterOutputImage("out", m_ScaleFilter->GetOutput()); - SetParameterOutputImage("out", m_ScaleFilter->GetOutput()); } + //Keep object references as a members of the class, else the pipeline will be broken after exiting DoExecute(). ImageToLuminanceImageFilterType ::Pointer m_ImageToLuminanceFilter; LuminanceToReflectanceImageFilterType::Pointer m_LuminanceToReflectanceFilter; + ReflectanceToLuminanceImageFilterType::Pointer m_ReflectanceToLuminanceFilter; + LuminanceToImageImageFilterType::Pointer m_LuminanceToImageFilter; ReflectanceToSurfaceReflectanceImageFilterType::Pointer m_ReflectanceToSurfaceReflectanceFilter; - ScaleFilterType::Pointer m_ScaleFilter; + ScaleFilterOutDoubleType::Pointer m_ScaleFilter; AtmosphericCorrectionParametersType::Pointer m_AtmosphericParam; ClampFilterType::Pointer m_ClampFilter; diff --git a/Code/Radiometry/otbLuminanceToImageImageFilter.h b/Code/Radiometry/otbLuminanceToImageImageFilter.h new file mode 100644 index 0000000000..b3c660c6b8 --- /dev/null +++ b/Code/Radiometry/otbLuminanceToImageImageFilter.h @@ -0,0 +1,225 @@ +/*========================================================================= + + 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 __otbLuminanceToImageImageFilter_h +#define __otbLuminanceToImageImageFilter_h + +#include "otb_6S.h" +#include "otbUnaryImageFunctorWithVectorImageFilter.h" +#include "otbVectorImage.h" +#include "itkNumericTraits.h" +#include "itkVariableLengthVector.h" +#include "otbMacro.h" +#include "otbOpticalImageMetadataInterface.h" +#include "otbOpticalImageMetadataInterfaceFactory.h" + +#include <fstream> + +namespace otb +{ +namespace Functor +{ +/** + * \class LuminanceToImageImageFunctor + * \brief Substract beta to the Input and multiply by alpha. + * + * \sa LuminanceToImageImageFilter + * \ingroup Functor + * \ingroup Radiometry + */ + +template <class TInput, class TOutput> +class LuminanceToImageImageFunctor +{ +public: + LuminanceToImageImageFunctor() : + m_Alpha(1.), + m_Beta(0.) + {} + + virtual ~LuminanceToImageImageFunctor() {} + + void SetAlpha(double alpha) + { + m_Alpha = alpha; + } + void SetBeta(double beta) + { + m_Beta = beta; + } + double GetAlpha() + { + return m_Alpha; + } + double GetBeta() + { + return m_Beta; + } + + inline TOutput operator ()(const TInput& inPixel) const + { + TOutput outPixel; + double temp; + temp = (static_cast<double>(inPixel) - m_Beta) * m_Alpha; + outPixel = static_cast<TOutput>(temp); + return outPixel; + } + +private: + double m_Alpha; + double m_Beta; +}; +} + +/** \class LuminanceToImageImageFilter + * \brief Convert a raw value into a luminance value + * + * Transform a luminance image into a classical image. For this it + * uses the functor LuminanceToImageImageFunctor calling for each component of each pixel. + * + * + * For Spot image in the dimap format, the correction parameters are + * retrieved automatically from the metadata + * + * \ingroup LuminanceToImageImageFunctor + * \ingroup Radiometry + * + * \example Radiometry/AtmosphericCorrectionSequencement.cxx + */ +template <class TInputImage, class TOutputImage> +class ITK_EXPORT LuminanceToImageImageFilter : + public UnaryImageFunctorWithVectorImageFilter<TInputImage, + TOutputImage, + typename Functor::LuminanceToImageImageFunctor<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::LuminanceToImageImageFunctor<typename InputImageType::InternalPixelType, + typename OutputImageType::InternalPixelType> FunctorType; + + /** "typedef" for standard classes. */ + typedef LuminanceToImageImageFilter 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(LuminanceToImageImageFilter, UnaryImageFunctorWithVectorImageFiltermageFilter); + + /** 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 typename itk::VariableLengthVector<double> VectorType; + + /** Image size "typedef" definition. */ + typedef typename InputImageType::SizeType SizeType; + + /** Set the absolute calibration gains. */ + itkSetMacro(Alpha, VectorType); + + /** Give the absolute calibration gains. */ + itkGetConstReferenceMacro(Alpha, VectorType); + + /** Set the absolute calibration bias. */ + itkSetMacro(Beta, VectorType); + /** Give the absolute calibration bias. */ + itkGetConstReferenceMacro(Beta, VectorType); + +protected: + /** Constructor */ + LuminanceToImageImageFilter() + { + m_Alpha.SetSize(0); + m_Beta.SetSize(0); + }; + + /** Destructor */ + virtual ~LuminanceToImageImageFilter() {} + + /** Update the functor list and input parameters */ + virtual void BeforeThreadedGenerateData(void) + { + OpticalImageMetadataInterface::Pointer imageMetadataInterface = OpticalImageMetadataInterfaceFactory::CreateIMI( + this->GetInput()->GetMetaDataDictionary()); + if (m_Alpha.GetSize() == 0) + { + m_Alpha = imageMetadataInterface->GetPhysicalGain(); + } + + if (m_Beta.GetSize() == 0) + { + m_Beta = imageMetadataInterface->GetPhysicalBias(); + } + + otbMsgDevMacro(<< "Dimension: "); + otbMsgDevMacro(<< "m_Alpha.GetSize(): " << m_Alpha.GetSize()); + otbMsgDevMacro(<< "m_Beta.GetSize() : " << m_Beta.GetSize()); + otbMsgDevMacro( + << "this->GetInput()->GetNumberOfComponentsPerPixel() : " << this->GetInput()->GetNumberOfComponentsPerPixel()); + + if ((m_Alpha.GetSize() != this->GetInput()->GetNumberOfComponentsPerPixel()) + || (m_Beta.GetSize() != this->GetInput()->GetNumberOfComponentsPerPixel())) + { + itkExceptionMacro(<< "Alpha and Beta parameters should have the same size as the number of bands"); + } + + otbMsgDevMacro(<< "Using correction parameters: "); + otbMsgDevMacro(<< "Alpha (gain): " << m_Alpha); + otbMsgDevMacro(<< "Beta (bias): " << m_Beta); + + this->GetFunctorVector().clear(); + for (unsigned int i = 0; i < this->GetInput()->GetNumberOfComponentsPerPixel(); ++i) + { + FunctorType functor; + functor.SetAlpha(m_Alpha[i]); + functor.SetBeta(m_Beta[i]); + this->GetFunctorVector().push_back(functor); + } + } + +private: + /** Ponderation declaration*/ + VectorType m_Alpha; + VectorType m_Beta; +}; + +} // end namespace otb + +#endif diff --git a/Code/Radiometry/otbReflectanceToImageImageFilter.h b/Code/Radiometry/otbReflectanceToImageImageFilter.h new file mode 100644 index 0000000000..8fe985bf64 --- /dev/null +++ b/Code/Radiometry/otbReflectanceToImageImageFilter.h @@ -0,0 +1,351 @@ +/*========================================================================= + + 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 __otbReflectanceToImageImageFilter_h +#define __otbReflectanceToImageImageFilter_h + +#include "otbUnaryImageFunctorWithVectorImageFilter.h" +#include "otbLuminanceToImageImageFilter.h" +#include "otbReflectanceToLuminanceImageFilter.h" +#include "otbVectorImage.h" +#include "itkNumericTraits.h" +#include "itkVariableLengthVector.h" + +namespace otb +{ +namespace Functor +{ +/** \class ReflectanceToImageImageFunctor + * \brief Call the ReflectanceToLuminanceFunctor over the input and the LuminanceToImageFunctor to this result. + * + * + * \sa ReflectanceToImageImageFilter + * + * \ingroup Functor + * \ingroup LuminanceToImageFunctor + * \ingroup ReflectanceToLuminanceFunctor + * \ingroup Radiometry + */ +template <class TInput, class TOutput> +class ReflectanceToImageImageFunctor +{ +public: + ReflectanceToImageImageFunctor() {} + virtual ~ReflectanceToImageImageFunctor() {} + + typedef Functor::LuminanceToImageImageFunctor<TInput, TOutput> LumToImFunctorType; + typedef Functor::ReflectanceToLuminanceImageFunctor<TInput, TOutput> ReflecToLumFunctorType; + + void SetAlpha(double alpha) + { + m_LumToImFunctor.SetAlpha(alpha); + } + void SetBeta(double beta) + { + m_LumToImFunctor.SetBeta(beta); + } + void SetSolarIllumination(double solarIllumination) + { + m_ReflecToLumFunctor.SetSolarIllumination(solarIllumination); + } + void SetIlluminationCorrectionCoefficient(double coef) + { + m_ReflecToLumFunctor.SetIlluminationCorrectionCoefficient(coef); + } + + double GetAlpha() + { + return m_LumToImFunctor.GetAlpha(); + } + double GetBeta() + { + return m_LumToImFunctor.GetBeta(); + } + double GetSolarIllumination() + { + return m_ReflecToLumFunctor.GetSolarIllumination(); + } + double GetIlluminationCorrectionCoefficient() + { + return m_ReflecToLumFunctor.GetIlluminationCorrectionCoefficient(); + } + + inline TOutput operator ()(const TInput& inPixel) const + { + TOutput outPixel; + TOutput tempPix; + tempPix = m_ReflecToLumFunctor(inPixel); + outPixel = m_LumToImFunctor(tempPix); + + return outPixel; + } + +private: + LumToImFunctorType m_LumToImFunctor; + ReflecToLumFunctorType m_ReflecToLumFunctor; + +}; +} + +/** \class ReflectanceToImageImageFilter + * \brief Convert a reflectance into a raw value value + * + * Transform a reflectance image into a classical image image. For this it uses the functor ReflectanceToImageFunctor calling for each component of each pixel. + * The flux normalization coefficient (that is the ratio solar distance over mean solar distance) can be directly set or the user can + * give the day and the mounth of the observation and the class will used a coefficient given by a 6S routine that will give the corresponding coefficient. + * To note that in the case, 6S gives the square of the distances ratio. + * + * + * For Spot image in the dimap format, the correction parameters are + * retrieved automatically from the metadata + * + * \ingroup ReflectanceToImageImageFunctor + * \ingroup LuminanceToImageImageFilter + * \ingroup ReflectanceToLuminanceImageFilter + * \ingroup Radiometry + */ +template <class TInputImage, class TOutputImage> +class ITK_EXPORT ReflectanceToImageImageFilter : + public UnaryImageFunctorWithVectorImageFilter<TInputImage, + TOutputImage, + typename Functor::ReflectanceToImageImageFunctor<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::ReflectanceToImageImageFunctor<typename InputImageType::InternalPixelType, + typename OutputImageType::InternalPixelType> FunctorType; + + /** "typedef" for standard classes. */ + typedef ReflectanceToImageImageFilter 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(ReflectanceToImageImageFilter, UnaryImageFunctorWithVectorImageFiltermageFilter); + + /** 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 typename itk::VariableLengthVector<double> VectorType; + + /** Image size "typedef" definition. */ + typedef typename InputImageType::SizeType SizeType; + + /** Set the absolute calibration gains. */ + itkSetMacro(Alpha, VectorType); + /** Give the absolute calibration gains. */ + itkGetConstReferenceMacro(Alpha, VectorType); + + /** Set the absolute calibration bias. */ + itkSetMacro(Beta, VectorType); + /** Give the absolute calibration bias. */ + itkGetConstReferenceMacro(Beta, VectorType); + /** Set the solar illumination value. */ + itkSetMacro(SolarIllumination, VectorType); + + /** Give the solar illumination value. */ + itkGetConstReferenceMacro(SolarIllumination, VectorType); + + /** Set the zenithal solar angle. */ + itkSetMacro(ZenithalSolarAngle, double); + /** Give the zenithal solar angle. */ + itkGetConstReferenceMacro(ZenithalSolarAngle, double); + + /** Set/Get the sun elevation angle (internally handled by the zenithal angle)*/ + virtual void SetElevationSolarAngle(double elevationAngle) + { + double zenithalAngle = 90.0 - elevationAngle; + if (this->m_ZenithalSolarAngle != zenithalAngle) + { + this->m_ZenithalSolarAngle = zenithalAngle; + this->Modified(); + } + } + + virtual double GetElevationSolarAngle() const + { + return 90.0 - this->m_ZenithalSolarAngle; + } + + /** Set the flux normalization coefficient. */ + void SetFluxNormalizationCoefficient(double coef) + { + m_FluxNormalizationCoefficient = coef; + m_IsSetFluxNormalizationCoefficient = true; + this->Modified(); + } + + /** Set the acquisition day. */ + itkSetClampMacro(Day, int, 1, 31); + /** Get the acquisition day. */ + itkGetConstReferenceMacro(Day, int); + /** Set the acquisition mounth. */ + itkSetClampMacro(Month, int, 1, 12); + /** Set the acquisition mounth. */ + itkGetConstReferenceMacro(Month, int); + +protected: + /** Constructor */ + ReflectanceToImageImageFilter() : + m_ZenithalSolarAngle(120.), //invalid value which will lead to negative radiometry + m_FluxNormalizationCoefficient(1.), + m_IsSetFluxNormalizationCoefficient(false), + m_Day(0), + m_Month(0) + { + m_Alpha.SetSize(0); + m_Beta.SetSize(0); + m_SolarIllumination.SetSize(0); + }; + + /** Destructor */ + virtual ~ReflectanceToImageImageFilter() {} + + /** Update the functor list and input parameters */ + virtual void BeforeThreadedGenerateData(void) + { + + OpticalImageMetadataInterface::Pointer imageMetadataInterface = OpticalImageMetadataInterfaceFactory::CreateIMI( + this->GetInput()->GetMetaDataDictionary()); + if (m_Alpha.GetSize() == 0) + { + m_Alpha = imageMetadataInterface->GetPhysicalGain(); + } + + if (m_Beta.GetSize() == 0) + { + m_Beta = imageMetadataInterface->GetPhysicalBias(); + } + + if ((m_Day == 0) && (!m_IsSetFluxNormalizationCoefficient)) + { + m_Day = imageMetadataInterface->GetDay(); + } + + if ((m_Month == 0) && (!m_IsSetFluxNormalizationCoefficient)) + { + m_Month = imageMetadataInterface->GetMonth(); + } + + if (m_SolarIllumination.GetSize() == 0) + { + m_SolarIllumination = imageMetadataInterface->GetSolarIrradiance(); + } + + if (m_ZenithalSolarAngle == 120.0) + { + //the zenithal angle is the complementary of the elevation angle + m_ZenithalSolarAngle = 90.0 - imageMetadataInterface->GetSunElevation(); + } + + otbMsgDevMacro(<< "Using correction parameters: "); + otbMsgDevMacro(<< "Alpha (gain): " << m_Alpha); + otbMsgDevMacro(<< "Beta (bias): " << m_Beta); + otbMsgDevMacro(<< "Day: " << m_Day); + otbMsgDevMacro(<< "Month: " << m_Month); + otbMsgDevMacro(<< "Solar irradiance: " << m_SolarIllumination); + otbMsgDevMacro(<< "Zenithal angle: " << m_ZenithalSolarAngle); + + if ((m_Alpha.GetSize() != this->GetInput()->GetNumberOfComponentsPerPixel()) + || (m_Beta.GetSize() != this->GetInput()->GetNumberOfComponentsPerPixel()) + || (m_SolarIllumination.GetSize() != this->GetInput()->GetNumberOfComponentsPerPixel())) + { + itkExceptionMacro( + << "Alpha, Beta and SolarIllumination parameters should have the same size as the number of bands"); + } + + this->GetFunctorVector().clear(); + for (unsigned int i = 0; i < this->GetInput()->GetNumberOfComponentsPerPixel(); ++i) + { + FunctorType functor; + double coefTemp = 0.; + if (!m_IsSetFluxNormalizationCoefficient) + { + if (m_Day * m_Month != 0 && m_Day < 32 && m_Month < 13) + { + otb_6s_doublereal dsol = 0.; + otb_6s_integer day = static_cast<otb_6s_integer>(m_Day); + otb_6s_integer month = static_cast<otb_6s_integer>(m_Month); + //int cr(0); + otb_6s_varsol_(&day, &month, &dsol); + coefTemp = vcl_cos(m_ZenithalSolarAngle * CONST_PI_180) * static_cast<double>(dsol); + } + else + { + itkExceptionMacro(<< "Day has to be included between 1 and 31, Month beetween 1 and 12."); + } + } + else + { + coefTemp = + vcl_cos(m_ZenithalSolarAngle * + CONST_PI_180) * m_FluxNormalizationCoefficient * m_FluxNormalizationCoefficient; + } + functor.SetIlluminationCorrectionCoefficient(coefTemp); + functor.SetAlpha(m_Alpha[i]); + functor.SetBeta(m_Beta[i]); + functor.SetSolarIllumination(m_SolarIllumination[i]); + this->GetFunctorVector().push_back(functor); + } + } + +private: + /** Ponderation declaration*/ + VectorType m_Alpha; + VectorType m_Beta; + /** Set the zenithal soalr angle. */ + double m_ZenithalSolarAngle; + /** Flux normalization coefficient. */ + double m_FluxNormalizationCoefficient; + /** Solar illumination value. */ + VectorType m_SolarIllumination; + /** Used to know if the user has set a value for the FluxNormalizationCoefficient parameter + * or if the class has to compute it */ + bool m_IsSetFluxNormalizationCoefficient; + /** Acquisition Day*/ + int m_Day; + /** Acquisition Month*/ + int m_Month; +}; + +} // end namespace otb + +#endif diff --git a/Code/Radiometry/otbReflectanceToLuminanceImageFilter.h b/Code/Radiometry/otbReflectanceToLuminanceImageFilter.h new file mode 100644 index 0000000000..992d783a80 --- /dev/null +++ b/Code/Radiometry/otbReflectanceToLuminanceImageFilter.h @@ -0,0 +1,329 @@ +/*========================================================================= + + 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 __otbReflectanceToLuminanceImageFilter_h +#define __otbReflectanceToLuminanceImageFilter_h + +#include "otb_6S.h" +#include "otbUnaryImageFunctorWithVectorImageFilter.h" +#include "otbVectorImage.h" +#include "otbMath.h" +#include "itkNumericTraits.h" +#include "itkVariableLengthVector.h" +#include "otbMacro.h" +#include "otbOpticalImageMetadataInterface.h" +#include "otbOpticalImageMetadataInterfaceFactory.h" +#include <iomanip> + +namespace otb +{ +namespace Functor +{ +/** + * \class ReflectanceToLuminanceImageFunctor + * \brief Compupute luminance from the reflectance value + * + * Divide by Pi and multiply by an illumination correction coefficient + * and the given solar illumination. + * + * + * \sa ReflectanceToLuminanceImageFilter + * + * \ingroup Functor + * \ingroup Radiometry + * + */ +template <class TInput, class TOutput> +class ReflectanceToLuminanceImageFunctor +{ +public: + ReflectanceToLuminanceImageFunctor() : + m_SolarIllumination(1.0), + m_IlluminationCorrectionCoefficient(1.0) + {} + + virtual ~ReflectanceToLuminanceImageFunctor() {} + + void SetSolarIllumination(double solarIllumination) + { + m_SolarIllumination = solarIllumination; + } + void SetIlluminationCorrectionCoefficient(double coef) + { + m_IlluminationCorrectionCoefficient = coef; + } + + double GetSolarIllumination() + { + return m_SolarIllumination; + } + double GetIlluminationCorrectionCoefficient() + { + return m_IlluminationCorrectionCoefficient; + } + + inline TOutput operator ()(const TInput& inPixel) const + { + TOutput outPixel; + double temp; + temp = static_cast<double>(inPixel) + / static_cast<double>(CONST_PI) + * m_IlluminationCorrectionCoefficient + * m_SolarIllumination; + + outPixel = static_cast<TOutput>(temp); + return outPixel; + } + +private: + double m_SolarIllumination; + double m_IlluminationCorrectionCoefficient; +}; +} + +/** \class ReflectanceToLuminanceImageFilter + * \brief Convert reflectance value into luminance value + * + * Transform a reflectance image into the luminance. For this it uses the + * functor ReflectanceToLuminanceImageFunctor calling for each component of each pixel. + * + * + * For Spot image in the dimap format, the correction parameters are + * retrieved automatically from the metadata + * + * \ingroup ImageToLuminanceImageFunctor + * \ingroup Radiometry + * + * \example Radiometry/AtmosphericCorrectionSequencement.cxx + */ +template <class TInputImage, class TOutputImage> +class ITK_EXPORT ReflectanceToLuminanceImageFilter : + public UnaryImageFunctorWithVectorImageFilter<TInputImage, + TOutputImage, + typename Functor::ReflectanceToLuminanceImageFunctor<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::ReflectanceToLuminanceImageFunctor<typename InputImageType::InternalPixelType, + typename OutputImageType::InternalPixelType> + FunctorType; + + /** "typedef" for standard classes. */ + typedef ReflectanceToLuminanceImageFilter 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(ReflectanceToLuminanceImageFilter, UnaryImageFunctorWithVectorImageFiltermageFilter); + + /** 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 typename itk::VariableLengthVector<double> VectorType; + + /** Image size "typedef" definition. */ + typedef typename InputImageType::SizeType SizeType; + + /** Set the solar illumination value. */ + itkSetMacro(SolarIllumination, VectorType); + /** Give the solar illumination value. */ + itkGetConstReferenceMacro(SolarIllumination, VectorType); + + /** Set the zenithal solar angle. */ + itkSetMacro(ZenithalSolarAngle, double); + /** Give the zenithal solar angle. */ + itkGetConstReferenceMacro(ZenithalSolarAngle, double); + + /** Set/Get the sun elevation angle (internally handled by the zenithal angle)*/ + virtual void SetElevationSolarAngle(double elevationAngle) + { + double zenithalAngle = 90.0 - elevationAngle; + if (this->m_ZenithalSolarAngle != zenithalAngle) + { + this->m_ZenithalSolarAngle = zenithalAngle; + this->Modified(); + } + } + + virtual double GetElevationSolarAngle() const + { + return 90.0 - this->m_ZenithalSolarAngle; + } + + /** Set the day. */ + itkSetClampMacro(Day, int, 1, 31); + /** Give the day. */ + itkGetConstReferenceMacro(Day, int); + + /** Set the mounth. */ + itkSetClampMacro(Month, int, 1, 12); + /** Give the mounth. */ + itkGetConstReferenceMacro(Month, int); + + /** Set the flux normalization coefficient. */ + void SetFluxNormalizationCoefficient(double coef) + { + m_FluxNormalizationCoefficient = coef; + m_IsSetFluxNormalizationCoefficient = true; + this->Modified(); + } + /** Give the flux normalization coefficient. */ + itkGetConstReferenceMacro(FluxNormalizationCoefficient, double); + + /** Set the IsSetFluxNormalizationCoefficient boolean. */ + itkSetMacro(IsSetFluxNormalizationCoefficient, bool); + /** Give the IsSetFluxNormalizationCoefficient boolean. */ + itkGetConstReferenceMacro(IsSetFluxNormalizationCoefficient, bool); + + /** Set the UseClamp boolean. */ + itkSetMacro(UseClamp, bool); + /** Give the UseClamp boolean. */ + itkGetConstReferenceMacro(UseClamp, bool); + +protected: + /** Constructor */ + ReflectanceToLuminanceImageFilter() : + m_ZenithalSolarAngle(120.0), //invalid value which will lead to negative radiometry + m_FluxNormalizationCoefficient(1.), + m_Day(0), + m_Month(0), + m_IsSetFluxNormalizationCoefficient(false) + { + m_SolarIllumination.SetSize(0); + }; + + /** Destructor */ + virtual ~ReflectanceToLuminanceImageFilter() {} + + /** Update the functor list and input parameters */ + virtual void BeforeThreadedGenerateData(void) + { + OpticalImageMetadataInterface::Pointer imageMetadataInterface = OpticalImageMetadataInterfaceFactory::CreateIMI( + this->GetInput()->GetMetaDataDictionary()); + if ((m_Day == 0) && (!m_IsSetFluxNormalizationCoefficient)) + { + m_Day = imageMetadataInterface->GetDay(); + } + + if ((m_Month == 0) && (!m_IsSetFluxNormalizationCoefficient)) + { + m_Month = imageMetadataInterface->GetMonth(); + } + + if (m_SolarIllumination.GetSize() == 0) + { + m_SolarIllumination = imageMetadataInterface->GetSolarIrradiance(); + } + + if (m_ZenithalSolarAngle == 120.0) + { + //the zenithal angle is the complementary of the elevation angle + m_ZenithalSolarAngle = 90.0 - imageMetadataInterface->GetSunElevation(); + } + + std::cout << "Using correction parameters: " << std::endl; + std::cout<< "Day: " << m_Day << std::endl; + std::cout<< "Month: " << m_Month << std::endl; + std::cout<< "Solar irradiance: " << m_SolarIllumination << std::endl; + std::cout<< "Zenithal angle: " << m_ZenithalSolarAngle << std::endl; + + if ((m_SolarIllumination.GetSize() != this->GetInput()->GetNumberOfComponentsPerPixel())) + { + itkExceptionMacro(<< "SolarIllumination parameter should have the same size as the number of bands"); + } + + this->GetFunctorVector().clear(); + + for (unsigned int i = 0; i < this->GetInput()->GetNumberOfComponentsPerPixel(); ++i) + { + FunctorType functor; + double coefTemp = 0.; + if (!m_IsSetFluxNormalizationCoefficient) + { + if (m_Day * m_Month != 0 && m_Day < 32 && m_Month < 13) + { + otb_6s_doublereal dsol = 0.; + otb_6s_integer day = static_cast<otb_6s_integer>(m_Day); + otb_6s_integer month = static_cast<otb_6s_integer>(m_Month); + //int cr(0); + otb_6s_varsol_(&day, &month, &dsol); + coefTemp = vcl_cos(m_ZenithalSolarAngle * CONST_PI_180) * static_cast<double>(dsol); + } + else + { + itkExceptionMacro(<< "Day has to be included between 1 and 31, Month beetween 1 and 12."); + } + } + else + { + coefTemp = + vcl_cos(m_ZenithalSolarAngle * + CONST_PI_180) * m_FluxNormalizationCoefficient * m_FluxNormalizationCoefficient; + } + functor.SetIlluminationCorrectionCoefficient(coefTemp); + functor.SetSolarIllumination(static_cast<double>(m_SolarIllumination[i])); + + this->GetFunctorVector().push_back(functor); + } + } + +private: + + /** Set the zenithal soalr angle. */ + double m_ZenithalSolarAngle; + /** Flux normalization coefficient. */ + double m_FluxNormalizationCoefficient; + /** Acquisition day. */ + int m_Day; + /** Acquisition mounth. */ + int m_Month; + /** Solar illumination value. */ + VectorType m_SolarIllumination; + /** Used to know if the user has set a value for the FluxNormalizationCoefficient parameter + * or if the class has to compute it */ + bool m_IsSetFluxNormalizationCoefficient; + /** Clamp values to [0,1] */ + bool m_UseClamp; + +}; + +} // end namespace otb +#endif diff --git a/Testing/Applications/Radiometry/CMakeLists.txt b/Testing/Applications/Radiometry/CMakeLists.txt index 5ecc9561ee..6f685b89de 100644 --- a/Testing/Applications/Radiometry/CMakeLists.txt +++ b/Testing/Applications/Radiometry/CMakeLists.txt @@ -1,3 +1,37 @@ +OTB_TEST_APPLICATION(NAME apTvRaOpticalCalibration_UnknownSensor + APP OpticalCalibration + OPTIONS + -in /home/cpalmann/Desktop/SPOT6/SPOT6_Primary_Pan-sharpened_12_bits/PROD_SPOT6_001/VOL_SPOT6_001_A/IMG_SPOT6_PMS_001_A/SPOT6_500x500.tif?&skipgeom=true + -out ${TEMP}/apTvRaOpticalCalibration_SPOT6.tif + -level toa + -acquisition.gainbias ${BASELINE}/../Files/apTvRaOpticalCalibrationUnknownSensorGainsBiases.txt + -acquisition.day 3 + -acquisition.month 6 + -acquisition.sunelevationangle 63.8861327469 + -acquisition.solarilluminations ${BASELINE}/../Files/apTvRaOpticalCalibrationUnknownSensorSolarIllumations.txt + -milli true + -clamp false + VALID --compare-image ${EPSILON_6} + ${BASELINE}/apTvRaOpticalCalibration_SPOT6.tif + ${TEMP}/apTvRaOpticalCalibration_SPOT6.tif ) + +OTB_TEST_APPLICATION(NAME apTvRaOpticalCalibration_Reverse_UnknownSensor + APP OpticalCalibration + OPTIONS + -in ${BASELINE}/apTvRaOpticalCalibration_SPOT6.tif + -out ${TEMP}/apTvRaOpticalCalibration_Rev_SPOT6.tif + -level toatoim + -acquisition.gainbias ${BASELINE}/../Files/apTvRaOpticalCalibrationUnknownSensorGainsBiases.txt + -acquisition.day 3 + -acquisition.month 6 + -acquisition.sunelevationangle 63.8861327469 + -acquisition.solarilluminations ${BASELINE}/../Files/apTvRaOpticalCalibrationUnknownSensorSolarIllumations.txt + -milli true + -clamp false + VALID --compare-image ${EPSILON_6} + /home/cpalmann/Desktop/SPOT6/SPOT6_Primary_Pan-sharpened_12_bits/PROD_SPOT6_001/VOL_SPOT6_001_A/IMG_SPOT6_PMS_001_A/SPOT6_500x500.tif + ${TEMP}/apTvRaOpticalCalibration_Rev_SPOT6.tif ) + OTB_TEST_APPLICATION(NAME apTvRaOpticalCalibration_QuickbirdPAN APP OpticalCalibration diff --git a/Testing/Code/Radiometry/CMakeLists.txt b/Testing/Code/Radiometry/CMakeLists.txt index 2faef1103b..39135effa2 100644 --- a/Testing/Code/Radiometry/CMakeLists.txt +++ b/Testing/Code/Radiometry/CMakeLists.txt @@ -95,7 +95,7 @@ add_test(raTvImageToLuminanceImageFilterAutoSpot5 ${RADIOMETRY_TESTS2} ${TEMP}/raTvImageToLuminanceImageFilterAutoSpot5.tif ) -add_test(raTvImageToLuminanceImageFilterAutoIkonos ${RADIOMETRY_TESTS2} +add_test(raTvImageToLuminanceImageFilterAutoIkonos ${RADIOMETRY_TESTS2} --compare-image ${EPSILON_12} ${BASELINE}/raTvImageToLuminanceImageFilterAutoIkonos.tif ${TEMP}/raTvImageToLuminanceImageFilterAutoIkonos.tif otbImageToLuminanceImageFilterAuto @@ -340,6 +340,407 @@ add_test(raTuImageToReflectanceImageFilterAutoFORMOSAT2 ${RADIOMETRY_TESTS2} endif() + +# ------- otb::LuminanceToImageImageFilter ------------------------------ +add_test(raTuLuminanceToImageImageFilterNew ${RADIOMETRY_TESTS2} + otbLuminanceToImageImageFilterNew +) + +add_test(raTvLuminanceToImageImageFilter ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} ${INPUTDATA}/verySmallFSATSW.tif + ${TEMP}/raTvverySmallFSATSWImageFilter.tif + otbLuminanceToImageImageFilter + ${BASELINE}/raTvImageToLuminanceImageFilter.tif + ${TEMP}/raTvverySmallFSATSWImageFilter.tif + 10 #channel 1 alpha + 20 #channel 2 alpha + 30 #channel 3 alpha + 40 #channel 4 alpha + 1 #channel 1 beta + 2 #channel 2 beta + 3 #channel 3 beta + 4 #channel 4 beta + ) + +if(OTB_DATA_USE_LARGEINPUT) + +set(TEHERANSPOT5DIR ${TEMP}/OpticalCalibLumToImSPOT5) +file(MAKE_DIRECTORY ${TEHERANSPOT5DIR}) +file(GLOB MTDATA ${LARGEINPUT}/SPOT5/TEHERAN/*.DIM) +foreach(f ${MTDATA}) + configure_file(${f} ${TEHERANSPOT5DIR} COPYONLY) +endforeach(f) +configure_file(${BASELINE}/raTvImageToLuminanceImageFilterAutoSpot5.tif ${TEHERANSPOT5DIR}/IMAGERY.TIF COPYONLY) +add_test(raTvLuminanceToImageImageFilterAutoSpot5 ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} ${TEHERANSPOT5DIR}/raTvSMALLTEHERANLuminanceToImageImageFilterAutoSpot5.tif + ${TEHERANSPOT5DIR}/raTvLuminanceToImageImageFilterAutoSpot5.tif + otbLuminanceToImageImageFilterAuto + ${TEHERANSPOT5DIR}/IMAGERY.TIF + ${LARGEINPUT}/SPOT5/TEHERAN/IMAGERY.TIF + ${TEHERANSPOT5DIR}/raTvLuminanceToImageImageFilterAutoSpot5.tif + ${TEHERANSPOT5DIR}/raTvSMALLTEHERANLuminanceToImageImageFilterAutoSpot5.tif + ) + +set(BLOSSEVILLEIKONOSDIR ${TEMP}/OpticalCalibLumToImIkonos) +file(MAKE_DIRECTORY ${BLOSSEVILLEIKONOSDIR}) +file(GLOB MTDATA ${LARGEINPUT}/IKONOS/BLOSSEVILLE/*metadata.txt ${LARGEINPUT}/IKONOS/BLOSSEVILLE/*pan*) +foreach(f ${MTDATA}) + configure_file(${f} ${BLOSSEVILLEIKONOSDIR} COPYONLY) +endforeach(f) +configure_file(${BASELINE}/raTvImageToLuminanceImageFilterAutoIkonos.tif ${BLOSSEVILLEIKONOSDIR}/po_2619900_pan_0000000.tif COPYONLY) +add_test(raTvLuminanceToImageImageFilterAutoIkonos ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} ${BLOSSEVILLEIKONOSDIR}/raTvSMALLBLOSSEVILLELuminanceToImageImageFilterAutoIkonos.tif + ${BLOSSEVILLEIKONOSDIR}/raTvLuminanceToImageImageFilterAutoIkonos.tif + otbLuminanceToImageImageFilterAuto + ${BLOSSEVILLEIKONOSDIR}/po_2619900_pan_0000000.tif + ${LARGEINPUT}/IKONOS/BLOSSEVILLE/po_2619900_pan_0000000.tif + ${BLOSSEVILLEIKONOSDIR}/raTvLuminanceToImageImageFilterAutoIkonos.tif + ${BLOSSEVILLEIKONOSDIR}/raTvSMALLBLOSSEVILLELuminanceToImageImageFilterAutoIkonos.tif + ) + +set(ROMEWV2DIR ${TEMP}/OpticalCalibLumToImWv2PAN) +file(MAKE_DIRECTORY ${ROMEWV2DIR}) +set(MTDATADIR ${LARGEINPUT}/WORLDVIEW2/ROME/WV-2_standard_8band_bundle_16bit/052298844010_01_P001_PAN) +file(GLOB MTDATA ${MTDATADIR}/*TIL ${MTDATADIR}/*RPB ${MTDATADIR}/*XML ${MTDATADIR}/*IMD) +foreach(f ${MTDATA}) + configure_file(${f} ${ROMEWV2DIR} COPYONLY) +endforeach(f) +configure_file(${BASELINE}/raTvImageToLuminanceImageFilterAutoWV2PAN.tif ${ROMEWV2DIR}/09DEC10103019-P2AS-052298844010_01_P001.TIF COPYONLY) +add_test(raTvLuminanceToImageImageFilterAutoWV2PAN ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} ${ROMEWV2DIR}/raTvSMALLROMELuminanceToImageImageFilterAutoWV2PAN.tif + ${ROMEWV2DIR}/raTvLuminanceToImageImageFilterAutoWV2PAN.tif + otbLuminanceToImageImageFilterAuto + ${ROMEWV2DIR}/09DEC10103019-P2AS-052298844010_01_P001.TIF + ${LARGEINPUT}/WORLDVIEW2/ROME/WV-2_standard_8band_bundle_16bit/052298844010_01_P001_PAN/09DEC10103019-P2AS-052298844010_01_P001.TIF + ${ROMEWV2DIR}/raTvLuminanceToImageImageFilterAutoWV2PAN.tif + ${ROMEWV2DIR}/raTvSMALLROMELuminanceToImageImageFilterAutoWV2PAN.tif + ) + +set(ROMEWV2DIR ${TEMP}/OpticalCalibLumToImWv2MULTI) +file(MAKE_DIRECTORY ${ROMEWV2DIR}) +set(MTDATADIR ${LARGEINPUT}/WORLDVIEW2/ROME/WV-2_standard_8band_bundle_16bit/052298844010_01_P001_MUL) +file(GLOB MTDATA ${MTDATADIR}/*TIL ${MTDATADIR}/*RPB ${MTDATADIR}/*XML ${MTDATADIR}/*IMD) +foreach(f ${MTDATA}) + configure_file(${f} ${ROMEWV2DIR} COPYONLY) +endforeach(f) +configure_file(${BASELINE}/raTvImageToLuminanceImageFilterAutoWV2Multi.tif ${ROMEWV2DIR}/09DEC10103019-M2AS-052298844010_01_P001.TIF COPYONLY) +add_test(raTvLuminanceToImageImageFilterAutoWV2MULTI ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} ${ROMEWV2DIR}/raTvSMALLROMELuminanceToImageImageFilterAutoWV2Multi.tif + ${ROMEWV2DIR}/raTvLuminanceToImageImageFilterAutoWV2Multi.tif + otbLuminanceToImageImageFilterAuto + ${ROMEWV2DIR}/09DEC10103019-M2AS-052298844010_01_P001.TIF + ${LARGEINPUT}/WORLDVIEW2/ROME/WV-2_standard_8band_bundle_16bit/052298844010_01_P001_MUL/09DEC10103019-M2AS-052298844010_01_P001.TIF + ${ROMEWV2DIR}/raTvLuminanceToImageImageFilterAutoWV2Multi.tif + ${ROMEWV2DIR}/raTvSMALLROMELuminanceToImageImageFilterAutoWV2Multi.tif + ) + +set(TOULOUSEQBDIR ${TEMP}/OpticalCalibLumToImQBPAN) +file(MAKE_DIRECTORY ${TOULOUSEQBDIR}) +set(MTDATADIR ${LARGEINPUT}/QUICKBIRD/TOULOUSE/000000128955_01_P001_PAN) +file(GLOB MTDATA ${MTDATADIR}/*TIL ${MTDATADIR}/*RPB ${MTDATADIR}/*XML ${MTDATADIR}/*IMD) +foreach(f ${MTDATA}) + configure_file(${f} ${TOULOUSEQBDIR} COPYONLY) +endforeach(f) +configure_file(${BASELINE}/raTvImageToLuminanceImageFilterAutoQuickbirdPAN.tif ${TOULOUSEQBDIR}/02APR01105228-P1BS-000000128955_01_P001.TIF COPYONLY) +add_test(raTvLuminanceToImageImageFilterAutoQuickbirdPAN ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} ${TOULOUSEQBDIR}/raTvSMALLTOULOUSELuminanceToImageImageFilterAutoQuickbirdPAN.tif + ${TOULOUSEQBDIR}/raTvLuminanceToImageImageFilterAutoQuickbirdPAN.tif + otbLuminanceToImageImageFilterAuto + ${TOULOUSEQBDIR}/02APR01105228-P1BS-000000128955_01_P001.TIF + ${LARGEINPUT}/QUICKBIRD/TOULOUSE/000000128955_01_P001_PAN/02APR01105228-P1BS-000000128955_01_P001.TIF + ${TOULOUSEQBDIR}/raTvLuminanceToImageImageFilterAutoQuickbirdPAN.tif + ${TOULOUSEQBDIR}/raTvSMALLTOULOUSELuminanceToImageImageFilterAutoQuickbirdPAN.tif + ) + +set(TOULOUSEQBDIR ${TEMP}/OpticalCalibLumToImQBXS) +file(MAKE_DIRECTORY ${TOULOUSEQBDIR}) +set(MTDATADIR ${LARGEINPUT}/QUICKBIRD/TOULOUSE/000000128955_01_P001_MUL) +file(GLOB MTDATA ${MTDATADIR}/*TIL ${MTDATADIR}/*RPB ${MTDATADIR}/*XML ${MTDATADIR}/*IMD) +foreach(f ${MTDATA}) + configure_file(${f} ${TOULOUSEQBDIR} COPYONLY) +endforeach(f) +configure_file(${BASELINE}/raTvImageToLuminanceImageFilterAutoQuickbirdXS.tif ${TOULOUSEQBDIR}/02APR01105228-M1BS-000000128955_01_P001.TIF COPYONLY) +add_test(raTvLuminanceToImageImageFilterAutoQuickbirdXS ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} ${TOULOUSEQBDIR}/raTvSMALLTOULOUSELuminanceToImageImageFilterAutoQuickbirdXS.tif + ${TOULOUSEQBDIR}/raTvLuminanceToImageImageFilterAutoQuickbirdXS.tif + otbLuminanceToImageImageFilterAuto + ${TOULOUSEQBDIR}/02APR01105228-M1BS-000000128955_01_P001.TIF + ${LARGEINPUT}/QUICKBIRD/TOULOUSE/000000128955_01_P001_MUL/02APR01105228-M1BS-000000128955_01_P001.TIF + ${TOULOUSEQBDIR}/raTvLuminanceToImageImageFilterAutoQuickbirdXS.tif + ${TOULOUSEQBDIR}/raTvSMALLTOULOUSELuminanceToImageImageFilterAutoQuickbirdXS.tif + ) + +set(SUDOUESTFORMOSATDIR ${TEMP}/OpticalCalibLumToImFormosat) +file(MAKE_DIRECTORY ${SUDOUESTFORMOSATDIR}) +set(MTDATADIR ${LARGEINPUT}/FORMOSAT/Sudouest_20071013_MS_fmsat) +file(GLOB MTDATA ${MTDATADIR}/*.DIM) +foreach(f ${MTDATA}) + configure_file(${f} ${SUDOUESTFORMOSATDIR} COPYONLY) +endforeach(f) +configure_file(${BASELINE}/raTvImageToLuminanceImageFilterAutoFormosat.tif ${SUDOUESTFORMOSATDIR}/IMAGERY.TIF COPYONLY) +add_test(raTvLuminanceToImageImageFilterAutoFORMOSAT ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} ${SUDOUESTFORMOSATDIR}/raTvSMALLSOLuminanceToImageImageFilterAutoFormosat.tif + ${SUDOUESTFORMOSATDIR}/raTvLuminanceToImageImageFilterAutoFormosat.tif + otbLuminanceToImageImageFilterAuto + ${SUDOUESTFORMOSATDIR}/IMAGERY.TIF + ${LARGEINPUT}/FORMOSAT/Sudouest_20071013_MS_fmsat/IMAGERY.TIF + ${SUDOUESTFORMOSATDIR}/raTvLuminanceToImageImageFilterAutoFormosat.tif + ${SUDOUESTFORMOSATDIR}/raTvSMALLSOLuminanceToImageImageFilterAutoFormosat.tif + ) + +endif() + + +# ------- otb::ReflectanceToLuminanceImageFilter ------------------------------ +add_test(raTuReflectanceToLuminanceImageFilterNew ${RADIOMETRY_TESTS2} + otbReflectanceToLuminanceImageFilterNew +) + +add_test(raTvReflectanceToLuminanceImageFilter ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} ${INPUTDATA}/verySmallFSATSW.tif + ${TEMP}/raTvverySmallFSATSWImageFilterDsol.tif + otbReflectanceToLuminanceImageFilter + ${BASELINE}/raTvLuminanceToReflectanceImageFilterDsol.tif + ${TEMP}/raTvverySmallFSATSWImageFilterDsol.tif + 0.2 #radius + 10 #channel 1 illumination + 20 #channel 2 illumination + 30 #channel 3 illumination + 40 #channel 4 illumination + 0.9923885328 +) + +add_test(raTvReflectanceToLuminanceImageFilterDayMonth ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} ${INPUTDATA}/verySmallFSATSW.tif + ${TEMP}/raTvverySmallFSATSWImageFilterDayMounth.tif + otbReflectanceToLuminanceImageFilter + ${BASELINE}/raTvLuminanceToReflectanceImageFilter.tif + ${TEMP}/raTvverySmallFSATSWImageFilterDayMounth.tif + 0.2 #radius + 10 #channel 1 illumination + 20 #channel 2 illumination + 30 #channel 3 illumination + 40 #channel 4 illumination + 3 #day + 5 #mounth + ) + +if(OTB_DATA_USE_LARGEINPUT) + +set(TEHERANSPOT5DIR ${TEMP}/OpticalCalibRefToLumSPOT5) +file(MAKE_DIRECTORY ${TEHERANSPOT5DIR}) +file(GLOB MTDATA ${LARGEINPUT}/SPOT5/TEHERAN/*.DIM) +foreach(f ${MTDATA}) + configure_file(${f} ${TEHERANSPOT5DIR} COPYONLY) +endforeach(f) +configure_file(${BASELINE}/raTvLuminanceToReflectanceImageFilterAutoSpot5.tif ${TEHERANSPOT5DIR}/IMAGERY.TIF COPYONLY) +add_test(raTvReflectanceToLuminanceImageFilterAutoSpot5 ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} ${BASELINE}/raTvImageToLuminanceImageFilterAutoSpot5.tif + ${TEHERANSPOT5DIR}/raTvReflectanceToLuminanceImageFilterAutoSpot5.tif + otbReflectanceToLuminanceImageFilterAuto + ${TEHERANSPOT5DIR}/IMAGERY.TIF + ${TEHERANSPOT5DIR}/raTvReflectanceToLuminanceImageFilterAutoSpot5.tif + ) + +set(BLOSSEVILLEIKONOSDIR ${TEMP}/OpticalCalibRefToLumIkonos) +file(MAKE_DIRECTORY ${BLOSSEVILLEIKONOSDIR}) +file(GLOB MTDATA ${LARGEINPUT}/IKONOS/BLOSSEVILLE/*metadata.txt ${LARGEINPUT}/IKONOS/BLOSSEVILLE/*pan*) +foreach(f ${MTDATA}) + configure_file(${f} ${BLOSSEVILLEIKONOSDIR} COPYONLY) +endforeach(f) +configure_file(${BASELINE}/raTvLuminanceToReflectanceImageFilterAutoIkonos.tif ${BLOSSEVILLEIKONOSDIR}/po_2619900_pan_0000000.tif COPYONLY) +add_test(raTvReflectanceToLuminanceImageFilterAutoIkonos ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} ${BASELINE}/raTvImageToLuminanceImageFilterAutoIkonos.tif + ${BLOSSEVILLEIKONOSDIR}/raTvReflectanceToLuminanceImageFilterAutoIkonos.tif + otbReflectanceToLuminanceImageFilterAuto + ${BLOSSEVILLEIKONOSDIR}/po_2619900_pan_0000000.tif + ${BLOSSEVILLEIKONOSDIR}/raTvReflectanceToLuminanceImageFilterAutoIkonos.tif + ) + +set(ROMEWV2DIR ${TEMP}/OpticalCalibRefToLumWv2PAN) +file(MAKE_DIRECTORY ${ROMEWV2DIR}) +set(MTDATADIR ${LARGEINPUT}/WORLDVIEW2/ROME/WV-2_standard_8band_bundle_16bit/052298844010_01_P001_PAN) +file(GLOB MTDATA ${MTDATADIR}/*TIL ${MTDATADIR}/*RPB ${MTDATADIR}/*XML ${MTDATADIR}/*IMD) +foreach(f ${MTDATA}) + configure_file(${f} ${ROMEWV2DIR} COPYONLY) +endforeach(f) +configure_file(${BASELINE}/raTvLuminanceToReflectanceImageFilterAutoWV2PAN.tif ${ROMEWV2DIR}/09DEC10103019-P2AS-052298844010_01_P001.TIF COPYONLY) +add_test(raTvReflectanceToLuminanceImageFilterAutoWV2PAN ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} ${BASELINE}/raTvImageToLuminanceImageFilterAutoWV2PAN.tif + ${ROMEWV2DIR}/raTvReflectanceToLuminanceImageFilterAutoWV2PAN.tif + otbReflectanceToLuminanceImageFilterAuto + ${ROMEWV2DIR}/09DEC10103019-P2AS-052298844010_01_P001.TIF + ${ROMEWV2DIR}/raTvReflectanceToLuminanceImageFilterAutoWV2PAN.tif + ) + +set(ROMEWV2DIR ${TEMP}/OpticalCalibRefToLumWv2MULTI) +file(MAKE_DIRECTORY ${ROMEWV2DIR}) +set(MTDATADIR ${LARGEINPUT}/WORLDVIEW2/ROME/WV-2_standard_8band_bundle_16bit/052298844010_01_P001_MUL) +file(GLOB MTDATA ${MTDATADIR}/*TIL ${MTDATADIR}/*RPB ${MTDATADIR}/*XML ${MTDATADIR}/*IMD) +foreach(f ${MTDATA}) + configure_file(${f} ${ROMEWV2DIR} COPYONLY) +endforeach(f) +configure_file(${BASELINE}/raTvLuminanceToReflectanceImageFilterAutoWV2Multi.tif ${ROMEWV2DIR}/09DEC10103019-M2AS-052298844010_01_P001.TIF COPYONLY) +add_test(raTvReflectanceToLuminanceImageFilterAutoWV2MULTI ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} ${BASELINE}/raTvImageToLuminanceImageFilterAutoWV2Multi.tif + ${ROMEWV2DIR}/raTvReflectanceToLuminanceImageFilterAutoWV2Multi.tif + otbReflectanceToLuminanceImageFilterAuto + ${ROMEWV2DIR}/09DEC10103019-M2AS-052298844010_01_P001.TIF + ${ROMEWV2DIR}/raTvReflectanceToLuminanceImageFilterAutoWV2Multi.tif + ) + +set(TOULOUSEQBDIR ${TEMP}/OpticalCalibRefToLumQBPAN) +file(MAKE_DIRECTORY ${TOULOUSEQBDIR}) +set(MTDATADIR ${LARGEINPUT}/QUICKBIRD/TOULOUSE/000000128955_01_P001_PAN) +file(GLOB MTDATA ${MTDATADIR}/*TIL ${MTDATADIR}/*RPB ${MTDATADIR}/*XML ${MTDATADIR}/*IMD) +foreach(f ${MTDATA}) + configure_file(${f} ${TOULOUSEQBDIR} COPYONLY) +endforeach(f) +configure_file(${BASELINE}/raTvLuminanceToReflectanceImageFilterAutoQuickbirdPAN.tif ${TOULOUSEQBDIR}/02APR01105228-P1BS-000000128955_01_P001.TIF COPYONLY) +add_test(raTvReflectanceToLuminanceImageFilterAutoQuickbirdPAN ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} ${BASELINE}/raTvImageToLuminanceImageFilterAutoQuickbirdPAN.tif + ${TOULOUSEQBDIR}/raTvReflectanceToLuminanceImageFilterAutoQuickbirdPAN.tif + otbReflectanceToLuminanceImageFilterAuto + ${TOULOUSEQBDIR}/02APR01105228-P1BS-000000128955_01_P001.TIF + ${TOULOUSEQBDIR}/raTvReflectanceToLuminanceImageFilterAutoQuickbirdPAN.tif + ) + +set(TOULOUSEQBDIR ${TEMP}/OpticalCalibRefToLumQBXS) +file(MAKE_DIRECTORY ${TOULOUSEQBDIR}) +set(MTDATADIR ${LARGEINPUT}/QUICKBIRD/TOULOUSE/000000128955_01_P001_MUL) +file(GLOB MTDATA ${MTDATADIR}/*TIL ${MTDATADIR}/*RPB ${MTDATADIR}/*XML ${MTDATADIR}/*IMD) +foreach(f ${MTDATA}) + configure_file(${f} ${TOULOUSEQBDIR} COPYONLY) +endforeach(f) +configure_file(${BASELINE}/raTvLuminanceToReflectanceImageFilterAutoQuickbirdXS.tif ${TOULOUSEQBDIR}/02APR01105228-M1BS-000000128955_01_P001.TIF COPYONLY) +add_test(raTvReflectanceToLuminanceImageFilterAutoQuickbirdXS ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} ${BASELINE}/raTvImageToLuminanceImageFilterAutoQuickbirdXS.tif + ${TOULOUSEQBDIR}/raTvReflectanceToLuminanceImageFilterAutoQuickbirdXS.tif + otbReflectanceToLuminanceImageFilterAuto + ${TOULOUSEQBDIR}/02APR01105228-M1BS-000000128955_01_P001.TIF + ${TOULOUSEQBDIR}/raTvReflectanceToLuminanceImageFilterAutoQuickbirdXS.tif + ) + +# ${BASELINE}/raTvLuminanceToReflectanceImageFilterAutoFormosat.tif : this image seems to be corrupted +#set(SUDOUESTFORMOSATDIR ${TEMP}/OpticalCalibRefToLumFormosat) +#file(MAKE_DIRECTORY ${SUDOUESTFORMOSATDIR}) +#set(MTDATADIR ${LARGEINPUT}/FORMOSAT/Sudouest_20071013_MS_fmsat) +#file(GLOB MTDATA ${MTDATADIR}/*.DIM) +#foreach(f ${MTDATA}) +# configure_file(${f} ${SUDOUESTFORMOSATDIR} COPYONLY) +#endforeach(f) +#configure_file(${BASELINE}/raTvLuminanceToReflectanceImageFilterAutoFormosat.tif ${SUDOUESTFORMOSATDIR}/IMAGERY.TIF COPYONLY) +#add_test(raTvReflectanceToLuminanceImageFilterAutoFORMOSAT ${RADIOMETRY_TESTS2} +# --compare-image ${EPSILON_12} ${BASELINE}/raTvImageToLuminanceImageFilterAutoFormosat.tif +# ${SUDOUESTFORMOSATDIR}/raTvReflectanceToLuminanceImageFilterAutoFormosat.tif +# otbReflectanceToLuminanceImageFilterAuto +# ${SUDOUESTFORMOSATDIR}/IMAGERY.TIF +# ${SUDOUESTFORMOSATDIR}/raTvReflectanceToLuminanceImageFilterAutoFormosat.tif +# ) + +endif() + +# ------- otb::ReflectanceToImageImageFilter ------------------------------ +add_test(raTuReflectanceToImageImageFilterNew ${RADIOMETRY_TESTS2} + otbReflectanceToImageImageFilterNew +) + +add_test(raTvReflectanceToImageImageFilter ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} ${INPUTDATA}/verySmallFSATSW.tif + ${TEMP}/raTvReflectanceToImageImageFilterDsol.tif + otbReflectanceToImageImageFilter + ${BASELINE}/raTvImageToReflectanceImageFilterDsol.tif + ${TEMP}/raTvReflectanceToImageImageFilterDsol.tif + 0.2 #radius + 1 #channel 1 alpha + 2 #channel 2 alpha + 3 #channel 3 alpha + 4 #channel 4 alpha + 10 #channel 1 beta + 11 #channel 2 beta + 12 #channel 3 beta + 13 #channel 4 beta + 10 #channel 1 illumination + 20 #channel 2 illumination + 30 #channel 3 illumination + 40 #channel 4 illumination + 0.9923885328 #d/d0 corresponding to the date 03/05 + ) + +add_test(raTvRomaniaReflectanceToImage ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} + ${INPUTDATA}/Romania_Extract.tif + ${TEMP}/raTvRomaniaImage.tif + otbReflectanceToImageImageFilter + ${BASELINE}/raTvRomania_Reflectance.tif + ${TEMP}/raTvRomaniaImage.tif + 27.3 # = 90-62.70 : elevation et azimuth solaire + 1.981247824 #channel 1 alpha = 0.881338*2.24800 + 4.332207085 #channel 2 alpha = 0.858713*5.04500 + 2.32064768 #channel 3 alpha = 0.685568*3.38500 + 9.3177861 #channel 4 alpha = 6.19122*1.50500 + 0 #channel 1 beta + 0 #channel 2 beta + 0 #channel 3 beta + 0 #channel 4 beta + 1057.56 #channel 1 illumination + 1570.23 #channel 2 illumination + 1842.94 #channel 3 illumination + 232.820 #channel 4 illumination + 0.9889145564708814 #= sqrt(0.977952) d/d0 corresponding to the date 03/05 + ) + +add_test(raTvReflectanceToImageImageFilterDayMounth ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} ${INPUTDATA}/verySmallFSATSW.tif + ${TEMP}/raTvReflectanceToImageImageFilterDayMounth.tif + otbReflectanceToImageImageFilter + ${BASELINE}/raTvImageToReflectanceImageFilter.tif + ${TEMP}/raTvReflectanceToImageImageFilterDayMounth.tif + 0.2 #radius + 1 #channel 1 alpha + 2 #channel 2 alpha + 3 #channel 3 alpha + 4 #channel 4 alpha + 10 #channel 1 beta + 11 #channel 2 beta + 12 #channel 3 beta + 13 #channel 4 beta + 10 #channel 1 illumination + 20 #channel 2 illumination + 30 #channel 3 illumination + 40 #channel 4 illumination + 3 #day + 5 #mounth + ) + +if(OTB_DATA_USE_LARGEINPUT) + +set(TEHERANSPOT5DIR ${TEMP}/OpticalCalibRefToImSPOT5) +file(MAKE_DIRECTORY ${TEHERANSPOT5DIR}) +file(GLOB MTDATA ${LARGEINPUT}/SPOT5/TEHERAN/*.DIM) +foreach(f ${MTDATA}) + configure_file(${f} ${TEHERANSPOT5DIR} COPYONLY) +endforeach(f) +configure_file(${BASELINE}/raTvImageToReflectanceImageFilterAuto.tif ${TEHERANSPOT5DIR}/IMAGERY.TIF COPYONLY) +add_test(raTvReflectanceToImageImageFilterAutoSpot5 ${RADIOMETRY_TESTS2} + --compare-image ${EPSILON_12} ${TEHERANSPOT5DIR}/raTvSMALLTEHERANReflectanceToImageImageFilterAutoSpot5.tif + ${TEHERANSPOT5DIR}/raTvReflectanceToImageImageFilterAutoSpot5.tif + otbReflectanceToImageImageFilterAuto + ${TEHERANSPOT5DIR}/IMAGERY.TIF + ${LARGEINPUT}/SPOT5/TEHERAN/IMAGERY.TIF + ${TEHERANSPOT5DIR}/raTvReflectanceToImageImageFilterAutoSpot5.tif + ${TEHERANSPOT5DIR}/raTvSMALLTEHERANReflectanceToImageImageFilterAutoSpot5.tif + ) + + +add_test(raTuReflectanceToImageImageFilterAutoFORMOSAT2 ${RADIOMETRY_TESTS2} + otbReflectanceToImageImageFilterAuto + ${LARGEINPUT}/FORMOSAT/Sudouest_20071013_MS_fmsat/IMAGERY.TIF + ${LARGEINPUT}/FORMOSAT/Sudouest_20071013_MS_fmsat/IMAGERY.TIF + ${TEMP}/raTvReflectanceToImageImageFilterAutoFORMOSAT2.tif + ${TEMP}/raTvSMALLReflectanceToImageImageFilterAutoFORMOSAT2.tif + ) +endif() # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ otbRADIOMETRY_TESTS3 ~~~~~~~~~~~~~~~~~~~~~~~~~~~ # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -1693,6 +2094,15 @@ otbLuminanceToReflectanceImageFilterAuto.cxx otbImageToReflectanceImageFilterNew.cxx otbImageToReflectanceImageFilter.cxx otbImageToReflectanceImageFilterAuto.cxx +otbLuminanceToImageImageFilterNew.cxx +otbLuminanceToImageImageFilter.cxx +otbLuminanceToImageImageFilterAuto.cxx +otbReflectanceToLuminanceImageFilterNew.cxx +otbReflectanceToLuminanceImageFilter.cxx +otbReflectanceToLuminanceImageFilterAuto.cxx +otbReflectanceToImageImageFilterNew.cxx +otbReflectanceToImageImageFilter.cxx +otbReflectanceToImageImageFilterAuto.cxx ) set(Radiometry_SRCS3 otbRadiometryTests3.cxx diff --git a/Testing/Code/Radiometry/otbLuminanceToImageImageFilter.cxx b/Testing/Code/Radiometry/otbLuminanceToImageImageFilter.cxx new file mode 100644 index 0000000000..b2aaa9f04d --- /dev/null +++ b/Testing/Code/Radiometry/otbLuminanceToImageImageFilter.cxx @@ -0,0 +1,68 @@ +/*========================================================================= + + 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 "otbLuminanceToImageImageFilter.h" +#include "otbVectorImage.h" +#include "otbImageFileReader.h" +#include "otbImageFileWriter.h" +#include "itkVariableLengthVector.h" + +int otbLuminanceToImageImageFilter(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::LuminanceToImageImageFilter<InputImageType, OutputImageType> LuminanceToImageImageFilterType; + typedef LuminanceToImageImageFilterType::VectorType VectorType; + + ReaderType::Pointer reader = ReaderType::New(); + WriterType::Pointer writer = WriterType::New(); + reader->SetFileName(inputFileName); + writer->SetFileName(outputFileName); + reader->UpdateOutputInformation(); + + unsigned int nbOfComponent = reader->GetOutput()->GetNumberOfComponentsPerPixel(); + + VectorType alpha(nbOfComponent); + VectorType beta(nbOfComponent); + alpha.Fill(0); + beta.Fill(0); + + for (unsigned int i = 0; i < nbOfComponent; ++i) + { + alpha[i] = static_cast<double>(atof(argv[i + 3])); + beta[i] = static_cast<double>(atof(argv[i + 7])); + } + + // Instantiating object + LuminanceToImageImageFilterType::Pointer filter = LuminanceToImageImageFilterType::New(); + filter->SetAlpha(alpha); + filter->SetBeta(beta); + filter->SetInput(reader->GetOutput()); + writer->SetInput(filter->GetOutput()); + writer->Update(); + + return EXIT_SUCCESS; +} diff --git a/Testing/Code/Radiometry/otbLuminanceToImageImageFilterAuto.cxx b/Testing/Code/Radiometry/otbLuminanceToImageImageFilterAuto.cxx new file mode 100644 index 0000000000..c62d4d559e --- /dev/null +++ b/Testing/Code/Radiometry/otbLuminanceToImageImageFilterAuto.cxx @@ -0,0 +1,74 @@ +/*========================================================================= + + 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 "otbLuminanceToImageImageFilter.h" +#include "otbVectorImage.h" +#include "otbImageFileReader.h" +#include "otbImageFileWriter.h" +#include "itkVariableLengthVector.h" +#include "otbMultiChannelExtractROI.h" + +//Test the retrieval of parameters from the image metadata +int otbLuminanceToImageImageFilterAuto(int argc, char * argv[]) +{ + const char * inputFileName = argv[1]; + const char * inputFileName2 = argv[2]; + const char * outputFileName = argv[3]; + const char * outputFileName2 = argv[4]; + + 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::LuminanceToImageImageFilter<InputImageType, OutputImageType> LuminanceToImageImageFilterType; + typedef LuminanceToImageImageFilterType::VectorType VectorType; + typedef otb::MultiChannelExtractROI<PixelType, PixelType> RoiFilterType; + + ReaderType::Pointer reader = ReaderType::New(); + WriterType::Pointer writer = WriterType::New(); + reader->SetFileName(inputFileName); + writer->SetFileName(outputFileName); + reader->UpdateOutputInformation(); + + ReaderType::Pointer reader2 = ReaderType::New(); + WriterType::Pointer writer2 = WriterType::New(); + reader2->SetFileName(inputFileName2); + writer2->SetFileName(outputFileName2); + reader2->UpdateOutputInformation(); + + // Instantiating object + LuminanceToImageImageFilterType::Pointer filter = LuminanceToImageImageFilterType::New(); + + filter->SetInput(reader->GetOutput()); + writer->SetInput(filter->GetOutput()); + writer->Update(); + + RoiFilterType::Pointer roiFilter = RoiFilterType::New(); + roiFilter->SetStartX(1000); + roiFilter->SetStartY(1000); + roiFilter->SetSizeX(100); + roiFilter->SetSizeY(100); + roiFilter->SetInput(reader2->GetOutput()); + writer2->SetInput(roiFilter->GetOutput()); + writer2->Update(); + + return EXIT_SUCCESS; +} diff --git a/Testing/Code/Radiometry/otbLuminanceToImageImageFilterNew.cxx b/Testing/Code/Radiometry/otbLuminanceToImageImageFilterNew.cxx new file mode 100644 index 0000000000..d1509e5f25 --- /dev/null +++ b/Testing/Code/Radiometry/otbLuminanceToImageImageFilterNew.cxx @@ -0,0 +1,37 @@ +/*========================================================================= + + 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 "otbLuminanceToImageImageFilter.h" +#include "otbVectorImage.h" + +int otbLuminanceToImageImageFilterNew(int argc, char * argv[]) +{ + const unsigned int Dimension = 2; + typedef double PixelType; + typedef otb::VectorImage<PixelType, Dimension> InputImageType; + + typedef otb::LuminanceToImageImageFilter<InputImageType, InputImageType> LuminanceToImageImageFilterType; + + // Instantiating object + LuminanceToImageImageFilterType::Pointer filter = LuminanceToImageImageFilterType::New(); + + std::cout << filter << std::endl; + + return EXIT_SUCCESS; +} diff --git a/Testing/Code/Radiometry/otbRadiometryTests2.cxx b/Testing/Code/Radiometry/otbRadiometryTests2.cxx index eb184961d1..b80e68ce73 100644 --- a/Testing/Code/Radiometry/otbRadiometryTests2.cxx +++ b/Testing/Code/Radiometry/otbRadiometryTests2.cxx @@ -37,4 +37,13 @@ void RegisterTests() REGISTER_TEST(otbImageToReflectanceImageFilterNew); REGISTER_TEST(otbImageToReflectanceImageFilter); REGISTER_TEST(otbImageToReflectanceImageFilterAuto); + REGISTER_TEST(otbLuminanceToImageImageFilterNew); + REGISTER_TEST(otbLuminanceToImageImageFilter); + REGISTER_TEST(otbLuminanceToImageImageFilterAuto); + REGISTER_TEST(otbReflectanceToLuminanceImageFilterNew); + REGISTER_TEST(otbReflectanceToLuminanceImageFilter); + REGISTER_TEST(otbReflectanceToLuminanceImageFilterAuto); + REGISTER_TEST(otbReflectanceToImageImageFilterNew); + REGISTER_TEST(otbReflectanceToImageImageFilter); + REGISTER_TEST(otbReflectanceToImageImageFilterAuto); } diff --git a/Testing/Code/Radiometry/otbReflectanceToImageImageFilter.cxx b/Testing/Code/Radiometry/otbReflectanceToImageImageFilter.cxx new file mode 100644 index 0000000000..4eab8f1c44 --- /dev/null +++ b/Testing/Code/Radiometry/otbReflectanceToImageImageFilter.cxx @@ -0,0 +1,183 @@ +/*========================================================================= + + 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 "otbReflectanceToImageImageFilter.h" + +//#include "otbReflectanceToLuminanceImageFilter.h" +//#include "otbLuminanceToImageImageFilter.h" + +#include "otbVectorImage.h" +#include "otbImageFileReader.h" +#include "otbImageFileWriter.h" +#include "itkVariableLengthVector.h" + +int otbReflectanceToImageImageFilter(int argc, char * argv[]) +{ + const char * inputFileName = argv[1]; + const char * outputFileName = argv[2]; + const double angle = static_cast<double>(atof(argv[3])); + double flux = 0.; + int day = 1; + int month = 1; + + if (argc == 17) + { + flux = static_cast<double>(atof(argv[16])); + } + else + { + day = atoi(argv[16]); + month = atoi(argv[17]); + } + + 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::ReflectanceToImageImageFilter<InputImageType, OutputImageType> ReflectanceToImageImageFilterType; + typedef ReflectanceToImageImageFilterType::VectorType VectorType; + + ReaderType::Pointer reader = ReaderType::New(); + WriterType::Pointer writer = WriterType::New(); + reader->SetFileName(inputFileName); + writer->SetFileName(outputFileName); + reader->UpdateOutputInformation(); + + unsigned int nbOfComponent = reader->GetOutput()->GetNumberOfComponentsPerPixel(); + + VectorType alpha(nbOfComponent); + VectorType beta(nbOfComponent); + VectorType solarIllumination(nbOfComponent); + alpha.Fill(0); + beta.Fill(0); + solarIllumination.Fill(0); + + for (unsigned int i = 0; i < nbOfComponent; ++i) + { + alpha[i] = static_cast<double>(atof(argv[i + 4])); + beta[i] = static_cast<double>(atof(argv[i + 8])); + solarIllumination[i] = static_cast<double>(atof(argv[i + 12])); + } + + // Instantiating object + ReflectanceToImageImageFilterType::Pointer filter = ReflectanceToImageImageFilterType::New(); + + filter->SetAlpha(alpha); + filter->SetBeta(beta); + filter->SetZenithalSolarAngle(angle); + filter->SetSolarIllumination(solarIllumination); + + if (argc == 17) + { + filter->SetFluxNormalizationCoefficient(flux); + } + else + { + filter->SetDay(day); + filter->SetMonth(month); + } + + filter->SetInput(reader->GetOutput()); + writer->SetInput(filter->GetOutput()); + writer->Update(); + + return EXIT_SUCCESS; + + + /*const char * inputFileName = argv[1]; + const char * outputFileName = argv[2]; + const double angle = static_cast<double>(atof(argv[3])); + double flux = 0.; + int day = 1; + int month = 1; + + if (argc == 17) + { + flux = static_cast<double>(atof(argv[16])); + } + else + { + day = atoi(argv[16]); + month = atoi(argv[17]); + } + + 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::ReflectanceToLuminanceImageFilter<InputImageType, InputImageType> ReflectanceToLuminanceImageFilterType; +typedef otb::LuminanceToImageImageFilter<InputImageType, OutputImageType> LuminanceToImageImageFilterType; + + typedef ReflectanceToLuminanceImageFilterType::VectorType VectorType; + + ReaderType::Pointer reader = ReaderType::New(); + WriterType::Pointer writer = WriterType::New(); + reader->SetFileName(inputFileName); + writer->SetFileName(outputFileName); + reader->UpdateOutputInformation(); + + unsigned int nbOfComponent = reader->GetOutput()->GetNumberOfComponentsPerPixel(); + + VectorType alpha(nbOfComponent); + VectorType beta(nbOfComponent); + VectorType solarIllumination(nbOfComponent); + alpha.Fill(0); + beta.Fill(0); + solarIllumination.Fill(0); + + for (unsigned int i = 0; i < nbOfComponent; ++i) + { + alpha[i] = static_cast<double>(atof(argv[i + 4])); + beta[i] = static_cast<double>(atof(argv[i + 8])); + solarIllumination[i] = static_cast<double>(atof(argv[i + 12])); + } + + // Instantiating object + +ReflectanceToLuminanceImageFilterType::Pointer filter = ReflectanceToLuminanceImageFilterType::New(); +LuminanceToImageImageFilterType::Pointer filter2 = LuminanceToImageImageFilterType::New(); + + filter2->SetAlpha(alpha); + filter2->SetBeta(beta); + filter->SetZenithalSolarAngle(angle); + filter->SetSolarIllumination(solarIllumination); + + if (argc == 17) + { + filter->SetFluxNormalizationCoefficient(flux); + } + else + { + filter->SetDay(day); + filter->SetMonth(month); + } + + filter->SetInput(reader->GetOutput()); + filter2->SetInput(filter->GetOutput()); + writer->SetInput(filter2->GetOutput()); + writer->Update(); + + return EXIT_SUCCESS;*/ + +} diff --git a/Testing/Code/Radiometry/otbReflectanceToImageImageFilterAuto.cxx b/Testing/Code/Radiometry/otbReflectanceToImageImageFilterAuto.cxx new file mode 100644 index 0000000000..a56588aff0 --- /dev/null +++ b/Testing/Code/Radiometry/otbReflectanceToImageImageFilterAuto.cxx @@ -0,0 +1,72 @@ +/*========================================================================= + + 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 "otbReflectanceToImageImageFilter.h" +#include "otbVectorImage.h" +#include "otbImageFileReader.h" +#include "otbImageFileWriter.h" +#include "itkVariableLengthVector.h" +#include "otbMultiChannelExtractROI.h" + +int otbReflectanceToImageImageFilterAuto(int argc, char * argv[]) +{ + const char * inputFileName = argv[1]; + const char * inputFileName2 = argv[2]; + const char * outputFileName = argv[3]; + const char * outputFileName2 = argv[4]; + + 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::ReflectanceToImageImageFilter<InputImageType, OutputImageType> ReflectanceToImageImageFilterType; + typedef ReflectanceToImageImageFilterType::VectorType VectorType; + typedef otb::MultiChannelExtractROI<PixelType, PixelType> RoiFilterType; + + ReaderType::Pointer reader = ReaderType::New(); + WriterType::Pointer writer = WriterType::New(); + reader->SetFileName(inputFileName); + writer->SetFileName(outputFileName); + reader->UpdateOutputInformation(); + + ReaderType::Pointer reader2 = ReaderType::New(); + WriterType::Pointer writer2 = WriterType::New(); + reader2->SetFileName(inputFileName2); + writer2->SetFileName(outputFileName2); + reader2->UpdateOutputInformation(); + + // Instantiating object + ReflectanceToImageImageFilterType::Pointer filter = ReflectanceToImageImageFilterType::New(); + filter->SetInput(reader->GetOutput()); + writer->SetInput(filter->GetOutput()); + writer->Update(); + + RoiFilterType::Pointer roiFilter = RoiFilterType::New(); + roiFilter->SetStartX(1000); + roiFilter->SetStartY(1000); + roiFilter->SetSizeX(100); + roiFilter->SetSizeY(100); + roiFilter->SetInput(reader2->GetOutput()); + writer2->SetInput(roiFilter->GetOutput()); + writer2->Update(); + + return EXIT_SUCCESS; +} diff --git a/Testing/Code/Radiometry/otbReflectanceToImageImageFilterNew.cxx b/Testing/Code/Radiometry/otbReflectanceToImageImageFilterNew.cxx new file mode 100644 index 0000000000..ec14c72636 --- /dev/null +++ b/Testing/Code/Radiometry/otbReflectanceToImageImageFilterNew.cxx @@ -0,0 +1,37 @@ +/*========================================================================= + + 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 "otbReflectanceToImageImageFilter.h" +#include "otbVectorImage.h" + +int otbReflectanceToImageImageFilterNew(int argc, char * argv[]) +{ + const unsigned int Dimension = 2; + typedef double PixelType; + typedef otb::VectorImage<PixelType, Dimension> InputImageType; + + typedef otb::ReflectanceToImageImageFilter<InputImageType, InputImageType> ReflectanceToImageImageFilterType; + + // Instantiating object + ReflectanceToImageImageFilterType::Pointer filter = ReflectanceToImageImageFilterType::New(); + + std::cout << filter << std::endl; + + return EXIT_SUCCESS; +} diff --git a/Testing/Code/Radiometry/otbReflectanceToLuminanceImageFilter.cxx b/Testing/Code/Radiometry/otbReflectanceToLuminanceImageFilter.cxx new file mode 100644 index 0000000000..e0088841f8 --- /dev/null +++ b/Testing/Code/Radiometry/otbReflectanceToLuminanceImageFilter.cxx @@ -0,0 +1,90 @@ +/*========================================================================= + + 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 "otbReflectanceToLuminanceImageFilter.h" +#include "otbVectorImage.h" +#include "otbImageFileReader.h" +#include "otbImageFileWriter.h" +#include "itkVariableLengthVector.h" + +int otbReflectanceToLuminanceImageFilter(int argc, char * argv[]) +{ + const char * inputFileName = argv[1]; + const char * outputFileName = argv[2]; + const double angle = static_cast<double>(atof(argv[3])); + double flux = 0.; + int day = 1; + int month = 1; + + if (argc == 9) + { + flux = static_cast<double>(atof(argv[8])); + } + else + { + day = atoi(argv[8]); + month = atoi(argv[9]); + } + + 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::ReflectanceToLuminanceImageFilter<InputImageType, OutputImageType> ReflectanceToLuminanceImageFilterType; + typedef ReflectanceToLuminanceImageFilterType::VectorType VectorType; + + ReaderType::Pointer reader = ReaderType::New(); + WriterType::Pointer writer = WriterType::New(); + reader->SetFileName(inputFileName); + writer->SetFileName(outputFileName); + reader->UpdateOutputInformation(); + + unsigned int nbOfComponent = reader->GetOutput()->GetNumberOfComponentsPerPixel(); + // WARNING : 1 ELEMENT FOR EACH PIXEL IMAGE COMPONENT + VectorType solarIllumination(nbOfComponent); + solarIllumination.Fill(0); + solarIllumination[0] = static_cast<double>(atof(argv[4])); + solarIllumination[1] = static_cast<double>(atof(argv[5])); + solarIllumination[2] = static_cast<double>(atof(argv[6])); + solarIllumination[3] = static_cast<double>(atof(argv[7])); + + // Instantiating object + ReflectanceToLuminanceImageFilterType::Pointer filter = ReflectanceToLuminanceImageFilterType::New(); + + filter->SetZenithalSolarAngle(angle); + filter->SetSolarIllumination(solarIllumination); + filter->SetUseClamp(false); + if (argc == 9) + { + filter->SetFluxNormalizationCoefficient(flux); + } + else + { + filter->SetDay(day); + filter->SetMonth(month); + } + + filter->SetInput(reader->GetOutput()); + writer->SetInput(filter->GetOutput()); + writer->Update(); + + return EXIT_SUCCESS; +} diff --git a/Testing/Code/Radiometry/otbReflectanceToLuminanceImageFilterAuto.cxx b/Testing/Code/Radiometry/otbReflectanceToLuminanceImageFilterAuto.cxx new file mode 100644 index 0000000000..16f42fef03 --- /dev/null +++ b/Testing/Code/Radiometry/otbReflectanceToLuminanceImageFilterAuto.cxx @@ -0,0 +1,58 @@ +/*========================================================================= + + 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 "otbReflectanceToLuminanceImageFilter.h" +#include "otbVectorImage.h" +#include "otbImageFileReader.h" +#include "otbImageFileWriter.h" +#include "itkVariableLengthVector.h" +#include "otbMultiChannelExtractROI.h" + +int otbReflectanceToLuminanceImageFilterAuto(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::ReflectanceToLuminanceImageFilter<OutputImageType, + OutputImageType> ReflectanceToLuminanceImageFilterType; + typedef ReflectanceToLuminanceImageFilterType::VectorType VectorType; + + ReaderType::Pointer reader = ReaderType::New(); + WriterType::Pointer writer = WriterType::New(); + reader->SetFileName(inputFileName); + writer->SetFileName(outputFileName); + reader->UpdateOutputInformation(); + + // Instantiating object + ReflectanceToLuminanceImageFilterType::Pointer filterToLuminance = ReflectanceToLuminanceImageFilterType::New(); + + filterToLuminance->SetInput(reader->GetOutput()); + writer->SetInput(filterToLuminance->GetOutput()); + writer->Update(); + + return EXIT_SUCCESS; + +} diff --git a/Testing/Code/Radiometry/otbReflectanceToLuminanceImageFilterNew.cxx b/Testing/Code/Radiometry/otbReflectanceToLuminanceImageFilterNew.cxx new file mode 100644 index 0000000000..6e86cc9ba3 --- /dev/null +++ b/Testing/Code/Radiometry/otbReflectanceToLuminanceImageFilterNew.cxx @@ -0,0 +1,37 @@ +/*========================================================================= + + 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 "otbReflectanceToLuminanceImageFilter.h" +#include "otbVectorImage.h" + +int otbReflectanceToLuminanceImageFilterNew(int argc, char * argv[]) +{ + const unsigned int Dimension = 2; + typedef double PixelType; + typedef otb::VectorImage<PixelType, Dimension> InputImageType; + + typedef otb::ReflectanceToLuminanceImageFilter<InputImageType, InputImageType> ReflectanceToLuminanceImageFilterType; + + // Instantiating object + ReflectanceToLuminanceImageFilterType::Pointer filter = ReflectanceToLuminanceImageFilterType::New(); + + std::cout << filter << std::endl; + + return EXIT_SUCCESS; +} -- GitLab