Orfeo Toolbox  4.0
otbGenericRSTransform.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12 
13  This software is distributed WITHOUT ANY WARRANTY; without even
14  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15  PURPOSE. See the above copyright notices for more information.
16 
17 =========================================================================*/
18 #ifndef __otbGenericRSTransform_txx
19 #define __otbGenericRSTransform_txx
20 
21 #include "otbGenericRSTransform.h"
22 #include "otbMacro.h"
23 #include "otbMetaDataKey.h"
24 #include "itkMetaDataObject.h"
25 
27 
28 #include "ogr_spatialref.h"
29 
30 namespace otb
31 {
32 
33 template<class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions>
36 {
37  m_InputProjectionRef.clear();
38  m_OutputProjectionRef.clear();
39  m_InputKeywordList.Clear();
40  m_OutputKeywordList.Clear();
41  m_InputSpacing.Fill(1);
42  m_InputOrigin.Fill(0);
43  m_OutputSpacing.Fill(1);
44  m_OutputOrigin.Fill(0);
45 
46  m_Transform = NULL;
47  m_InputTransform = NULL;
48  m_OutputTransform = NULL;
49  m_TransformUpToDate = false;
50  m_TransformAccuracy = Projection::UNKNOWN;
51 }
52 
53 template<class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions>
57 {
58  itkDebugMacro("returning MapProjection address " << this->m_Transform);
59  if ((!m_TransformUpToDate) || (m_Transform.IsNull()))
60  {
61  itkExceptionMacro(<< "m_Transform not up-to-date, call InstanciateTransform() first");
62  }
63 
64  return this->m_Transform;
65 }
66 
70 template<class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions>
71 void
74 {
75  m_Transform = TransformType::New();
76 
77  if (m_InputKeywordList.GetSize() == 0)
78  {
79  itk::ExposeMetaData<ImageKeywordlist>(m_InputDictionary, MetaDataKey::OSSIMKeywordlistKey, m_InputKeywordList);
80  }
81  if (m_InputProjectionRef.empty())
82  {
83  itk::ExposeMetaData<std::string>(m_InputDictionary, MetaDataKey::ProjectionRefKey, m_InputProjectionRef);
84  }
85 
86  otbMsgDevMacro(<< "Information to instanciate transform: ");
87  otbMsgDevMacro(<< " * Input Origin: " << m_InputOrigin);
88  otbMsgDevMacro(<< " * Input Spacing: " << m_InputSpacing);
89  otbMsgDevMacro(<< " * Input keyword list: "
90  << ((m_InputKeywordList.GetSize() == 0) ? "Empty" : "Full"));
91  otbMsgDevMacro(<< " * Input projection: " << m_InputProjectionRef);
92  otbMsgDevMacro(<< " * Output keyword list: "
93  << ((m_OutputKeywordList.GetSize() == 0) ? "Empty" : "Full"));
94  otbMsgDevMacro(<< " * Output projection: " << m_OutputProjectionRef);
95  otbMsgDevMacro(<< " * Output Origin: " << m_OutputOrigin);
96  otbMsgDevMacro(<< " * Output Spacing: " << m_OutputSpacing);
97 
98  //Make sure that the state is clean:
99  m_InputTransform = NULL;
100  m_OutputTransform = NULL;
101 
102  bool firstTransformGiveGeo = true;
103  bool inputTransformIsSensor = false;
104  bool inputTransformIsMap = false;
105  bool outputTransformIsSensor = false;
106  bool outputTransformIsMap = false;
107 
108  // Prepare the projection ref (eventually convert the EPSG code into full WKT)
109  //
110  // Note that we do that at the GenericRSTransform level and not in the member
111  // class for several reasons:
112  // - at the GenericMapProjection and MapProjectionAdapter the method are
113  // called SetWkt and thus should not take a SRID.
114  // - we do not want to mix the GeoInformationConversion (which uses gdal) in
115  // the MapProjectionAdapter to keep ossim and gdal dependencies as separated
116  // as possible.
117  m_InputProjectionRef = GeoInformationConversion::ToWKT(m_InputProjectionRef);
118  m_OutputProjectionRef = GeoInformationConversion::ToWKT(m_OutputProjectionRef);
119 
120  //*****************************
121  //Set the input transformation
122  //*****************************
123 
124  // First, try to make a geo transform
125  if (!m_InputProjectionRef.empty()) //map projection
126  {
128  InverseMapProjectionType;
129  typename InverseMapProjectionType::Pointer mapTransform = InverseMapProjectionType::New();
130 
131  mapTransform->SetWkt(m_InputProjectionRef);
132  if (mapTransform->IsProjectionDefined())
133  {
134  m_InputTransform = mapTransform.GetPointer();
135  inputTransformIsMap = true;
136  otbMsgDevMacro(<< "Input projection set to map transform: " << m_InputTransform);
137  }
138  }
139 
140  // If not, try to make a sensor model
141  if ((m_InputTransform.IsNull()) && (m_InputKeywordList.GetSize() > 0))
142  {
144  typename ForwardSensorModelType::Pointer sensorModel = ForwardSensorModelType::New();
145 
146  sensorModel->SetImageGeometry(m_InputKeywordList);
147 
148  if (sensorModel->IsValidSensorModel())
149  {
150  m_InputTransform = sensorModel.GetPointer();
151  inputTransformIsSensor = true;
152  otbMsgDevMacro(<< "Input projection set to sensor model.");
153  }
154  }
155 
156  if (m_InputTransform.IsNull()) //default if we didn't manage to instantiate it before
157  {
159 // firstTransformGiveGeo = false;
160 
161  OGRSpatialReferenceH hSRS = NULL;
162  hSRS = OSRNewSpatialReference(NULL);
163  const char * wktString = m_InputProjectionRef.c_str();
164  if (OSRImportFromWkt(hSRS, (char **) &wktString) != OGRERR_NONE)
165  {
166  firstTransformGiveGeo = false;
167  otbMsgDevMacro(<< "- Considering that the first transform does not give geo (WKT)")
168  }
169 
170  else if ( OSRIsGeographic(hSRS) )
171  {
172  firstTransformGiveGeo = true;
173  otbMsgDevMacro(<< "- Considering that the first transform gives geo")
174  }
175  else
176  {
177  firstTransformGiveGeo = false;
178  otbMsgDevMacro(<< "- Considering that the first transform does not give geo (fallback)")
179  }
180  OSRRelease(hSRS);
181  otbMsgDevMacro(<< "Input projection set to identity")
182  }
183 
184  //*****************************
185  //Set the output transformation
186  //*****************************
187  if (!m_OutputProjectionRef.empty()) //map projection
188  {
190  OutputSpaceDimension> ForwardMapProjectionType;
191  typename ForwardMapProjectionType::Pointer mapTransform = ForwardMapProjectionType::New();
192  mapTransform->SetWkt(m_OutputProjectionRef);
193  if (mapTransform->IsProjectionDefined())
194  {
195  m_OutputTransform = mapTransform.GetPointer();
196  outputTransformIsMap = true;
197  otbMsgDevMacro(<< "Output projection set to map transform: " << m_OutputTransform);
198  }
199  }
200 
201  // If not, try to make a sensor model
202  if ((m_OutputTransform.IsNull()) && (m_OutputKeywordList.GetSize() > 0))
203  {
205  typename InverseSensorModelType::Pointer sensorModel = InverseSensorModelType::New();
206 
207  sensorModel->SetImageGeometry(m_OutputKeywordList);
208 
209  if (sensorModel->IsValidSensorModel())
210  {
211  m_OutputTransform = sensorModel.GetPointer();
212  outputTransformIsSensor = true;
213  otbMsgDevMacro(<< "Output projection set to sensor model");
214  }
215  }
216 
217 
218  if (m_OutputTransform.IsNull()) //default if we didn't manage to instantiate it before
219  {
221  if (firstTransformGiveGeo)
222  {
223  m_OutputProjectionRef =
224  "GEOGCS[\"GCS_WGS_1984\", DATUM[\"D_WGS_1984\", SPHEROID[\"WGS_1984\", 6378137, 298.257223563]], PRIMEM[\"Greenwich\", 0], UNIT[\"Degree\", 0.017453292519943295]]";
225  }
226  otbMsgDevMacro(<< "Output projection set to identity");
227  }
228 
229  m_Transform->SetFirstTransform(m_InputTransform);
230  m_Transform->SetSecondTransform(m_OutputTransform);
231  m_TransformUpToDate = true;
232 
233  //The accuracy information is a simplistic model for now and should be refined
234  if ((inputTransformIsSensor || outputTransformIsSensor))
235  {
236  //Sensor model
237  m_TransformAccuracy = Projection::ESTIMATE;
238  }
239  else if (firstTransformGiveGeo && !outputTransformIsSensor && !outputTransformIsMap)
240  {
241  //The original image was in lon/lat and we did not change anything
242  m_TransformAccuracy = Projection::PRECISE;
243  }
244  else if (!inputTransformIsSensor && !outputTransformIsSensor && !inputTransformIsMap && !outputTransformIsMap)
245  {
246  //no transform
247  m_TransformAccuracy = Projection::UNKNOWN;
248  }
249  else
250  {
251  m_TransformAccuracy = Projection::PRECISE;
252  }
253 
254 }
255 
256 template<class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions>
259 ::TransformPoint(const InputPointType& point) const
260 {
261  InputPointType inputPoint = point;
262 
263  // Apply input origin/spacing
264  inputPoint[0] = inputPoint[0] * m_InputSpacing[0] + m_InputOrigin[0];
265  inputPoint[1] = inputPoint[1] * m_InputSpacing[1] + m_InputOrigin[1];
266 
267  // Transform point
268  OutputPointType outputPoint;
269  outputPoint = this->GetTransform()->TransformPoint(inputPoint);
270 
271  // Apply output origin/spacing
272  outputPoint[0] = (outputPoint[0] - m_OutputOrigin[0]) / m_OutputSpacing[0];
273  outputPoint[1] = (outputPoint[1] - m_OutputOrigin[1]) / m_OutputSpacing[1];
274 
275 // otbMsgDevMacro("GenericRSTransform: " << point << " -> " << outputPoint);
276 
277  return outputPoint;
278 }
279 
280 template<class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions>
281 bool
283 ::GetInverse(Self * inverseTransform) const
284 {
285  // Test the inverseTransform pointer
286  if (inverseTransform == NULL)
287  {
288  return false;
289  }
290 
291  // Swich projection refs
292  inverseTransform->SetInputProjectionRef(m_OutputProjectionRef);
293  inverseTransform->SetOutputProjectionRef(m_InputProjectionRef);
294 
295  // Switch keywordlists
296  inverseTransform->SetInputKeywordList(m_OutputKeywordList);
297  inverseTransform->SetOutputKeywordList(m_InputKeywordList);
298 
299  // Switch dictionnaries
300  inverseTransform->SetInputDictionary(m_OutputDictionary);
301  inverseTransform->SetOutputDictionary(m_InputDictionary);
302 
303  // Switch spacings
304  inverseTransform->SetInputSpacing(m_OutputSpacing);
305  inverseTransform->SetOutputSpacing(m_InputSpacing);
306 
307  // Switch origins
308  inverseTransform->SetInputOrigin(m_OutputOrigin);
309  inverseTransform->SetOutputOrigin(m_InputOrigin);
310 
311  // Instantiate transform
312  inverseTransform->InstanciateTransform();
313 
314  return true;
315 }
316 
317 template<class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions>
319 ::InverseTransformBasePointer
322 {
323  Self * inverseTransform = Self::New();
324 
325  bool success = this->GetInverse(inverseTransform);
326 
327  if (!success)
328  {
329  itkExceptionMacro(<< "Failed to create inverse transform");
330  }
331 
332  return inverseTransform;
333 }
334 
335 template<class TScalarType, unsigned int NInputDimensions, unsigned int NOutputDimensions>
336 void
338 ::PrintSelf(std::ostream& os, itk::Indent indent) const
339 {
340  Superclass::PrintSelf(os, indent);
341  os << indent << "Up to date: " << m_TransformUpToDate << std::endl;
342  if (m_TransformUpToDate)
343  {
344  os << indent << "Input transform: "<< std::endl;
345  m_InputTransform->Print(os, indent.GetNextIndent());
346  os << indent << "Output transform: " << std::endl;
347  m_OutputTransform->Print(os, indent.GetNextIndent());
348  }
349  else
350  {
351  os << indent << "Input transform: NULL" << std::endl;
352  os << indent << "Output transform: NULL" << std::endl;
353  }
354  os << indent << "Accuracy: "
355  << (m_TransformAccuracy == Projection::PRECISE ?
356  "PRECISE" : (m_TransformAccuracy == Projection::ESTIMATE ?
357  "ESTIMATE" : "UNKNOWN")) << std::endl;
358 }
359 
360 } // namespace otb
361 
362 #endif

Generated at Sat Mar 8 2014 15:57:17 for Orfeo Toolbox with doxygen 1.8.3.1