otbGDALRPCTransformer.cxx 6.18 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
/*
 * Copyright (C) 2005-2020 Centre National d'Etudes Spatiales (CNES)
 *
 * This file is part of Orfeo Toolbox
 *
 *     https://www.orfeo-toolbox.org/
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include "otbGDALRPCTransformer.h"
22
#include <assert.h>
23
#include "cpl_string.h"
24
#include "otbDEMHandler.h"
25
26
27
28

namespace otb
{
GDALRPCTransformer::GDALRPCTransformer(double LineOffset, double SampleOffset, double LatOffset, double LonOffset, double HeightOffset,
29
				       double LineScale, double SampleScale, double LatScale, double LonScale, double HeightScale,
30
31
32
				       const double (&LineNum)[20], const double (&LineDen)[20], const double (&SampleNum)[20], const double (&SampleDen)[20],
               bool useDEM)
  : m_UseDEM(useDEM)
33
34
35
36
37
38
39
40
41
42
43
44
45
46
{
  // Offsets
  this->m_GDALRPCInfo.dfLINE_OFF = LineOffset;
  this->m_GDALRPCInfo.dfSAMP_OFF = SampleOffset;
  this->m_GDALRPCInfo.dfLAT_OFF = LatOffset;
  this->m_GDALRPCInfo.dfLONG_OFF = LonOffset;
  this->m_GDALRPCInfo.dfHEIGHT_OFF = HeightOffset;
  // Scales
  this->m_GDALRPCInfo.dfLINE_SCALE = LineScale;
  this->m_GDALRPCInfo.dfSAMP_SCALE = SampleScale;
  this->m_GDALRPCInfo.dfLAT_SCALE = LatScale;
  this->m_GDALRPCInfo.dfLONG_SCALE = LonScale;
  this->m_GDALRPCInfo.dfHEIGHT_SCALE = HeightScale;
  // Coefficients
47
48
49
50
  std::copy_n(LineNum, 20, this->m_GDALRPCInfo.adfLINE_NUM_COEFF);
  std::copy_n(LineDen, 20, this->m_GDALRPCInfo.adfLINE_DEN_COEFF);
  std::copy_n(SampleNum, 20, this->m_GDALRPCInfo.adfSAMP_NUM_COEFF);
  std::copy_n(SampleDen, 20, this->m_GDALRPCInfo.adfSAMP_DEN_COEFF);
51

52
53
54
55
56
57
  // min/max longitude and latitude (default values)
  this->m_GDALRPCInfo.dfMIN_LONG = -180.0;
  this->m_GDALRPCInfo.dfMIN_LAT = -90.0;
  this->m_GDALRPCInfo.dfMAX_LONG = 180.0;
  this->m_GDALRPCInfo.dfMAX_LAT = 90.0;

58
  otb::DEMHandler::GetInstance().AttachObserver(this);
59
60
61

  this->SetOption("RPC_MAX_ITERATIONS", "40");
  this->SetOption("RPC_PIXEL_ERROR_THRESHOLD",  "0.000001");
62
63
}

64
65
66
67
68
69
70
71
GDALRPCTransformer::GDALRPCTransformer(const Projection::RPCParam & param, bool useDEM)
  : GDALRPCTransformer( param.LineOffset, param.SampleOffset, param.LatOffset, param.LonOffset, param.HeightOffset,
               param.LineScale, param.SampleScale, param.LatScale, param.LonScale, param.HeightScale,
               param.LineNum, param.LineDen, param.SampleNum, param.SampleDen,
               useDEM)
{
}

72
73
74
75
GDALRPCTransformer::~GDALRPCTransformer()
{
  if(m_TransformArg != nullptr)
    GDALDestroyTransformer(m_TransformArg);
76
  CSLDestroy(m_Options);
77
78

  otb::DEMHandler::GetInstance().DetachObserver(this);
79
80
}

81
void GDALRPCTransformer::SetOption(const std::string& Name, const std::string& Value)
82
{
83
  this->m_Options = CSLSetNameValue(m_Options, Name.c_str(), Value.c_str());
84
85
86
87
88
89
90
91
92
93
94
  this->m_Modified = true;
}

void GDALRPCTransformer::SetPixErrThreshold(double PixErrThreshold)
{
  this->m_PixErrThreshold = PixErrThreshold;
  this->m_Modified = true;
}

void GDALRPCTransformer::Update()
{
95
96
97
98
99
  // We need a lock here because Update() is not called until the first call
  // to Forward/InverseTransform(), which might be done in a thread.
  const std::lock_guard<std::mutex> lock(m_Mutex);

  auto & demHandler = otb::DEMHandler::GetInstance();
100
101

  if (m_UseDEM)
102
  {
103
    if (demHandler.GetDEMCount() > 0)
104
    {
105
106
107
108
109
110
111
112
113
      if (demHandler.GetGeoidFile().empty())
      {
        this->SetOption("RPC_DEM", demHandler.DEM_DATASET_PATH);
      }
      else
      {
        this->SetOption("RPC_DEM", demHandler.DEM_SHIFTED_DATASET_PATH);
      }
      this->SetOption("RPC_DEM_MISSING_VALUE", std::to_string(demHandler.GetDefaultHeightAboveEllipsoid()));
114
115
116
    }
    else
    {
117
118
      // RPC height is used as a constant height offset applied to all points in case no DEM is set.
      this->SetOption("RPC_HEIGHT", std::to_string(demHandler.GetDefaultHeightAboveEllipsoid()));
119
120
121
    }
  }

122
123
  if(m_TransformArg != nullptr)
    GDALDestroyTransformer(m_TransformArg);
124
125
126
127
128
129
  this->m_TransformArg = GDALCreateRPCTransformer(&this->m_GDALRPCInfo, false, this->m_PixErrThreshold, this->m_Options);
  this->m_Modified = false;
}

bool GDALRPCTransformer::ForwardTransform(double* x, double* y, double* z, int nPointCount)
{
130
131
132
  assert(x);
  assert(y);
  assert(z);
133
134
  if (this->m_Modified)
    this->Update();
135
  std::vector<int> success(nPointCount);
136
  const std::lock_guard<std::mutex> lock(m_Mutex);
137
138
  GDALRPCTransform(this->m_TransformArg, false, nPointCount, x, y, z, success.data());
  bool finalSuccess = std::all_of(success.begin(), success.end(), [](int i){return i;});
139
  return finalSuccess;
140
141
}

142
143
144
145
146
GDALRPCTransformer::PointType GDALRPCTransformer::ForwardTransform(GDALRPCTransformer::PointType p)
{
  if (m_Modified)
    this->Update();
  int success;
147
  const std::lock_guard<std::mutex> lock(m_Mutex);
148
149
150
151
  GDALRPCTransform(this->m_TransformArg, false, 1, &p[0], &p[1], &p[2], &success);
  if (!success)
    throw std::runtime_error("GDALRPCTransform was not able to process the ForwardTransform.");
  return p;
152
153
154
}


155
bool GDALRPCTransformer::InverseTransform(double* x, double* y, double* z, int nPointCount)
156
{
157
158
159
  assert(x);
  assert(y);
  assert(z);
160
161
  if (this->m_Modified)
    this->Update();
162
  std::vector<int> success(nPointCount);
163
  const std::lock_guard<std::mutex> lock(m_Mutex);
164
165
  GDALRPCTransform(this->m_TransformArg, true, nPointCount, x, y, z, success.data());
  bool finalSuccess = std::all_of(success.begin(), success.end(), [](int i){return i;});
166
167
168
169
170
171
172
173
  return finalSuccess;
}

GDALRPCTransformer::PointType GDALRPCTransformer::InverseTransform(GDALRPCTransformer::PointType p)
{
  if (m_Modified)
    this->Update();
  int success;
174
  const std::lock_guard<std::mutex> lock(m_Mutex);
175
176
177
178
  GDALRPCTransform(this->m_TransformArg, true, 1, &p[0], &p[1], &p[2], &success);
  if (!success)
    throw std::runtime_error("GDALRPCTransform was not able to process the InverseTransform.");
  return p;
179
180
}
}