From 0d6569b25c9dbee5d50656f30fb5bb787fb84d0d Mon Sep 17 00:00:00 2001
From: Otmane Lahlou <otmane.lahlou@c-s.fr>
Date: Thu, 13 Oct 2011 11:09:17 +0200
Subject: [PATCH] ENH: port orthorectification to new application api
---
.../Projections/otbOrthoRectification.cxx | 906 +++++++-----------
1 file changed, 326 insertions(+), 580 deletions(-)
diff --git a/Applications/Projections/otbOrthoRectification.cxx b/Applications/Projections/otbOrthoRectification.cxx
index 7763395df8..522f8b7981 100644
--- a/Applications/Projections/otbOrthoRectification.cxx
+++ b/Applications/Projections/otbOrthoRectification.cxx
@@ -15,646 +15,392 @@
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
-#include "otbOrthoRectification.h"
-
-#include "otbVectorImage.h"
-
-#include "otbImageFileReader.h"
-#include "otbStreamingImageFileWriter.h"
-#include "otbStandardWriterWatcher.h"
-
-#include "otbOrthoRectificationFilter.h"
+#include "otbWrapperApplication.h"
+#include "otbWrapperApplicationFactory.h"
+#include "otbWrapperNumericalParameter.h"
+#include "otbGenericRSResampleImageFilter.h"
+#include "otbImageToGenericRSOutputParameters.h"
#include "otbMapProjections.h"
-#include "itkExceptionObject.h"
-#include "otbMacro.h"
#include "itkLinearInterpolateImageFunction.h"
#include "otbBCOInterpolateImageFunction.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
-#include "otbImageToGenericRSOutputParameters.h"
+
+#include "otbMacro.h"
namespace otb
{
-template<typename TMapProjection>
-int generic_main(otb::ApplicationOptionsResult* parseResult,
- TMapProjection* mapProjection)
+enum
{
- try
- {
- typedef otb::VectorImage<unsigned short int, 2> ImageType;
- typedef otb::ImageFileReader<ImageType> ReaderType;
- typedef otb::StreamingImageFileWriter<ImageType> WriterType;
-
- typedef TMapProjection MapProjectionType;
- typedef otb::OrthoRectificationFilter< ImageType, ImageType, MapProjectionType > OrthorectifFilterType;
-
- typedef itk::LinearInterpolateImageFunction<ImageType, double> LinearInterpolationType;
- typedef itk::NearestNeighborInterpolateImageFunction<ImageType, double> NearestNeighborInterpolationType;
- typedef otb::BCOInterpolateImageFunction<ImageType> BCOInterpolationType;
-
- typedef otb::PipelineMemoryPrintCalculator MemoryCalculatorType;
- typedef itk::ExtractImageFilter<ImageType, ImageType> ExtractFilterType;
-
- // Read input image information
- ReaderType::Pointer reader=ReaderType::New();
- reader->SetFileName(parseResult->GetInputImage().c_str());
- reader->GenerateOutputInformation();
- ImageType::SizeType lsize = reader->GetOutput()->GetLargestPossibleRegion().GetSize();
-
- // Orthorectification filter
- typename OrthorectifFilterType::Pointer orthoFilter = OrthorectifFilterType::New();
- orthoFilter->SetInput(reader->GetOutput());
-
- // If activated, generate RPC model
- if(parseResult->IsOptionPresent("RPC"))
- {
- orthoFilter->EstimateInputRpcModelOn();
- orthoFilter->SetInputRpcGridSize(parseResult->GetParameterUInt("RPC"));
- }
-
- // Compute the output parameters
- typedef otb::ImageToGenericRSOutputParameters<ImageType> OutputParametersEstimatorType;
- typename OutputParametersEstimatorType::Pointer genericRSEstimator = OutputParametersEstimatorType::New();
- genericRSEstimator->SetInput(reader->GetOutput());
- genericRSEstimator->SetOutputProjectionRef(mapProjection->GetWkt());
-
- // Set up output image informations
-
- ImageType::SpacingType spacing;
- if(parseResult->IsOptionPresent("OutputSpacing"))
- {
- spacing[0]= parseResult->GetParameterDouble("OutputSpacing", 0);
- spacing[1]= parseResult->GetParameterDouble("OutputSpacing", 1);
- }
- else
- {
- typedef otb::ImageMetadataInterfaceBase ImageMetadataInterfaceType;
- ImageMetadataInterfaceType::Pointer metadataInterface = ImageMetadataInterfaceFactory::CreateIMI(
- reader->GetOutput()->GetMetaDataDictionary());
- double isotropicSpacing = std::max(metadataInterface->GetXPixelSpacing(), metadataInterface->GetYPixelSpacing());
+ Map_Utm,
+ Map_Lambert2,
+ Map_Epsg
+};
- spacing[0] = isotropicSpacing;
- spacing[1] = -isotropicSpacing;
- }
- genericRSEstimator->ForceSpacingTo(spacing);
-
- otbMsgDevMacro(<< "Output image spacing: " << spacing);
-
- genericRSEstimator->Compute();
- // set the start index
- ImageType::IndexType start;
- start[0]=0;
- start[1]=0;
-
- //Update now the orthoFilter with output image parameters
- ImageType::PixelType defaultValue;
- itk::PixelBuilder<ImageType::PixelType>::Zero(defaultValue,
- reader->GetOutput()->GetNumberOfComponentsPerPixel());
-
- orthoFilter->SetInput(reader->GetOutput());
-
- ImageType::PointType origin;
- if(parseResult->IsOptionPresent("UpperLeft"))
- {
- origin[0]= parseResult->GetParameterDouble("UpperLeft", 0);
- origin[1]= parseResult->GetParameterDouble("UpperLeft", 1);
- }
- else
- {
- origin = genericRSEstimator->GetOutputOrigin();
- }
-
- otbMsgDevMacro(<< "Output image origin: " << origin);
+enum
+{
+ Interpolator_Linear,
+ Interpolator_NNeighbor,
+ Interpolator_BCO
+};
- ImageType::SizeType size;
- if(parseResult->IsOptionPresent("OutputSize"))
- {
- size[0]= parseResult->GetParameterULong("OutputSize", 0);
- size[1]= parseResult->GetParameterULong("OutputSize", 1);
- genericRSEstimator->ForceSizeTo(size);
- }
- else
- {
- size = genericRSEstimator->GetOutputSize();
- }
+namespace Wrapper
+{
- otbMsgDevMacro(<< "Output image size: " << size);
+class OrthoRectification : public Application
+{
+public:
+/** Standard class typedefs. */
+ typedef OrthoRectification Self;
+ typedef Application Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Standard macro */
+ itkNewMacro(Self);
+
+ itkTypeMacro(OrthoRectification, otb::Application);
+
+ /** Generic Remote Sensor Resampler */
+ typedef otb::GenericRSResampleImageFilter<FloatVectorImageType,
+ FloatVectorImageType> ResampleFilterType;
+
+ /** Interpolators typedefs*/
+ typedef itk::LinearInterpolateImageFunction<FloatVectorImageType,
+ double> LinearInterpolationType;
+ typedef itk::NearestNeighborInterpolateImageFunction<FloatVectorImageType,
+ double> NearestNeighborInterpolationType;
+ typedef otb::BCOInterpolateImageFunction<FloatVectorImageType> BCOInterpolationType;
+
+ typedef otb::PipelineMemoryPrintCalculator MemoryCalculatorType;
+
+private:
+ OrthoRectification()
+ {
+ SetName("OrthoRectification");
+ SetDescription("Using available image metadata to determine the sensor model, \ncomputes a cartographic projection of the image");
+ }
- orthoFilter->SetOutputOrigin(origin);
- orthoFilter->SetOutputSpacing(genericRSEstimator->GetOutputSpacing());
- orthoFilter->SetOutputSize(size);
- orthoFilter->SetOutputStartIndex(start);
- orthoFilter->SetEdgePaddingValue(defaultValue);
+ void DoCreateParameters()
+ {
+ // Set the parameters
+ AddParameter(ParameterType_InputImage, "in", "Input Image");
+
+ AddParameter(ParameterType_OutputImage, "out", "Output Image");
+ SetParameterDescription("out","Projected image");
+
+
+ // Add the output paramters in a group
+ AddParameter(ParameterType_Group, "outputs", "Output Parameters Group");
+ MandatoryOn("outputs");
+
+ // UserDefined values
+ AddParameter(ParameterType_Empty, "outputs.uservalues", "User defined Parameters");
+
+ // Upper left point coordinates
+ AddParameter(ParameterType_Float, "outputs.ulx", "Upper Left X");
+ SetParameterDescription("outputs.ulx","Cartographic X coordinate of upper left corner");
+ //MandatoryOn("outputs.ulx");
+
+ AddParameter(ParameterType_Float, "outputs.uly", "Upper Left Y");
+ SetParameterDescription("outputs.uly","Cartographic Y coordinate of upper left corner");
+ //MandatoryOn("outputs.uly");
+
+ // Size of the output image
+ AddParameter(ParameterType_Int, "outputs.sizex", "Size X");
+ SetParameterDescription("outputs.sizex","Size of projected image along X");
+ MandatoryOn("outputs.sizex");
+
+ AddParameter(ParameterType_Int, "outputs.sizey", "Size Y");
+ SetParameterDescription("outputs.sizey","Size of projected image along Y");
+ MandatoryOn("outputs.sizey");
+
+ // Spacing of the output image
+ AddParameter(ParameterType_Float, "outputs.spacingx", "Pixel Size X");
+ SetParameterDescription("outputs.spacingx","Size of each pixel along X axis");
+ MandatoryOn("outputs.spacingx");
+
+ AddParameter(ParameterType_Float, "outputs.spacingy", "Pixel Size Y");
+ SetParameterDescription("outputs.spacingy","Size of each pixel along Y axis");
+ MandatoryOn("outputs.spacingy");
+
+ // DEM
+ AddParameter(ParameterType_Directory, "dem", "DEM directory");
+ MandatoryOff("dem");
+
+ // Estimate a RPC model (for spot image for instance)
+ AddParameter(ParameterType_Group, "rpc", "Estimate RPC model");
+ AddParameter(ParameterType_Int, "rpc.ncp", "Nb control Points");
+ SetParameterDescription("rpc","Activate RPC sensor model estimation. Parameter is the number of control points per axis");
+ MandatoryOff("rpc");
+ MandatoryOff("rpc.ncp");
+
+ // Interpolators
+ AddParameter(ParameterType_Choice, "interpolator", "Interpolator");
+ AddChoice("interpolator.linear", "Linear");
+ AddChoice("interpolator.nn", "Neareast Neighbor");
+ AddChoice("interpolator.bco", "BCO");
+ AddParameter(ParameterType_Radius, "interpolator.bco.radius", "Radius");
+ SetParameterInt("interpolator.bco.radius", 2);
+
+ // Built the Output Map Projection
+ AddParameter(ParameterType_Choice, "map", "Map Projection");
+
+ AddChoice("map.utm", "UTM"); // OK
+ AddParameter(ParameterType_Int, "map.utm.zone", "Zone number");
+ AddParameter(ParameterType_Empty, "map.utm.hem", "Hemisphere North");
+ MandatoryOff("map.utm.zone");
+ MandatoryOff("map.utm.hem");
+
+
+ AddChoice("map.lambert2", "Lambert II Etendu"); // OK
+
+ AddChoice("map.epsg","EPSG");
+ AddParameter(ParameterType_Int, "map.epsg.code", "EPSG Code");
+ SetParameterInt("map.epsg.code",32631);
+ SetParameterString("map", "epsg");
+
+ // descriptor->AddOption("LocMapSpacing",
+ // "Generate a coarser deformation field with the given spacing.","lmSpacing",
+ // 1, false, otb::ApplicationDescriptor::Real);
+ }
- // Configure DEM directory
- if(parseResult->IsOptionPresent("DEMDirectory"))
- {
- orthoFilter->SetDEMDirectory(parseResult->GetParameterString("DEMDirectory", 0));
- }
- else
+ void DoUpdateParameters()
+ {
+ if (HasValue("in"))
{
- if ( otb::ConfigurationFile::GetInstance()->IsValid() )
- {
- orthoFilter->SetDEMDirectory(otb::ConfigurationFile::GetInstance()->GetDEMDirectory());
- }
+ // input image
+ FloatVectorImageType* inImage = GetParameterImage("in");
+
+ // Get the output projection Ref
+ this->UpdateOutputProjectionRef();
+
+
+
+ // Compute the output image spacing and size
+ typedef otb::ImageToGenericRSOutputParameters<FloatVectorImageType> OutputParametersEstimatorType;
+ OutputParametersEstimatorType::Pointer genericRSEstimator = OutputParametersEstimatorType::New();
+ genericRSEstimator->SetInput(inImage);
+ genericRSEstimator->SetOutputProjectionRef(m_OutputProjectionRef);
+ genericRSEstimator->Compute();
+
+ // TODO : Add a ParameterGroup with the output Parameters
+ // Add an option user defined param or automatic params
+
+ dynamic_cast< NumericalParameter<float> * >(GetParameterByKey("outputs.spacingy"))->SetMinimumValue(-50000.0);
+ dynamic_cast< NumericalParameter<float> * >(GetParameterByKey("outputs.spacingy"))->SetMaximumValue(10000.0);
+
+ // Fill the Gui with the computed parameters
+ if (!HasUserValue("outputs.ulx"))
+ SetParameterFloat("outputs.ulx", genericRSEstimator->GetOutputOrigin()[0]);
+
+ if (!HasUserValue("outputs.uly"))
+ SetParameterFloat("outputs.uly", genericRSEstimator->GetOutputOrigin()[1]);
+
+ if (!HasUserValue("outputs.sizex"))
+ SetParameterInt("outputs.sizex", genericRSEstimator->GetOutputSize()[0]);
+
+ if (!HasUserValue("outputs.sizey"))
+ SetParameterInt("outputs.sizey", genericRSEstimator->GetOutputSize()[1]);
+
+ if (!HasUserValue("outputs.spacingx"))
+ SetParameterFloat("outputs.spacingx", genericRSEstimator->GetOutputSpacing()[0]);
+
+ if (!HasUserValue("outputs.spacingy"))
+ SetParameterFloat("outputs.spacingy", genericRSEstimator->GetOutputSpacing()[1]);
+
+// // Mode : user defined
+// if (!IsParameterEnabled("outputs.uservalues"))
+// {
+// // Enable UL corners widget
+// DisableParameter("outputs.ulx");
+// DisableParameter("outputs.uly");
+// // Force size or spacing to user defined params
+// if (HasUserValue("outputs.sizex") || HasUserValue("outputs.sizey"))
+// {
+// ResampleFilterType::SizeType size;
+// size[0] = GetParameterInt("outputs.sizex");
+// size[1] = GetParameterInt("outputs.sizey");
+// genericRSEstimator->ForceSizeTo(size);
+// genericRSEstimator->Compute();
+
+// std::cout <<"Size Forced to "<<size << " --> implies Sapcing : "
+// <<genericRSEstimator->GetOutputSpacing() << std::endl;
+
+// // Set the processed spacing relative to this forced size
+// SetParameterFloat("outputs.spacingx", genericRSEstimator->GetOutputSpacing()[0]);
+// SetParameterFloat("outputs.spacingy", genericRSEstimator->GetOutputSpacing()[1]);
+// }
+
+// if (HasUserValue("outputs.spacingy") || HasUserValue("outputs.spacingx"))
+// {
+// ResampleFilterType::SpacingType spacing;
+// spacing[0] = GetParameterFloat("outputs.spacingx");
+// spacing[1] = GetParameterFloat("outputs.spacingy");
+
+// genericRSEstimator->ForceSpacingTo(spacing);
+// genericRSEstimator->Compute();
+
+// std::cout <<"Spacing Forced to "<< spacing<< " --> implies Size : "
+// <<genericRSEstimator->GetOutputSize() << std::endl;
+
+// // Set the processed size relative to this forced spacing
+// SetParameterInt("outputs.sizex", genericRSEstimator->GetOutputSize()[0]);
+// SetParameterInt("outputs.sizey", genericRSEstimator->GetOutputSize()[1]);
+// }
+// }
+// else
+// {
+// // Enable UL corners widget
+// EnableParameter("outputs.ulx");
+// EnableParameter("outputs.uly");
+
+// }
}
+ }
- // Set the output map projection
- orthoFilter->SetMapProjection(mapProjection);
-
- if (parseResult->IsOptionPresent("LocMapSpacing"))
+ void UpdateOutputProjectionRef()
+ {
+ switch ( GetParameterInt("map") )
{
- double defScalarSpacing = parseResult->GetParameterFloat("LocMapSpacing");
- ImageType::SpacingType defSpacing;
- defSpacing[0] = defScalarSpacing;
- defSpacing[1] = -defScalarSpacing;
-
- orthoFilter->SetDeformationFieldSpacing(defSpacing);
- }
- else
+ case Map_Utm:
{
- // By default, generate a 10 times coarser deformation field
- ImageType::SpacingType defSpacing;
- defSpacing[0] = 10*spacing[0];
- defSpacing[1] = 10*spacing[1];
- orthoFilter->SetDeformationFieldSpacing(defSpacing);
+ typedef UtmInverseProjection UtmProjectionType;
+ UtmProjectionType::Pointer utmProjection = UtmProjectionType::New();
+
+ // Set the zone
+ utmProjection->SetZone(GetParameterInt("map.utm.zone"));
+
+ // Set the hem
+ std::string hem = "N";
+ if (!IsParameterEnabled("map.utm.hem"))
+ hem = "S";
+ utmProjection->SetHemisphere(hem[0]);
+
+ // Get the projection ref
+ m_OutputProjectionRef = utmProjection->GetWkt();
}
-
- // Set the interpolator type
- if (parseResult->IsOptionPresent("InterpolatorType"))
+ break;
+ case Map_Lambert2:
{
- int nbInterpParams = parseResult->GetNumberOfParameters("InterpolatorType") - 1;
- std::string typeInterpolator = parseResult->GetParameterString("InterpolatorType");
-
- if (typeInterpolator == "BCO")
- {
- BCOInterpolationType::Pointer interpolator = BCOInterpolationType::New();
- if( nbInterpParams == 1)
- interpolator->SetRadius(parseResult->GetParameterUInt("InterpolatorType", 1));
- orthoFilter->SetInterpolator(interpolator);
- }
- else if (typeInterpolator == "NEARESTNEIGHBOR" && nbInterpParams == 0)
- {
- NearestNeighborInterpolationType::Pointer interpolator = NearestNeighborInterpolationType::New();
- orthoFilter->SetInterpolator(interpolator);
- }
- else if (typeInterpolator == "LINEAR" && nbInterpParams == 0)
- {
- LinearInterpolationType::Pointer interpolator = LinearInterpolationType::New();
- orthoFilter->SetInterpolator(interpolator);
- }
- else
- {
- itkGenericExceptionMacro(<< "Interpolator type not recognized, choose one with (parameters) : BCO(0/1), NEARESTNEIGHBOR(0), LINEAR(0)");
- }
+ typedef Lambert2EtenduForwardProjection Lambert2ProjectionType;
+ Lambert2ProjectionType::Pointer lambert2Projection = Lambert2ProjectionType::New();
+ m_OutputProjectionRef = lambert2Projection->GetWkt();
}
- else
+ break;
+ case Map_Epsg:
{
- // Default is BCO
- BCOInterpolationType::Pointer interpolator = BCOInterpolationType::New();
- orthoFilter->SetInterpolator(interpolator);
+ m_OutputProjectionRef = otb::GeoInformationConversion::ToWKT(GetParameterInt("map.epsg.code"));
}
-
- //Update informations from orthorectification filter
- orthoFilter->GetOutput()->UpdateOutputInformation();
-
- //Instantiate the writer
- WriterType::Pointer writer = WriterType::New();
- writer->SetFileName(parseResult->GetOutputImage());
- writer->SetInput(orthoFilter->GetOutput());
-
- unsigned int ram = 256;
- if (parseResult->IsOptionPresent("AvailableMemory"))
- {
- ram = parseResult->GetParameterUInt("AvailableMemory");
+ break;
}
- writer->SetAutomaticTiledStreaming(ram);
-
- otb::StandardWriterWatcher watcher(writer, orthoFilter,"Orthorectification");
- writer->Update();
- }
- catch ( itk::ExceptionObject & err )
- {
- std::cout << "Exception itk::ExceptionObject raised !" << std::endl;
- std::cout << err << std::endl;
- return EXIT_FAILURE;
- }
- catch ( std::bad_alloc & err )
- {
- std::cout << "Exception bad_alloc : "<<(char*)err.what()<< std::endl;
- return EXIT_FAILURE;
}
- catch ( ... )
- {
- std::cout << "Unknown exception raised !" << std::endl;
- return EXIT_FAILURE;
- }
- return EXIT_SUCCESS;
-
-}
-
-int OrthoRectification::Describe(ApplicationDescriptor* descriptor)
-{
- descriptor->SetName("Orthorectification");
- descriptor->SetDescription("Using available image metadata to determine the sensor model, computes a cartographic projection of the image");
- descriptor->AddInputImage();
- descriptor->AddOutputImage();
- descriptor->AddOption("UpperLeft","Cartographic X/Y coordinate of upper left corner ","ul", 2, false, otb::ApplicationDescriptor::Real);
- descriptor->AddOption("OutputSize","Size of result image in X/Y","size", 2, false, otb::ApplicationDescriptor::Integer);
- descriptor->AddOption("OutputSpacing","Spacing resolution in meters on X/Y Axis","spacing", 2, false, otb::ApplicationDescriptor::Real);
- descriptor->AddOption("DEMDirectory","Directory where to find the DEM tiles","dem", 1, false, otb::ApplicationDescriptor::DirectoryName);
- descriptor->AddOptionNParams("MapProjectionType",
- "Type (UTM/LAMBERT/LAMBERT2/LAMBERT93/SINUS/ECKERT4/TRANSMERCATOR/MOLLWEID/EPSG) and parameters of map projection used. If EPSG is used, their is only one parameter corresponding to the EPSG code (see http://spatialreference.org).",
- "mapProj", false, otb::ApplicationDescriptor::String);
- descriptor->AddOption("RPC","Activate RPC sensor model estimation. Parameter is the number of control points per axis.","rpc", 1, false, otb::ApplicationDescriptor::Integer);
- descriptor->AddOption("LocMapSpacing","Generate a coarser deformation field with the given spacing.","lmSpacing", 1, false, otb::ApplicationDescriptor::Real);
- descriptor->AddOptionNParams("InterpolatorType",
- "Type LINEAR/BCO/NEARESTNEIGHBOR (optional, BCO by default)","interp", false, otb::ApplicationDescriptor::String);
- descriptor->AddOption("AvailableMemory","Set the maximum of available memory for the pipeline execution in mega bytes (optional, 256 by default)","ram", 1, false, otb::ApplicationDescriptor::Integer);
-
- return EXIT_SUCCESS;
-}
-
-int epsg_main(otb::ApplicationOptionsResult* parseResult)
-{
- try
+ void DoExecute()
{
- typedef otb::VectorImage<unsigned short int, 2> ImageType;
- typedef otb::ImageFileReader<ImageType> ReaderType;
- typedef otb::StreamingImageFileWriter<ImageType> WriterType;
-
- typedef otb::GenericRSResampleImageFilter< ImageType, ImageType > OrthorectifFilterType;
-
- typedef itk::LinearInterpolateImageFunction<ImageType, double> LinearInterpolationType;
- typedef itk::NearestNeighborInterpolateImageFunction<ImageType, double> NearestNeighborInterpolationType;
- typedef otb::BCOInterpolateImageFunction<ImageType> BCOInterpolationType;
-
- typedef otb::PipelineMemoryPrintCalculator MemoryCalculatorType;
- typedef itk::ExtractImageFilter<ImageType, ImageType> ExtractFilterType;
-
- // Read input image information
- ReaderType::Pointer reader=ReaderType::New();
- reader->SetFileName(parseResult->GetInputImage().c_str());
- reader->GenerateOutputInformation();
- ImageType::SizeType lsize = reader->GetOutput()->GetLargestPossibleRegion().GetSize();
-
- // Orthorectification filter
- OrthorectifFilterType::Pointer orthoFilter = OrthorectifFilterType::New();
- orthoFilter->SetInput(reader->GetOutput());
-
- // If activated, generate RPC model
- if(parseResult->IsOptionPresent("RPC"))
- {
- orthoFilter->EstimateInputRpcModelOn();
- orthoFilter->SetInputRpcGridSize(parseResult->GetParameterUInt("RPC"));
- }
-
- std::string outputProjectionRef = otb::GeoInformationConversion::ToWKT(parseResult->GetParameterInt("MapProjectionType", 1));
-
- otbMsgDevMacro(<<"Output projection ref: "<<outputProjectionRef);
-
- // Compute the output parameters
- typedef otb::ImageToGenericRSOutputParameters<ImageType> OutputParametersEstimatorType;
- OutputParametersEstimatorType::Pointer genericRSEstimator = OutputParametersEstimatorType::New();
- genericRSEstimator->SetInput(reader->GetOutput());
- genericRSEstimator->SetOutputProjectionRef(outputProjectionRef);
-
- // Set up output image informations
-
- ImageType::SpacingType spacing;
- if(parseResult->IsOptionPresent("OutputSpacing"))
+ GetLogger()->Debug("Entering DoExecute");
+
+ // Get the input image
+ FloatVectorImageType* inImage = GetParameterImage("in");
+
+ // Resampler Instanciation
+ m_ResampleFilter = ResampleFilterType::New();
+ m_ResampleFilter->SetInput(inImage);
+
+ // Set the output projection Ref
+ m_ResampleFilter->SetInputProjectionRef(inImage->GetProjectionRef());
+ m_ResampleFilter->SetInputKeywordList(inImage->GetImageKeywordlist());
+ m_ResampleFilter->SetOutputProjectionRef(m_OutputProjectionRef);
+
+ // Get Interpolator
+ switch ( GetParameterInt("interpolator") )
{
- spacing[0]= parseResult->GetParameterDouble("OutputSpacing", 0);
- spacing[1]= parseResult->GetParameterDouble("OutputSpacing", 1);
- }
- else
+ case Interpolator_Linear:
{
- typedef otb::ImageMetadataInterfaceBase ImageMetadataInterfaceType;
- ImageMetadataInterfaceType::Pointer metadataInterface = ImageMetadataInterfaceFactory::CreateIMI(
- reader->GetOutput()->GetMetaDataDictionary());
- double isotropicSpacing = std::max(metadataInterface->GetXPixelSpacing(), metadataInterface->GetYPixelSpacing());
-
- spacing[0] = isotropicSpacing;
- spacing[1] = -isotropicSpacing;
+ typedef itk::LinearInterpolateImageFunction<FloatVectorImageType, double> LinearInterpolationType;
+ LinearInterpolationType::Pointer interpolator = LinearInterpolationType::New();
+ m_ResampleFilter->SetInterpolator(interpolator);
}
- genericRSEstimator->ForceSpacingTo(spacing);
-
- otbMsgDevMacro(<< "Output image spacing: " << spacing);
-
- genericRSEstimator->Compute();
- // set the start index
- ImageType::IndexType start;
- start[0]=0;
- start[1]=0;
-
- //Update now the orthoFilter with output image parameters
- ImageType::PixelType defaultValue;
- itk::PixelBuilder<ImageType::PixelType>::Zero(defaultValue,
- reader->GetOutput()->GetNumberOfComponentsPerPixel());
-
- orthoFilter->SetInput(reader->GetOutput());
-
- ImageType::PointType origin;
- if(parseResult->IsOptionPresent("UpperLeft"))
+ break;
+ case Interpolator_NNeighbor:
{
- origin[0]= parseResult->GetParameterDouble("UpperLeft", 0);
- origin[1]= parseResult->GetParameterDouble("UpperLeft", 1);
+ typedef itk::NearestNeighborInterpolateImageFunction<FloatVectorImageType, double> NearestNeighborInterpolationType;
+ NearestNeighborInterpolationType::Pointer interpolator = NearestNeighborInterpolationType::New();
+ m_ResampleFilter->SetInterpolator(interpolator);
}
- else
- {
- origin = genericRSEstimator->GetOutputOrigin();
- }
-
- otbMsgDevMacro(<< "Output image origin: " << origin);
-
- ImageType::SizeType size;
- if(parseResult->IsOptionPresent("OutputSize"))
+ break;
+ case Interpolator_BCO:
{
- size[0]= parseResult->GetParameterDouble("OutputSize", 0);
- size[1]= parseResult->GetParameterDouble("OutputSize", 1);
- genericRSEstimator->ForceSizeTo(size);
+ typedef otb::BCOInterpolateImageFunction<FloatVectorImageType> BCOInterpolationType;
+ BCOInterpolationType::Pointer interpolator = BCOInterpolationType::New();
+ interpolator->SetRadius(GetParameterInt("interpolator.bco.radius"));
+ m_ResampleFilter->SetInterpolator(interpolator);
}
- else
- {
- size = genericRSEstimator->GetOutputSize();
+ break;
}
- otbMsgDevMacro(<< "Output image size: " << size);
-
- orthoFilter->SetOutputOrigin(origin);
- orthoFilter->SetOutputSpacing(genericRSEstimator->GetOutputSpacing());
- orthoFilter->SetOutputSize(size);
- orthoFilter->SetOutputStartIndex(start);
- orthoFilter->SetEdgePaddingValue(defaultValue);
-
- // Configure DEM directory
- if(parseResult->IsOptionPresent("DEMDirectory"))
+ // DEM
+ if (IsParameterEnabled("dem") && HasValue("dem"))
{
- orthoFilter->SetDEMDirectory(parseResult->GetParameterString("DEMDirectory", 0));
+ m_ResampleFilter->SetDEMDirectory(GetParameterString("dem"));
}
else
{
if ( otb::ConfigurationFile::GetInstance()->IsValid() )
{
- orthoFilter->SetDEMDirectory(otb::ConfigurationFile::GetInstance()->GetDEMDirectory());
+ m_ResampleFilter->SetDEMDirectory(otb::ConfigurationFile::GetInstance()->GetDEMDirectory());
}
}
- // Set the output map projection
- orthoFilter->SetOutputProjectionRef(outputProjectionRef);
-
- if (parseResult->IsOptionPresent("LocMapSpacing"))
- {
- double defScalarSpacing = parseResult->GetParameterFloat("LocMapSpacing");
- ImageType::SpacingType defSpacing;
- defSpacing[0] = defScalarSpacing;
- defSpacing[1] = -defScalarSpacing;
-
- orthoFilter->SetDeformationFieldSpacing(defSpacing);
- }
- else
- {
- // By default, generate a 10 times coarser deformation field
- ImageType::SpacingType defSpacing;
- defSpacing[0] = 10*spacing[0];
- defSpacing[1] = 10*spacing[1];
- orthoFilter->SetDeformationFieldSpacing(defSpacing);
- }
-
- // Set the interpolator type
- if (parseResult->IsOptionPresent("InterpolatorType"))
+ // If activated, generate RPC model
+ if(IsParameterEnabled("rpc.ncp"))
{
- int nbInterpParams = parseResult->GetNumberOfParameters("InterpolatorType") - 1;
- std::string typeInterpolator = parseResult->GetParameterString("InterpolatorType");
-
- if (typeInterpolator == "BCO")
- {
- BCOInterpolationType::Pointer interpolator = BCOInterpolationType::New();
- if( nbInterpParams == 1)
- interpolator->SetRadius(parseResult->GetParameterUInt("InterpolatorType", 1));
- orthoFilter->SetInterpolator(interpolator);
- }
- else if (typeInterpolator == "NEARESTNEIGHBOR" && nbInterpParams == 0)
- {
- NearestNeighborInterpolationType::Pointer interpolator = NearestNeighborInterpolationType::New();
- orthoFilter->SetInterpolator(interpolator);
- }
- else if (typeInterpolator == "LINEAR" && nbInterpParams == 0)
- {
- LinearInterpolationType::Pointer interpolator = LinearInterpolationType::New();
- orthoFilter->SetInterpolator(interpolator);
- }
- else
- {
- itkGenericExceptionMacro(<< "Interpolator type not recognized, choose one with (parameters) : BCO(0/1), NEARESTNEIGHBOR(0), LINEAR(0)");
- }
+ //std::cout <<"Estimating the rpc Model " <<std::endl;
+ m_ResampleFilter->EstimateInputRpcModelOn();
+ m_ResampleFilter->SetInputRpcGridSize( GetParameterInt("rpc.ncp"));
}
- else
+
+ // Set Output information
+ ResampleFilterType::SizeType size;
+ if (IsParameterEnabled("outputs.sizex") && IsParameterEnabled("outputs.sizey"))
{
- // Default is BCO
- BCOInterpolationType::Pointer interpolator = BCOInterpolationType::New();
- orthoFilter->SetInterpolator(interpolator);
+ size[0] = GetParameterInt("outputs.sizex");
+ size[1] = GetParameterInt("outputs.sizey");
+ m_ResampleFilter->SetOutputSize(size);
}
- //Update informations from orthorectification filter
- orthoFilter->GetOutput()->UpdateOutputInformation();
-
- //Instantiate the writer
- WriterType::Pointer writer = WriterType::New();
- writer->SetFileName(parseResult->GetOutputImage());
- writer->SetInput(orthoFilter->GetOutput());
-
- unsigned int ram = 256;
- if (parseResult->IsOptionPresent("AvailableMemory"))
+ ResampleFilterType::SpacingType spacing;
+ if (IsParameterEnabled("outputs.spacingx") && IsParameterEnabled("outputs.spacingy"))
{
- ram = parseResult->GetParameterUInt("AvailableMemory");
+ spacing[0] = GetParameterFloat("outputs.spacingx");
+ spacing[1] = GetParameterFloat("outputs.spacingy");
+ m_ResampleFilter->SetOutputSpacing(spacing);
}
- writer->SetAutomaticTiledStreaming(ram);
-
- otb::StandardWriterWatcher watcher(writer, orthoFilter,"Orthorectification");
- writer->Update();
- }
- catch ( itk::ExceptionObject & err )
- {
- std::cout << "Exception itk::ExceptionObject raised !" << std::endl;
- std::cout << err << std::endl;
- return EXIT_FAILURE;
- }
- catch ( std::bad_alloc & err )
- {
- std::cout << "Exception bad_alloc : "<<(char*)err.what()<< std::endl;
- return EXIT_FAILURE;
- }
- catch ( ... )
- {
- std::cout << "Unknown exception raised !" << std::endl;
- return EXIT_FAILURE;
- }
- return EXIT_SUCCESS;
-}
-
-
-int OrthoRectification::Execute(otb::ApplicationOptionsResult* parseResult)
-{
- if ( parseResult->IsOptionPresent("MapProjectionType"))
- {
- std::string typeMap = parseResult->GetParameterString("MapProjectionType", 0);
- int nbParams = parseResult->GetNumberOfParameters("MapProjectionType");
- nbParams--;
-
- if ((typeMap == "UTM")&&(nbParams==2))
+
+ ResampleFilterType::OriginType ul;
+ if (IsParameterEnabled("outputs.ulx") && IsParameterEnabled("outputs.uly"))
{
- int numZone = parseResult->GetParameterUInt("MapProjectionType", 1);
- char hemisphere = parseResult->GetParameterChar("MapProjectionType", 2);
-
- typedef otb::UtmInverseProjection UtmProjectionType;
- UtmProjectionType::Pointer utmProjection = UtmProjectionType::New();
-
- utmProjection->SetZone(numZone);
- utmProjection->SetHemisphere(hemisphere);
-
- return generic_main<UtmProjectionType>(parseResult, utmProjection );
+ ul[0] = GetParameterFloat("outputs.ulx");
+ ul[1] = GetParameterFloat("outputs.uly");
+ m_ResampleFilter->SetOutputOrigin(ul);
}
- else
- {
- std::vector<double> parameters;
-
- for (int i=1; i<nbParams+1; i++)
- {
- parameters.push_back(parseResult->GetParameterDouble("MapProjectionType", i));
- }
-
- if ((typeMap == "LAMBERT")&&(nbParams==4))
- {
- typedef otb::LambertConformalConicInverseProjection LambertProjectionType;
- LambertProjectionType::Pointer lambertProjection = LambertProjectionType::New();
- lambertProjection->SetParameters(parameters[0], parameters[1], parameters[2], parameters[3]);
-
- return generic_main<LambertProjectionType>(parseResult, lambertProjection);
- }
- else if ((typeMap == "LAMBERT2")&&(nbParams==0))
- {
- typedef otb::Lambert2EtenduInverseProjection Lambert2ProjectionType;
- Lambert2ProjectionType::Pointer lambert2Projection = Lambert2ProjectionType::New();
-
- return generic_main<Lambert2ProjectionType>(parseResult, lambert2Projection);
- }
- else if ((typeMap == "LAMBERT93")&&(nbParams==0))
- {
- typedef otb::Lambert93InverseProjection Lambert93ProjectionType;
- Lambert93ProjectionType::Pointer lambert93Projection = Lambert93ProjectionType::New();
+// std::cout <<"Resampler : Output Origin "<<m_ResampleFilter->GetOutputOrigin() << std::endl;
+// std::cout <<"Resampler : Output Spacing "<<m_ResampleFilter->GetOutputSpacing() << std::endl;
+// std::cout <<"Resampler : Output Size "<< m_ResampleFilter->GetOutputSize()<< std::endl;
- return generic_main<Lambert93ProjectionType>(parseResult, lambert93Projection);
- }
- else if ((typeMap == "SINUS")&&(nbParams==2))
- {
- typedef otb::SinusoidalInverseProjection SinusoidalProjectionType;
- SinusoidalProjectionType::Pointer sinusoidalProjection = SinusoidalProjectionType::New();
-
- sinusoidalProjection->SetParameters(parameters[0], parameters[1]);
-
- return generic_main<SinusoidalProjectionType>(parseResult, sinusoidalProjection);
- }
- else if ((typeMap == "ECKERT4")&&(nbParams==2))
- {
- typedef otb::Eckert4InverseProjection Eckert4ProjectionType;
- Eckert4ProjectionType::Pointer eckert4Projection = Eckert4ProjectionType::New();
-
- eckert4Projection->SetParameters(parameters[0], parameters[1]);
-
- return generic_main<Eckert4ProjectionType>(parseResult, eckert4Projection);
- }
- else if ((typeMap == "TRANSMERCATOR")&&(nbParams==3))
- {
- typedef otb::TransMercatorInverseProjection TransMercatorProjectionType;
- TransMercatorProjectionType::Pointer transMercatorProjection = TransMercatorProjectionType::New();
-
- transMercatorProjection->SetParameters(parameters[0], parameters[1], parameters[2]);
+ // Output Image
+ SetParameterOutputImage("out", m_ResampleFilter->GetOutput());
+ }
- return generic_main<TransMercatorProjectionType>(parseResult, transMercatorProjection);
- }
- else if ((typeMap == "MOLLWEID")&&(nbParams==2))
- {
- typedef otb::MollweidInverseProjection MollweidProjectionType;
- MollweidProjectionType::Pointer mollweidProjection = MollweidProjectionType::New();
+ ResampleFilterType::Pointer m_ResampleFilter;
+ std::string m_OutputProjectionRef;
+ };
- mollweidProjection->SetParameters(parameters[0], parameters[1]);
- return generic_main<MollweidProjectionType>(parseResult, mollweidProjection);
- }
- else if ((typeMap == "EPSG")&&(nbParams==1))
- {
- return epsg_main(parseResult);
- }
- else
- {
- itkGenericExceptionMacro(<< "TypeMap not recognized, choose one with (parameters) : UTM(2), LAMBERT(4), LAMBERT2(0), LAMBERT93(0), SINUS(2), ECKERT4(2), TRANSMERCATOR(3), MOLLWEID(2), EPSG(1)");
- }
- }
- }
- else
- {
- //TODO get utm zone and hemispher and build wkt
- typedef otb::VectorImage<unsigned short int, 2> ImageType;
- typedef otb::ImageFileReader<ImageType> ReaderType;
-
- typedef otb::UtmInverseProjection UtmProjectionType;
- typedef otb::GenericRSResampleImageFilter< ImageType, ImageType> ResamplerImageType;
- typedef ResamplerImageType::GenericRSTransformType TransformType;
-
- // Read input image information
- ReaderType::Pointer reader=ReaderType::New();
- reader->SetFileName(parseResult->GetInputImage().c_str());
- reader->GenerateOutputInformation();
-
- ImageType::Pointer input = reader->GetOutput();
- TransformType::Pointer transform=TransformType::New();
-
- transform->SetOutputDictionary(input->GetMetaDataDictionary());
- transform->SetOutputProjectionRef(input->GetProjectionRef());
- transform->SetOutputKeywordList(input->GetImageKeywordlist());
- //TODO is DEM informations mandatory here?
- //transform->SetDEMDirectory();
-
- // Needed variable
- std::string projectionRef;
- // The inverse transform is need here
- TransformType::Pointer invTransform = TransformType::New();
- transform->GetInverse(invTransform);
-
- // Build the UTM transform : Need the zone & the hemisphere
- // For this we us the geographic coordinate of the input UL corner
- typedef itk::Point<double, 2> GeoPointType;
-
- // get the utm zone and hemisphere using the input UL corner
- // geographic coordinates
- ImageType::PointType pSrc;
- ImageType::IndexType index;
- GeoPointType geoPoint;
- index[0] = input->GetLargestPossibleRegion().GetIndex()[0];
- index[1] = input->GetLargestPossibleRegion().GetIndex()[1];
- input->TransformIndexToPhysicalPoint(index, pSrc);
-
- // The first transform of the inverse transform : input -> WGS84
- geoPoint = invTransform->GetTransform()->GetFirstTransform()->TransformPoint(pSrc);
-
- // Guess the zone and the hemisphere
- int zone = otb::Utils::GetZoneFromGeoPoint(geoPoint[0], geoPoint[1]);
- bool hem = (geoPoint[1]>1e-10)?true:false;
-
- typedef otb::UtmInverseProjection UtmProjectionType;
- UtmProjectionType::Pointer utmProjection = UtmProjectionType::New();
-
- utmProjection->SetZone(zone);
- utmProjection->SetHemisphere(hem);
-
- otbMsgDevMacro(<< "Guess the UTM parameters, zone: " << zone << " hemisphere: " << hem);
-
- return generic_main<UtmProjectionType>(parseResult, utmProjection);
- }
-}
+ } // namespace Wrapper
} // namespace otb
+
+OTB_APPLICATION_EXPORT(otb::Wrapper::OrthoRectification)
--
GitLab