OTB  8.1.0
Orfeo Toolbox
otbGCPsToRPCSensorModelImageFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2022 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 #ifndef otbGCPsToRPCSensorModelImageFilter_hxx
22 #define otbGCPsToRPCSensorModelImageFilter_hxx
23 
25 
26 #include "otbRPCSolver.h"
27 #include "otbGenericRSTransform.h"
28 
29 #include "itkMetaDataObject.h"
30 #include "otbMetaDataKey.h"
31 
32 namespace otb
33 {
34 
35 template <class TImage>
37  : m_UseImageGCPs(false),
38  m_RMSGroundError(0.),
39  m_ErrorsContainer(),
40  m_MeanError(0.),
41  m_UseDEM(false),
42  m_MeanElevation(0.),
43  m_GCPsContainer(),
44  m_ModelUpToDate(false)
45 {
46  // This filter does not modify the image buffer, but only its
47  // metadata.Therefore, it can be run inplace to reduce memory print.
48  // CastImageFilter has InPlaceOff by default (see UnaryFunctorImgeFilter constructor)
49  this->InPlaceOn();
50 
51  // Clear the GCPs container
52  this->ClearGCPs();
53 }
54 
55 template <class TImage>
57 {
58  // Clear the GCPs container
59  this->ClearGCPs();
60 }
61 
62 
63 template <class TImage>
65 {
66  // Call superclass implementation
67  this->Superclass::Modified();
68 
69  // Deactivate up-to-date flag
70  m_ModelUpToDate = false;
71 }
72 
73 
74 template <class TImage>
76 {
77  // return the GCPs container
78  return m_GCPsContainer;
79 }
80 
81 template <class TImage>
83 {
84  // return the GCPs container
85  return m_ErrorsContainer;
86 }
87 
88 template <class TImage>
90 {
91  // Set the GCPs container
92  m_GCPsContainer = container;
93  // Call modified method
94  this->Modified();
95 }
96 
97 template <class TImage>
98 void GCPsToRPCSensorModelImageFilter<TImage>::AddGCP(const Point2DType& sensorPoint, const Point3DType& groundPoint)
99 {
100  // Create a new GCP
101  GCPType newGCP(sensorPoint, groundPoint);
102  // Add it to the container
103  m_GCPsContainer.push_back(newGCP);
104  // Call the modified method */
105  this->Modified();
106 }
107 
108 template <class TImage>
109 void GCPsToRPCSensorModelImageFilter<TImage>::AddGCP(const Point2DType& sensorPoint, const Point2DType& groundPoint)
110 {
111  // Create the ground point with elevation
112  Point3DType groundPointWithElevation;
113  // Fill it with ground point data
114  groundPointWithElevation[0] = groundPoint[0];
115  groundPointWithElevation[1] = groundPoint[1];
116 
117  // Check whether we are using a DEM or not
118  if (m_UseDEM)
119  {
120  // If so, use it to get the elevation
121  double height = DEMHandler::GetInstance().GetHeightAboveEllipsoid(groundPoint);
122  // To avoid nan value
123  if (height != height)
124  height = 0;
125  groundPointWithElevation[2] = height;
126  }
127  else
128  {
129  // Use the MeanElevation value
130  groundPointWithElevation[2] = m_MeanElevation;
131  }
132 
133  // Call the 3D version of the method
134  this->AddGCP(sensorPoint, groundPointWithElevation);
135 }
136 
137 template <class TImage>
139 {
140  if (id < m_GCPsContainer.size())
141  {
142  m_GCPsContainer.erase(m_GCPsContainer.begin() + id);
143  }
144 
145  this->Modified();
146 }
147 
148 template <class TImage>
150 {
151  m_GCPsContainer.clear();
152  this->Modified();
153 }
154 
155 template <class TImage>
157 {
158  // Set the internal value
159  m_UseImageGCPs = use;
160 
161  this->LoadImageGCPs();
162 }
163 
164 template <class TImage>
166 {
167  // First, retrieve the image pointer
168  typename TImage::Pointer imagePtr = this->GetOutput();
169 
170  // Iterate on the image GCPs
171  for (unsigned int i = 0; i < imagePtr->GetGCPCount(); ++i)
172  {
173  // Store row/col in an index
174  typename TImage::IndexType index;
175  index[0] = static_cast<long int>(imagePtr->GetGCPCol(i));
176  index[1] = static_cast<long int>(imagePtr->GetGCPRow(i));
177 
178  // Transform to physical point
179  Point2DType sensorPoint;
180  imagePtr->TransformIndexToPhysicalPoint(index, sensorPoint);
181 
182  // Sensor and Ground points
183  Point3DType groundPoint;
184 
185  groundPoint[0] = imagePtr->GetGCPX(i); // Lon
186  groundPoint[1] = imagePtr->GetGCPY(i); // Lat
187  groundPoint[2] = imagePtr->GetGCPZ(i);
188 
189  // Search image GCPs and remove them from container
190  bool found = false;
191 
192  unsigned int currentGCP = 0;
193 
194  while ((currentGCP < m_GCPsContainer.size()) && !found)
195  {
196  if ((m_GCPsContainer[currentGCP].first == sensorPoint) && (m_GCPsContainer[currentGCP].second == groundPoint))
197  {
198  // If we don't use image GCPs, erase the GCP
199  if (!m_UseImageGCPs)
200  m_GCPsContainer.erase(m_GCPsContainer.begin() + currentGCP);
201 
202  found = true;
203  }
204 
205  currentGCP++;
206  }
207 
208  // If we use image GCPs, add it to the container
209  if (m_UseImageGCPs && !found)
210  this->AddGCP(sensorPoint, groundPoint);
211  }
212 
213  this->Modified();
214 }
215 
216 template <class TImage>
218 {
219  using RSTransformType = GenericRSTransform<double, 2, 2>;
220 
221  RSTransformType::Pointer rsTransform = RSTransformType::New();
222  rsTransform->SetInputImageMetadata(&m_ImageMetadata);
223 
224  rsTransform->InstantiateTransform();
225 
226  double sum = 0.;
227  m_MeanError = 0.;
228 
229  // Clear Error container
230  m_ErrorsContainer.clear();
231 
232  for (unsigned int i = 0; i < m_GCPsContainer.size(); ++i)
233  {
234  // GCP value
235  const auto & sensorPoint = m_GCPsContainer[i].first;
236  const auto & groundPoint = m_GCPsContainer[i].second;
237 
238  // Compute Transform
239  Point3DType groundPointTemp, sensorPointTemp;
240  sensorPointTemp[0] = sensorPoint[0];
241  sensorPointTemp[1] = sensorPoint[1];
242 
243  auto outPoint = rsTransform->TransformPoint(sensorPoint);
244 
245  Point2DType groundPoint2D;
246 
247  groundPoint2D[0] = groundPoint[0];
248  groundPoint2D[1] = groundPoint[1];
249 
250  double error = outPoint.EuclideanDistanceTo(outPoint);
251 
252  // Add error to the container
253  m_ErrorsContainer.push_back(error);
254 
255  // Compute mean error
256  sum += error;
257  }
258 
259  m_MeanError = sum / static_cast<double>(m_ErrorsContainer.size());
260 }
261 
262 template <class TImage>
264 {
265  // First, call the superclass implementation
266  Superclass::GenerateOutputInformation();
267 
268  // First, retrieve the image pointer
269  auto imagePtr = this->GetOutput();
270 
271  if (!m_ModelUpToDate)
272  {
273  double rmsError;
274 
275  Projection::RPCParam rpcParams;
276  otb::RPCSolver::Solve(m_GCPsContainer, rmsError, rpcParams);
277 
278  // Retrieve the residual ground error
279  m_RMSGroundError = rmsError;
280 
281  m_ImageMetadata = imagePtr->GetImageMetadata();
282  m_ImageMetadata.Add(MDGeom::RPC, rpcParams);
283 
284  // Compute errors
285  this->ComputeErrors();
286 
287  m_ModelUpToDate = true;
288  }
289 
290  imagePtr->SetImageMetadata(m_ImageMetadata);
291 }
292 
293 template <class TImage>
294 void GCPsToRPCSensorModelImageFilter<TImage>::PrintSelf(std::ostream& os, itk::Indent indent) const
295 {
296  Superclass::PrintSelf(os, indent);
297  os << indent << "UseImageGCPs: " << (m_UseImageGCPs ? "yes" : "no") << std::endl;
298  os << indent << "UseDEM: " << (m_UseDEM ? "yes" : "no") << std::endl;
299  if (!m_UseDEM)
300  {
301  os << indent << "MeanElevation: " << m_MeanElevation << std::endl;
302  }
303  os << indent << "RMS ground error: " << m_RMSGroundError << std::endl;
304 }
305 
306 } // end of namespace otb
307 
308 #endif
itk::Index< Monteverdi_DIMENSION > IndexType
Definition: mvdTypes.h:133
double GetHeightAboveEllipsoid(double lon, double lat) const
void PrintSelf(std::ostream &os, itk::Indent indent) const override
void SetGCPsContainer(const GCPsContainerType &container)
void Add(const MDNum &, const MetaDataKey::VariableLengthVectorType)
std::pair< Point2DType, Point3DType > GCPType
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
void Solve(const GCPsContainerType &gcpContainer, double &rmsError, Projection::RPCParam &outputParams)
void AddGCP(const Point2DType &sensorPoint, const Point3DType &groundPoint)
This is the class to handle generic remote sensing transform.
Coefficients for RPC model (quite similar to GDALRPCInfo)
static DEMHandler & GetInstance()