otbOrthoRectification.cxx 33.6 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
/*=========================================================================

 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.

 =========================================================================*/
18 19 20
#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"
#include "otbWrapperNumericalParameter.h"
21

22
#include "otbImageToGenericRSOutputParameters.h"
23
#include "otbGenericRSResampleImageFilter.h"
24 25 26 27

#include "itkLinearInterpolateImageFunction.h"
#include "otbBCOInterpolateImageFunction.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
28

OTB Bot's avatar
STYLE  
OTB Bot committed
29
// MapProjection handler
30 31
#include "otbWrapperMapProjectionParametersHandler.h"

32 33 34
// Elevation handler
#include "otbWrapperElevationParametersHandler.h"

35 36
#include "otbGeographicalDistance.h"

37
#include "otbMacro.h"
38

39 40 41
namespace otb
{

42

43 44
enum
{
45
  Interpolator_BCO,
46
  Interpolator_NNeighbor,
47
  Interpolator_Linear
48
};
49

50 51 52 53
enum
{
  Mode_UserDefined,
  Mode_AutomaticSize,
54 55 56
  Mode_AutomaticSpacing,
  Mode_OutputROI,
  Mode_OrthoFit
57 58
};

59 60
const float DefaultGridSpacingMeter = 4.0;

61 62
namespace Wrapper
{
63

64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
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;
81

82
  /** Interpolators typedefs*/
Otmane Lahlou's avatar
STYLE  
Otmane Lahlou committed
83
  typedef itk::LinearInterpolateImageFunction<FloatVectorImageType,
84 85 86 87 88 89
                                              double>              LinearInterpolationType;
  typedef itk::NearestNeighborInterpolateImageFunction<FloatVectorImageType,
                                                       double>     NearestNeighborInterpolationType;
  typedef otb::BCOInterpolateImageFunction<FloatVectorImageType>   BCOInterpolationType;

private:
90
  void DoInit()
91 92
  {
    SetName("OrthoRectification");
Otmane Lahlou's avatar
STYLE  
Otmane Lahlou committed
93
    std::ostringstream oss;
94
    oss << "This application allows to ortho-rectify optical images from supported sensors." << std::endl;
Otmane Lahlou's avatar
STYLE  
Otmane Lahlou committed
95
    SetDescription(oss.str());
96
    // Documentation
97
    SetDocName("Ortho-rectification");
98 99 100 101 102
    oss.str("");
    oss<<"An inverse sensor model is built from the input image metadata to convert geographical to raw geometry coordinates. ";
    oss<<"This inverse sensor model is then combined with the chosen map projection to build a global coordinate mapping grid. Last, this grid is used to resample using the chosen interpolation algorithm. ";
    oss<<"A Digital Elevation Model can be specified to account for terrain deformations. "<<std::endl;
    oss<<"In case of SPOT5 images, the sensor model can be approximated by an RPC model in order to speed-up computation.";
103
    SetDocLongDescription(oss.str());
104
    SetDocLimitations("Supported sensors are Pleiades, SPOT5 (TIF format), Ikonos, Quickbird, Worldview2, GeoEye.");
105
    SetDocAuthors("OTB-Team");
106
    SetDocSeeAlso("Ortho-rectification chapter from the OTB Software Guide");
107

108
    AddDocTag(Tags::Geometry);
109

110
    // Set the parameters
111 112 113 114 115 116
    AddParameter(ParameterType_Group,"io","Input and output data");
    SetParameterDescription("io","This group of parameters allows to set the input and output images.");
    AddParameter(ParameterType_InputImage, "io.in", "Input Image");
    SetParameterDescription("io.in","The input image to ortho-rectify");
    AddParameter(ParameterType_OutputImage, "io.out", "Output Image");
    SetParameterDescription("io.out","The ortho-rectified output image");
117

118 119
    // Build the Output Map Projection
    MapProjectionParametersHandler::AddMapProjectionParameters(this, "map");
120

121
    // Add the output paramters in a group
122 123
    AddParameter(ParameterType_Group, "outputs", "Output Image Grid");
    SetParameterDescription("outputs","This group of parameters allows to define the grid on which the input image will be resampled.");
124 125

    // UserDefined values
126
    AddParameter(ParameterType_Choice, "outputs.mode", "Parameters estimation modes");
Otmane Lahlou's avatar
STYLE  
Otmane Lahlou committed
127
    AddChoice("outputs.mode.auto", "User Defined");
128 129 130 131 132
    SetParameterDescription("outputs.mode.auto","This mode allows you to fully modify default values.");
    AddChoice("outputs.mode.autosize", "Automatic Size from Spacing");
    SetParameterDescription("outputs.mode.autosize","This mode allows you to automatically compute the optimal image size from given spacing (pixel size) values");
    AddChoice("outputs.mode.autospacing", "Automatic Spacing from Size");
    SetParameterDescription("outputs.mode.autospacing","This mode allows you to automatically compute the optimal image spacing (pixel size) from the given size");
133 134 135 136
    AddChoice("outputs.mode.outputroi", "Automatic Size from Spacing and output corners");
    SetParameterDescription("outputs.mode.outputroi","This mode allows you to automatically compute the optimal image size from spacing (pixel size) and output corners");
    AddChoice("outputs.mode.orthofit", "Fit to ortho");
    SetParameterDescription("outputs.mode.orthofit", "Fit the size, origin and spacing to an existing ortho image (uses the value of outputs.ortho)");
137

138
    // Upper left point coordinates
Otmane Lahlou's avatar
STYLE  
Otmane Lahlou committed
139
    AddParameter(ParameterType_Float, "outputs.ulx", "Upper Left X");
140
    SetParameterDescription("outputs.ulx","Cartographic X coordinate of upper-left corner (meters for cartographic projections, degrees for geographic ones)");
141 142

    AddParameter(ParameterType_Float, "outputs.uly", "Upper Left Y");
143
    SetParameterDescription("outputs.uly","Cartographic Y coordinate of the upper-left corner (meters for cartographic projections, degrees for geographic ones)");
144 145 146

    // Size of the output image
    AddParameter(ParameterType_Int, "outputs.sizex", "Size X");
147
    SetParameterDescription("outputs.sizex","Size of projected image along X (in pixels)");
148 149

    AddParameter(ParameterType_Int, "outputs.sizey", "Size Y");
150
    SetParameterDescription("outputs.sizey","Size of projected image along Y (in pixels)");
151

152 153
    // Spacing of the output image
    AddParameter(ParameterType_Float, "outputs.spacingx", "Pixel Size X");
154
    SetParameterDescription("outputs.spacingx","Size of each pixel along X axis (meters for cartographic projections, degrees for geographic ones)");
155

156 157

    AddParameter(ParameterType_Float, "outputs.spacingy", "Pixel Size Y");
158
    SetParameterDescription("outputs.spacingy","Size of each pixel along Y axis (meters for cartographic projections, degrees for geographic ones)");
159

160 161 162 163 164 165
    // Lower right point coordinates
    AddParameter(ParameterType_Float, "outputs.lrx", "Lower right X");
    SetParameterDescription("outputs.lrx","Cartographic X coordinate of the lower-right corner (meters for cartographic projections, degrees for geographic ones)");

    AddParameter(ParameterType_Float, "outputs.lry", "Lower right Y");
    SetParameterDescription("outputs.lry","Cartographic Y coordinate of the lower-right corner (meters for cartographic projections, degrees for geographic ones)");
166

167
    // Existing ortho image that can be used to compute size, origin and spacing of the output
168
    AddParameter(ParameterType_InputImage, "outputs.ortho", "Model ortho-image");
169
    SetParameterDescription("outputs.ortho","A model ortho-image that can be used to compute size, origin and spacing of the output");
170

171 172 173 174 175 176 177
    // Setup parameters for initial UserDefined mode
    DisableParameter("outputs.lrx");
    DisableParameter("outputs.lry");
    DisableParameter("outputs.ortho");
    MandatoryOff("outputs.lrx");
    MandatoryOff("outputs.lry");
    MandatoryOff("outputs.ortho");
178

179
    AddParameter(ParameterType_Empty,"outputs.isotropic","Force isotropic spacing by default");
180
    std::ostringstream isotropOss;
181
    isotropOss << "Default spacing (pixel size) values are estimated from the sensor modeling of the image. It can therefore result in a non-isotropic spacing. ";
182 183
    isotropOss << "This option allows you to force default values to be isotropic (in this case, the minimum of spacing in both direction is applied. ";
    isotropOss << "Values overriden by user are not affected by this option.";
OTB Bot's avatar
STYLE  
OTB Bot committed
184
    SetParameterDescription("outputs.isotropic", isotropOss.str());
185
    EnableParameter("outputs.isotropic");
186

187
    AddParameter(ParameterType_Float, "outputs.default", "Default pixel value");
188
    SetParameterDescription("outputs.default","Default value to write when outside of input image.");
189
    SetDefaultParameterFloat("outputs.default",0.);
190
    MandatoryOff("outputs.default");
191

192 193
    // Elevation
    ElevationParametersHandler::AddElevationParameters(this, "elev");
194 195

    // Interpolators
196
    AddParameter(ParameterType_Choice,   "interpolator", "Interpolation");
197 198 199
    AddChoice("interpolator.bco",    "Bicubic interpolation");
    AddParameter(ParameterType_Radius, "interpolator.bco.radius", "Radius for bicubic interpolation");
    SetParameterDescription("interpolator.bco.radius","This parameter allows to control the size of the bicubic interpolation filter. If the target pixel size is higher than the input pixel size, increasing this parameter will reduce aliasing artefacts.");
200 201 202 203 204
    SetParameterDescription("interpolator","This group of parameters allows to define how the input image will be interpolated during resampling.");
    AddChoice("interpolator.nn",     "Nearest Neighbor interpolation");
    SetParameterDescription("interpolator.nn","Nearest neighbor interpolation leads to poor image quality, but it is very fast.");
    AddChoice("interpolator.linear", "Linear interpolation");
    SetParameterDescription("interpolator.linear","Linear interpolation leads to average image quality but is quite fast");
205
    SetDefaultParameterInt("interpolator.bco.radius", 2);
206 207
    AddParameter(ParameterType_Group,"opt","Speed optimization parameters");
    SetParameterDescription("opt","This group of parameters allows to optimize processing time.");
208

209 210
    // Estimate a RPC model (for spot image for instance)
    AddParameter(ParameterType_Int, "opt.rpc", "RPC modeling (points per axis)");
OTB Bot's avatar
STYLE  
OTB Bot committed
211
    SetDefaultParameterInt("opt.rpc", 10);
212 213 214
    SetParameterDescription("opt.rpc","Enabling RPC modeling allows to speed-up SPOT5 ortho-rectification. Value is the number of control points per axis for RPC estimation");
    DisableParameter("opt.rpc");
    MandatoryOff("opt.rpc");
215

216
    // RAM available
217
    AddRAMParameter("opt.ram");
218
    SetParameterDescription("opt.ram","This allows to set the maximum amount of RAM available for processing. As the writing task is time consuming, it is better to write large pieces of data, which can be achieved by increasing this parameter (pay attention to your system capabilities)");
219

220
    // Deformation Field Spacing
221
    AddParameter(ParameterType_Float, "opt.gridspacing", "Resampling grid spacing");
222
    SetDefaultParameterFloat("opt.gridspacing", DefaultGridSpacingMeter);
223 224
    SetParameterDescription("opt.gridspacing",
                            "Resampling is done according to a coordinate mapping deformation grid, "
225 226
                            "whose pixel size is set by this parameter, and "
                            "expressed in the coordinate system of the output image "
227 228 229
                            "The closer to the output spacing this parameter is, "
                            "the more precise will be the ortho-rectified image,"
                            "but increasing this parameter will reduce processing time.");
230
    MandatoryOff("opt.gridspacing");
231

232 233 234
    // Doc example parameter settings
    SetDocExampleParameterValue("io.in", "QB_TOULOUSE_MUL_Extract_500_500.tif");
    SetDocExampleParameterValue("io.out","QB_Toulouse_ortho.tif");
235
  }
236

237 238
  void DoUpdateParameters()
  {
239
    if (HasValue("io.in"))
240
      {
OTB Bot's avatar
STYLE  
OTB Bot committed
241
      // input image
242
      FloatVectorImageType::Pointer inImage = GetParameterImage("io.in");
243

244
      // Update the UTM zone params
245 246
      MapProjectionParametersHandler::InitializeUTMParameters(this, "io.in", "map");

247
      // Get the output projection Ref
248
      m_OutputProjectionRef = MapProjectionParametersHandler::GetProjectionRefFromChoice(this, "map");
249 250 251 252

      // Compute the output image spacing and size
      typedef otb::ImageToGenericRSOutputParameters<FloatVectorImageType> OutputParametersEstimatorType;
      OutputParametersEstimatorType::Pointer genericRSEstimator = OutputParametersEstimatorType::New();
253

254 255 256 257 258 259 260 261
      if(IsParameterEnabled("outputs.isotropic"))
        {
        genericRSEstimator->EstimateIsotropicSpacingOn();
        }
      else
        {
        genericRSEstimator->EstimateIsotropicSpacingOff();
        }
262

263 264 265
      genericRSEstimator->SetInput(inImage);
      genericRSEstimator->SetOutputProjectionRef(m_OutputProjectionRef);
      genericRSEstimator->Compute();
266

OTB Bot's avatar
STYLE  
OTB Bot committed
267
      // Fill the Gui with the computed parameters
268
      if (!HasUserValue("outputs.ulx"))
269
        {
270
        SetParameterFloat("outputs.ulx", genericRSEstimator->GetOutputOrigin()[0]);
271
        }
272

273 274 275 276 277 278 279 280 281 282 283 284 285 286 287
      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]);

288 289 290 291 292
      if (!HasUserValue("outputs.lrx"))
       SetParameterFloat("outputs.lrx", GetParameterFloat("outputs.ulx") + GetParameterFloat("outputs.spacingx") * static_cast<double>(GetParameterInt("outputs.sizex")));

      if (!HasUserValue("outputs.lry"))
       SetParameterFloat("outputs.lry", GetParameterFloat("outputs.uly") + GetParameterFloat("outputs.spacingy") * static_cast<double>(GetParameterInt("outputs.sizey")));
293

OTB Bot's avatar
STYLE  
OTB Bot committed
294
      // Handle the spacing and size field following the mode
295 296 297 298 299
      // choosed by the user
      switch (GetParameterInt("outputs.mode") )
        {
        case Mode_UserDefined:
        {
300 301 302
        // Automatic set to off except lower right coordinates
        AutomaticValueOff("outputs.ulx");
        AutomaticValueOff("outputs.uly");
303 304 305 306
        AutomaticValueOff("outputs.sizex");
        AutomaticValueOff("outputs.sizey");
        AutomaticValueOff("outputs.spacingx");
        AutomaticValueOff("outputs.spacingy");
307 308
        AutomaticValueOn("outputs.lrx");
        AutomaticValueOn("outputs.lry");
309

310 311 312
        // Enable all the parameters except lower right coordinates
        EnableParameter("outputs.ulx");
        EnableParameter("outputs.uly");
313 314 315 316
        EnableParameter("outputs.spacingx");
        EnableParameter("outputs.spacingy");
        EnableParameter("outputs.sizex");
        EnableParameter("outputs.sizey");
317 318 319
        DisableParameter("outputs.lrx");
        DisableParameter("outputs.lry");
        DisableParameter("outputs.ortho");
320

321 322 323
        // Make all the parameters in this mode mandatory except lower right coordinates
        MandatoryOn("outputs.ulx");
        MandatoryOn("outputs.uly");
324 325 326 327
        MandatoryOn("outputs.spacingx");
        MandatoryOn("outputs.spacingy");
        MandatoryOn("outputs.sizex");
        MandatoryOn("outputs.sizey");
328 329 330
        MandatoryOff("outputs.lrx");
        MandatoryOff("outputs.lry");
        MandatoryOff("outputs.ortho");
331

332 333 334
        // Update lower right
        SetParameterFloat("outputs.lrx", GetParameterFloat("outputs.ulx") + GetParameterFloat("outputs.spacingx") * static_cast<double>(GetParameterInt("outputs.sizex")));
        SetParameterFloat("outputs.lry", GetParameterFloat("outputs.uly") + GetParameterFloat("outputs.spacingy") * static_cast<double>(GetParameterInt("outputs.sizey")));
335 336 337 338 339
        }
        break;
        case Mode_AutomaticSize:
        {
        // Disable the size fields
340 341
        DisableParameter("outputs.ulx");
        DisableParameter("outputs.uly");
342 343 344 345
        DisableParameter("outputs.sizex");
        DisableParameter("outputs.sizey");
        EnableParameter("outputs.spacingx");
        EnableParameter("outputs.spacingy");
346 347 348
        DisableParameter("outputs.lrx");
        DisableParameter("outputs.lry");
        DisableParameter("outputs.ortho");
349

350
        // Update the automatic value mode of each filed
351 352
        AutomaticValueOn("outputs.ulx");
        AutomaticValueOn("outputs.uly");
353 354 355 356
        AutomaticValueOn("outputs.sizex");
        AutomaticValueOn("outputs.sizey");
        AutomaticValueOff("outputs.spacingx");
        AutomaticValueOff("outputs.spacingy");
357 358
        AutomaticValueOn("outputs.lrx");
        AutomaticValueOn("outputs.lry");
359

360
        // Adapat the status of the param to this mode
361 362
        MandatoryOff("outputs.ulx");
        MandatoryOff("outputs.uly");
363 364 365 366
        MandatoryOn("outputs.spacingx");
        MandatoryOn("outputs.spacingy");
        MandatoryOff("outputs.sizex");
        MandatoryOff("outputs.sizey");
367 368 369
        MandatoryOff("outputs.lrx");
        MandatoryOff("outputs.lry");
        MandatoryOff("outputs.ortho");
370

371 372 373
        ResampleFilterType::SpacingType spacing;
        spacing[0] = GetParameterFloat("outputs.spacingx");
        spacing[1] = GetParameterFloat("outputs.spacingy");
374

375 376
        genericRSEstimator->ForceSpacingTo(spacing);
        genericRSEstimator->Compute();
377

378 379 380
        // Set the  processed size relative to this forced spacing
        SetParameterInt("outputs.sizex", genericRSEstimator->GetOutputSize()[0]);
        SetParameterInt("outputs.sizey", genericRSEstimator->GetOutputSize()[1]);
381

382 383 384
        // Reset Origin to default
        SetParameterFloat("outputs.ulx", genericRSEstimator->GetOutputOrigin()[0]);
        SetParameterFloat("outputs.uly", genericRSEstimator->GetOutputOrigin()[1]);
385

386 387 388
        // Update lower right
        SetParameterFloat("outputs.lrx", GetParameterFloat("outputs.ulx") + GetParameterFloat("outputs.spacingx") * static_cast<double>(GetParameterInt("outputs.sizex")));
        SetParameterFloat("outputs.lry", GetParameterFloat("outputs.uly") + GetParameterFloat("outputs.spacingy") * static_cast<double>(GetParameterInt("outputs.sizey")));
389 390 391 392
        }
        break;
        case Mode_AutomaticSpacing:
        {
393
        // Disable the spacing fields and enable the size fields
394 395
        DisableParameter("outputs.ulx");
        DisableParameter("outputs.uly");
396 397 398 399
        DisableParameter("outputs.spacingx");
        DisableParameter("outputs.spacingy");
        EnableParameter("outputs.sizex");
        EnableParameter("outputs.sizey");
400 401 402
        DisableParameter("outputs.lrx");
        DisableParameter("outputs.lry");
        DisableParameter("outputs.ortho");
403

404
        // Update the automatic value mode of each filed
405 406
        AutomaticValueOn("outputs.ulx");
        AutomaticValueOn("outputs.uly");
407 408 409 410
        AutomaticValueOn("outputs.spacingx");
        AutomaticValueOn("outputs.spacingy");
        AutomaticValueOff("outputs.sizex");
        AutomaticValueOff("outputs.sizey");
411 412
        AutomaticValueOn("outputs.lrx");
        AutomaticValueOn("outputs.lry");
413

414
        // Adapat the status of the param to this mode
415 416
        MandatoryOff("outputs.ulx");
        MandatoryOff("outputs.uly");
417 418 419 420
        MandatoryOff("outputs.spacingx");
        MandatoryOff("outputs.spacingy");
        MandatoryOn("outputs.sizex");
        MandatoryOn("outputs.sizey");
421 422 423
        MandatoryOff("outputs.lrx");
        MandatoryOff("outputs.lry");
        MandatoryOff("outputs.ortho");
424

OTB Bot's avatar
STYLE  
OTB Bot committed
425
        ResampleFilterType::SizeType size;
426 427
        size[0] = GetParameterInt("outputs.sizex");
        size[1] = GetParameterInt("outputs.sizey");
428

429 430
        genericRSEstimator->ForceSizeTo(size);
        genericRSEstimator->Compute();
431

432 433 434
        // Set the  processed spacing relative to this forced size
        SetParameterFloat("outputs.spacingx", genericRSEstimator->GetOutputSpacing()[0]);
        SetParameterFloat("outputs.spacingy", genericRSEstimator->GetOutputSpacing()[1]);
435

436 437 438
        // Reset Origin to default
        SetParameterFloat("outputs.ulx", genericRSEstimator->GetOutputOrigin()[0]);
        SetParameterFloat("outputs.uly", genericRSEstimator->GetOutputOrigin()[1]);
439

440 441 442
        // Update lower right
        SetParameterFloat("outputs.lrx", GetParameterFloat("outputs.ulx") + GetParameterFloat("outputs.spacingx") * static_cast<double>(GetParameterInt("outputs.sizex")));
        SetParameterFloat("outputs.lry", GetParameterFloat("outputs.uly") + GetParameterFloat("outputs.spacingy") * static_cast<double>(GetParameterInt("outputs.sizey")));
443

444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482
        }
        break;
        case Mode_OutputROI:
        {
          // Set the size according to the output ROI
          EnableParameter("outputs.ulx");
          EnableParameter("outputs.uly");
          DisableParameter("outputs.sizex");
          DisableParameter("outputs.sizey");
          EnableParameter("outputs.spacingx");
          EnableParameter("outputs.spacingy");
          EnableParameter("outputs.lrx");
          EnableParameter("outputs.lry");
          DisableParameter("outputs.ortho");

          // Update the automatic value mode of each filed
          AutomaticValueOff("outputs.ulx");
          AutomaticValueOff("outputs.uly");
          AutomaticValueOn("outputs.sizex");
          AutomaticValueOn("outputs.sizey");
          AutomaticValueOff("outputs.spacingx");
          AutomaticValueOff("outputs.spacingy");
          AutomaticValueOff("outputs.lrx");
          AutomaticValueOff("outputs.lry");

          // Adapt the status of the param to this mode
          MandatoryOn("outputs.ulx");
          MandatoryOn("outputs.uly");
          MandatoryOn("outputs.spacingx");
          MandatoryOn("outputs.spacingy");
          MandatoryOff("outputs.sizex");
          MandatoryOff("outputs.sizey");
          MandatoryOn("outputs.lrx");
          MandatoryOn("outputs.lry");
          MandatoryOff("outputs.ortho");

          ResampleFilterType::SpacingType spacing;
          spacing[0] = GetParameterFloat("outputs.spacingx");
          spacing[1] = GetParameterFloat("outputs.spacingy");
483

484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502
          // Set the  processed size relative to this forced spacing
          if (vcl_abs(spacing[0]) > 0.0)
            SetParameterInt("outputs.sizex", static_cast<int>(vcl_ceil((GetParameterFloat("outputs.lrx")-GetParameterFloat("outputs.ulx"))/spacing[0])));
          if (vcl_abs(spacing[1]) > 0.0)
            SetParameterInt("outputs.sizey", static_cast<int>(vcl_ceil((GetParameterFloat("outputs.lry")-GetParameterFloat("outputs.uly"))/spacing[1])));
        }
        break;
        case Mode_OrthoFit:
        {
          // Make all the parameters in this mode mandatory
          MandatoryOff("outputs.ulx");
          MandatoryOff("outputs.uly");
          MandatoryOff("outputs.spacingx");
          MandatoryOff("outputs.spacingy");
          MandatoryOff("outputs.sizex");
          MandatoryOff("outputs.sizey");
          MandatoryOff("outputs.lrx");
          MandatoryOff("outputs.lry");
          MandatoryOn("outputs.ortho");
503

504 505 506 507 508 509 510 511 512 513
          // Disable the parameters
          DisableParameter("outputs.ulx");
          DisableParameter("outputs.uly");
          DisableParameter("outputs.spacingx");
          DisableParameter("outputs.spacingy");
          DisableParameter("outputs.sizex");
          DisableParameter("outputs.sizey");
          DisableParameter("outputs.lrx");
          DisableParameter("outputs.lry");
          EnableParameter("outputs.ortho");
514

515 516 517 518 519 520 521 522 523 524 525
          if (HasValue("outputs.ortho"))
          {
            // Automatic set to on
            AutomaticValueOn("outputs.ulx");
            AutomaticValueOn("outputs.uly");
            AutomaticValueOn("outputs.sizex");
            AutomaticValueOn("outputs.sizey");
            AutomaticValueOn("outputs.spacingx");
            AutomaticValueOn("outputs.spacingy");
            AutomaticValueOn("outputs.lrx");
            AutomaticValueOn("outputs.lry");
526

527 528
            // input image
            FloatVectorImageType::Pointer inOrtho = GetParameterImage("outputs.ortho");
529

530 531 532
            ResampleFilterType::OriginType orig = inOrtho->GetOrigin();
            ResampleFilterType::SpacingType spacing = inOrtho->GetSpacing();
            ResampleFilterType::SizeType size = inOrtho->GetLargestPossibleRegion().GetSize();
533

534 535 536 537 538 539 540 541 542 543
            SetParameterInt("outputs.sizex",size[0]);
            SetParameterInt("outputs.sizey",size[1]);
            SetParameterFloat("outputs.spacingx",spacing[0]);
            SetParameterFloat("outputs.spacingy",spacing[1]);
            SetParameterFloat("outputs.ulx", orig[0]);
            SetParameterFloat("outputs.uly", orig[1]);
            // Update lower right
            SetParameterFloat("outputs.lrx", GetParameterFloat("outputs.ulx") + GetParameterFloat("outputs.spacingx") * static_cast<double>(GetParameterInt("outputs.sizex")));
            SetParameterFloat("outputs.lry", GetParameterFloat("outputs.uly") + GetParameterFloat("outputs.spacingy") * static_cast<double>(GetParameterInt("outputs.sizey")));
          }
544 545
        }
        break;
546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589
        } // switch (GetParameterInt("outputs.mode") )

      if (!HasUserValue("opt.gridspacing"))
        {
        // Update opt.gridspacing
        // In case output coordinate system is WG84,
        if (m_OutputProjectionRef == otb::GeoInformationConversion::ToWKT(4326))
          {
          // How much is 4 meters in degrees ?
          typedef itk::Point<float,2> FloatPointType;
          FloatPointType point1, point2;

          typedef otb::GeographicalDistance<FloatPointType> GeographicalDistanceType;
          GeographicalDistanceType::Pointer geoDistanceCalculator = GeographicalDistanceType::New();

          // center
          point1[0] = GetParameterFloat("outputs.ulx") + GetParameterFloat("outputs.spacingx") * GetParameterInt("outputs.sizex") / 2;
          point1[1] = GetParameterFloat("outputs.uly") + GetParameterFloat("outputs.spacingy") * GetParameterInt("outputs.sizey") / 2;

          // center + [1,0]
          point2[0] = point1[0] + GetParameterFloat("outputs.spacingx");
          point2[1] = point1[1];
          double xgroundspacing = geoDistanceCalculator->Evaluate(point1, point2);
          otbAppLogINFO( "Output X ground spacing in meter = " << xgroundspacing );

          // center + [0,1]
          point2[0] = point1[0];
          point2[1] = point1[1] + GetParameterFloat("outputs.spacingy");
          double ygroundspacing = geoDistanceCalculator->Evaluate(point1, point2);
          otbAppLogINFO( "Output Y ground spacing in meter = " << ygroundspacing );

          double xgridspacing = DefaultGridSpacingMeter * GetParameterFloat("outputs.spacingx") / xgroundspacing;
          double ygridspacing = DefaultGridSpacingMeter * GetParameterFloat("outputs.spacingy") / ygroundspacing;

          otbAppLogINFO( << DefaultGridSpacingMeter << " meters in X direction correspond roughly to "
                         << xgridspacing << " degrees" );
          otbAppLogINFO( << DefaultGridSpacingMeter << " meters in Y direction correspond roughly to "
                         << ygridspacing << " degrees" );

          // Use the smallest spacing (more precise grid)
          double optimalSpacing = std::min( vcl_abs(xgridspacing), vcl_abs(ygridspacing) );
          otbAppLogINFO( "Setting grid spacing to " << optimalSpacing );
          SetParameterFloat("opt.gridspacing", optimalSpacing);
          }
590 591 592 593
        else // if (m_OutputProjectionRef == otb::GeoInformationConversion::ToWKT(4326))
          {
          SetParameterFloat("opt.gridspacing", DefaultGridSpacingMeter);
          } // if (m_OutputProjectionRef == otb::GeoInformationConversion::ToWKT(4326))
594 595
        } // if (!HasUserValue("opt.gridspacing"))
      } // if (HasValue("io.in"))
596
  }
597

598
  void DoExecute()
599
    {
600
    // Get the input image
601
    FloatVectorImageType* inImage = GetParameterImage("io.in");
602 603 604 605 606

    // Resampler Instanciation
    m_ResampleFilter = ResampleFilterType::New();
    m_ResampleFilter->SetInput(inImage);

OTB Bot's avatar
STYLE  
OTB Bot committed
607
    // Set the output projection Ref
608 609 610
    m_ResampleFilter->SetInputProjectionRef(inImage->GetProjectionRef());
    m_ResampleFilter->SetInputKeywordList(inImage->GetImageKeywordlist());
    m_ResampleFilter->SetOutputProjectionRef(m_OutputProjectionRef);
611

612 613 614 615 616
    // Check size
    if (GetParameterInt("outputs.sizex") <= 0 || GetParameterInt("outputs.sizey") <= 0)
      {
      otbAppLogCRITICAL("Wrong value : negative size : ("<<GetParameterInt("outputs.sizex")<<" , "<<GetParameterInt("outputs.sizey")<<")");
      }
617

618 619 620 621 622 623
    //Check spacing sign
    if (GetParameterFloat("outputs.spacingy") > 0.)
      {
      otbAppLogWARNING(<<"Wrong value for outputs.spacingy: Pixel size along Y axis should be negative, (outputs.spacingy=" <<GetParameterFloat("outputs.spacingy") << ")" )
      }

624 625
    // Get Interpolator
    switch ( GetParameterInt("interpolator") )
626
      {
627
      case Interpolator_Linear:
628
      {
Otmane Lahlou's avatar
STYLE  
Otmane Lahlou committed
629 630
      typedef itk::LinearInterpolateImageFunction<FloatVectorImageType,
        double>          LinearInterpolationType;
631 632
      LinearInterpolationType::Pointer interpolator = LinearInterpolationType::New();
      m_ResampleFilter->SetInterpolator(interpolator);
633
      }
634 635
      break;
      case Interpolator_NNeighbor:
636
      {
Otmane Lahlou's avatar
STYLE  
Otmane Lahlou committed
637 638
      typedef itk::NearestNeighborInterpolateImageFunction<FloatVectorImageType,
        double> NearestNeighborInterpolationType;
639 640
      NearestNeighborInterpolationType::Pointer interpolator = NearestNeighborInterpolationType::New();
      m_ResampleFilter->SetInterpolator(interpolator);
641
      }
642 643
      break;
      case Interpolator_BCO:
644
      {
Otmane Lahlou's avatar
STYLE  
Otmane Lahlou committed
645
      typedef otb::BCOInterpolateImageFunction<FloatVectorImageType>     BCOInterpolationType;
646 647 648
      BCOInterpolationType::Pointer interpolator = BCOInterpolationType::New();
      interpolator->SetRadius(GetParameterInt("interpolator.bco.radius"));
      m_ResampleFilter->SetInterpolator(interpolator);
649
      }
650
      break;
651 652
      }

653 654
    // Setup the DEM Handler
    otb::Wrapper::ElevationParametersHandler::SetupDEMHandlerFromElevationParameters(this,"elev");
655

656
    // If activated, generate RPC model
657
    if(IsParameterEnabled("opt.rpc"))
658
      {
659
      m_ResampleFilter->EstimateInputRpcModelOn();
660 661
      m_ResampleFilter->SetInputRpcGridSize( GetParameterInt("opt.rpc") );
      otbAppLogINFO("Generating RPC modeling with " << GetParameterInt("opt.rpc") << " points per axis");
662
      }
663

664 665
    // Set Output information
    ResampleFilterType::SizeType size;
666 667 668
    size[0] = GetParameterInt("outputs.sizex");
    size[1] = GetParameterInt("outputs.sizey");
    m_ResampleFilter->SetOutputSize(size);
669
    otbAppLogINFO("Generating output with size = " << size);
670

671
    ResampleFilterType::SpacingType spacing;
672 673 674
    spacing[0] = GetParameterFloat("outputs.spacingx");
    spacing[1] = GetParameterFloat("outputs.spacingy");
    m_ResampleFilter->SetOutputSpacing(spacing);
675
    otbAppLogINFO("Generating output with pixel spacing = " << spacing);
676

677
    ResampleFilterType::OriginType ul;
678 679 680
    ul[0] = GetParameterFloat("outputs.ulx");
    ul[1] = GetParameterFloat("outputs.uly");
    m_ResampleFilter->SetOutputOrigin(ul);
681
    otbAppLogINFO("Generating output with origin = " << ul);
682

683 684 685 686 687
    // Build the default pixel
    FloatVectorImageType::PixelType defaultValue;
    defaultValue.SetSize(inImage->GetNumberOfComponentsPerPixel());
    defaultValue.Fill(GetParameterFloat("outputs.default"));
    m_ResampleFilter->SetEdgePaddingValue(defaultValue);
688
    otbAppLogINFO("Area outside input image bounds will have a pixel value of " << defaultValue);
689 690


691 692
    // Deformation Field spacing
    ResampleFilterType::SpacingType gridSpacing;
693
    if (IsParameterEnabled("opt.gridspacing"))
694
      {
695 696
      gridSpacing[0] = GetParameterFloat("opt.gridspacing");
      gridSpacing[1] = -GetParameterFloat("opt.gridspacing");
697 698 699 700 701

      otbAppLogINFO("Using a deformation grid with a physical spacing of " << GetParameterFloat("opt.gridspacing"));

      // Predict size of deformation grid
      ResampleFilterType::SizeType deformationGridSize;
702 703 704 705
      deformationGridSize[0] = static_cast<ResampleFilterType::SizeType::SizeValueType>(vcl_abs(
          GetParameterInt("outputs.sizex") * GetParameterFloat("outputs.spacingx") / GetParameterFloat("opt.gridspacing") ));
      deformationGridSize[1] = static_cast<ResampleFilterType::SizeType::SizeValueType>(vcl_abs(
          GetParameterInt("outputs.sizey") * GetParameterFloat("outputs.spacingy") / GetParameterFloat("opt.gridspacing") ));
706
      otbAppLogINFO("Using a deformation grid of size " << deformationGridSize);
707 708 709

      if (deformationGridSize[0] * deformationGridSize[1] == 0)
        {
710
        otbAppLogFATAL("Deformation grid degenerated (size of 0). "
711 712
            "You shall set opt.gridspacing appropriately. "
            "opt.gridspacing units are the same as outputs.spacing units");
713 714 715 716 717
        }

      if (vcl_abs(GetParameterFloat("opt.gridspacing")) < vcl_abs(GetParameterFloat("outputs.spacingx"))
           || vcl_abs(GetParameterFloat("opt.gridspacing")) < vcl_abs(GetParameterFloat("outputs.spacingy")) )
        {
718 719 720 721 722
        otbAppLogWARNING("Spacing of deformation grid should be at least equal to "
            "spacing of output image. Otherwise, computation time will be slow, "
            "and precision of output will not be better. "
            "You shall set opt.gridspacing appropriately. "
            "opt.gridspacing units are the same as outputs.spacing units");
723 724
        }

725 726 727
      m_ResampleFilter->SetDeformationFieldSpacing(gridSpacing);
      }

OTB Bot's avatar
STYLE  
OTB Bot committed
728
    // Output Image
729
    SetParameterOutputImage("io.out", m_ResampleFilter->GetOutput());
730
    }
731

732 733
  ResampleFilterType::Pointer     m_ResampleFilter;
  std::string                     m_OutputProjectionRef;
734
  };
735

736
  } // namespace Wrapper
737
} // namespace otb
738 739

OTB_APPLICATION_EXPORT(otb::Wrapper::OrthoRectification)