OTB  9.0.0
Orfeo Toolbox
otbGenericRSResampleImageFilter.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 otbGenericRSResampleImageFilter_hxx
22 #define otbGenericRSResampleImageFilter_hxx
23 
25 
26 #include "itkMetaDataObject.h"
27 
28 #include "itkProgressAccumulator.h"
29 
30 #include "itkPoint.h"
31 
32 #include "ogr_spatialref.h"
33 #include "cpl_conv.h"
34 
35 #include "otbSpatialReference.h"
37 
38 namespace otb
39 {
40 
41 template <class TInputImage, class TOutputImage>
43  :
44  // flags initialization
45  m_EstimateInputRpcModel(false),
46  m_EstimateOutputRpcModel(false),
47  m_RpcEstimationUpdated(false),
48 
49  // internal filters instantiation
50  m_Resampler(ResamplerType::New()),
51  m_InputRpcEstimator(InputRpcModelEstimatorType::New()),
52  m_OutputRpcEstimator(OutputRpcModelEstimatorType::New()),
53  m_Transform(GenericRSTransformType::New())
54 {
58 }
59 
60 template <class TInputImage, class TOutputImage>
62 {
63  // Set up progress reporting
64  typename itk::ProgressAccumulator::Pointer progress = itk::ProgressAccumulator::New();
65  progress->SetMiniPipelineFilter(this);
66  progress->RegisterInternalFilter(m_Resampler, 1.f);
67 
68  m_Resampler->GraftOutput(this->GetOutput());
69  m_Resampler->UpdateOutputData(m_Resampler->GetOutput());
70  this->GraftOutput(m_Resampler->GetOutput());
71 }
72 
73 
79 template <class TInputImage, class TOutputImage>
81 {
82  // Estimate the output rpc Model if needed
83  if (m_EstimateOutputRpcModel)
84  this->EstimateOutputRpcModel();
85 
86  // Estimate the input rpc model if it is needed
87  if (m_EstimateInputRpcModel && !m_RpcEstimationUpdated)
88  {
89  this->EstimateInputRpcModel();
90  }
91 
92  // Instantiate the RS transform
93  this->UpdateTransform();
94 
95  m_Resampler->SetInput(this->GetInput());
96  m_Resampler->SetTransform(m_Transform);
97  m_Resampler->SetDisplacementFieldSpacing(this->GetDisplacementFieldSpacing());
98  m_Resampler->GraftOutput(this->GetOutput());
99  m_Resampler->UpdateOutputInformation();
100  this->GraftOutput(m_Resampler->GetOutput());
101 
102  // Encapsulate output projRef and metadata
103  if (this->GetOutputImageMetadata() != nullptr)
104  this->GetOutput()->m_Imd.Merge(*(this->GetOutputImageMetadata()));
105  this->GetOutput()->m_Imd.Add(MDGeom::ProjectionWKT, this->GetOutputProjectionRef());
106 }
107 
111 template <class TInputImage, class TOutputImage>
113 {
114  // Temp image : not allocated but with the same metadata than the
115  // output
116  typename OutputImageType::Pointer tempPtr = OutputImageType::New();
117 
118  typename OutputImageType::RegionType region;
119  region.SetSize(this->GetOutputSize());
120  region.SetIndex(this->GetOutputStartIndex());
121  tempPtr->SetRegions(region);
122 
123  // Encapsulate the output metadata in the temp image
124  tempPtr->m_Imd.Add(MDGeom::ProjectionWKT, this->GetOutputProjectionRef());
125  tempPtr->SetImageMetadata(*(this->GetOutputImageMetadata()));
126 
127  // Estimate the rpc model from the temp image
128  m_OutputRpcEstimator->SetInput(tempPtr);
129  m_OutputRpcEstimator->UpdateOutputInformation();
130 
131  // Fill the transform with the right metadata
132  m_Transform->SetInputImageMetadata(&(m_OutputRpcEstimator->GetOutput()->GetImageMetadata()));
133 }
134 
139 template <class TInputImage, class TOutputImage>
141 {
142  if (!m_EstimateInputRpcModel)
143  {
144  m_Transform->SetOutputImageMetadata(&(this->GetInput()->GetImageMetadata()));
145  m_Transform->SetOutputProjectionRef(this->GetInput()->GetProjectionRef());
146  }
147  m_Transform->InstantiateTransform();
148 }
150 
151 template <class TInputImage, class TOutputImage>
153 {
154  if (this->m_Updating)
155  return;
156 
157  // Retrieve output requested region
158  m_Resampler->GetOutput()->SetRequestedRegion(output);
159  m_Resampler->GetOutput()->PropagateRequestedRegion();
160 }
161 
162 
168 template <class TInputImage, class TOutputImage>
170 {
171  // Temp image : not allocated but with the sampe metadata as the
172  // output
173  typename InputImageType::Pointer tempPtr = InputImageType::New();
174  tempPtr->SetRegions(this->GetInput()->GetLargestPossibleRegion());
175  tempPtr->CopyInformation(this->GetInput());
177 
178  // Estimate the rpc model with the temp image
179  m_InputRpcEstimator->SetInput(tempPtr);
180  m_InputRpcEstimator->UpdateOutputInformation();
181 
182  // setup the transform with the estimated RPC model
183  m_Transform->SetOutputImageMetadata(&(m_InputRpcEstimator->GetOutput()->GetImageMetadata()));
184 
185  // Update the flag for input rpcEstimation in order to not compute
186  // the rpc model for each stream
187  m_RpcEstimationUpdated = true;
188 }
189 
194 template <class TInputImage, class TOutputImage>
196 {
197  const InputImageType* src = dynamic_cast<const InputImageType*>(image);
198 
199  this->SetOutputOrigin(src->GetOrigin());
200  this->SetOutputSpacing(src->GetSignedSpacing());
201  this->SetOutputStartIndex(src->GetLargestPossibleRegion().GetIndex());
202  this->SetOutputSize(src->GetLargestPossibleRegion().GetSize());
203  this->SetOutputProjectionRef(src->GetProjectionRef());
204  this->GetOutput()->SetImageMetadata(src->GetImageMetadata());
205 }
206 
211 template <class TInputImage, class TOutputImage>
212 template <class TImageType>
214 {
215  this->SetOutputOrigin(image->GetOrigin());
216  this->SetOutputSpacing(image->GetSignedSpacing());
217  this->SetOutputStartIndex(image->GetLargestPossibleRegion().GetIndex());
218  this->SetOutputSize(image->GetLargestPossibleRegion().GetSize());
219  this->SetOutputProjectionRef(image->GetProjectionRef());
220  this->GetOutput()->SetImageMetadata(image->GetImageMetadata());
221 }
223 
230 template <class TInputImage, class TOutputImage>
232 {
233  // Get the input Image
234  const InputImageType* input = this->GetInput();
235 
236  // Update the transform with input information
237  // Done here because the transform is not instantiated
238  // yet
239  this->UpdateTransform();
240 
241  // Needed variable
242  std::string projectionRef;
243  // The inverse transform is need here
244  GenericRSTransformPointerType invTransform = GenericRSTransformType::New();
245  m_Transform->GetInverse(invTransform);
246 
247  if (map == "UTM")
248  {
249  // Build the UTM transform : Need the zone & the hemisphere
250  // For this we us the geographic coordinate of the input UL corner
251  typedef itk::Point<double, 2> GeoPointType;
252 
253  // get the utm zone and hemisphere using the input UL corner
254  // geographic coordinates
255  typename InputImageType::PointType pSrc;
256  IndexType index;
257  GeoPointType geoPoint;
258  index[0] = input->GetLargestPossibleRegion().GetIndex()[0];
259  index[1] = input->GetLargestPossibleRegion().GetIndex()[1];
260  input->TransformIndexToPhysicalPoint(index, pSrc);
261 
262  // The first transform of the inverse transform : input -> WGS84
263  geoPoint = invTransform->GetTransform()->GetFirstTransform()->TransformPoint(pSrc);
264 
265  // Guess the zone and the hemisphere
266  unsigned int zone(0);
268  otb::SpatialReference::UTMFromGeoPoint(geoPoint[0], geoPoint[1], zone, hem);
269 
271 
272  projectionRef = oSRS.ToWkt();
273  }
274  else if (map == "WGS84")
275  {
276  projectionRef = otb::SpatialReference::FromWGS84().ToWkt(); // WGS84
277  }
278  else
279  {
280  itkExceptionMacro("The output map " << map << "is not supported, please try UTM or WGS84");
281  }
282 
283  // Compute the output parameters
284  typedef otb::ImageToGenericRSOutputParameters<InputImageType> OutputParametersEstimatorType;
285  typename OutputParametersEstimatorType::Pointer genericRSEstimator = OutputParametersEstimatorType::New();
286 
287  genericRSEstimator->SetInput(input);
288  genericRSEstimator->SetOutputProjectionRef(projectionRef);
289  genericRSEstimator->ForceSpacingTo(spacing);
290  genericRSEstimator->Compute();
291 
292  // Update the Output Parameters
293  this->SetOutputProjectionRef(projectionRef);
294  this->SetOutputOrigin(genericRSEstimator->GetOutputOrigin());
295  this->SetOutputSpacing(genericRSEstimator->GetOutputSpacing());
296  this->SetOutputSize(genericRSEstimator->GetOutputSize());
297  this->UpdateTransform();
298 }
299 
306 template <class TInputImage, class TOutputImage>
308 {
309  const InputImageType* input = this->GetInput();
310 
311  // Compute the output parameters
312  typedef otb::ImageToGenericRSOutputParameters<InputImageType> OutputParametersEstimatorType;
313  typename OutputParametersEstimatorType::Pointer genericRSEstimator = OutputParametersEstimatorType::New();
314 
315  genericRSEstimator->SetInput(input);
316  genericRSEstimator->SetOutputProjectionRef(projectionRef);
317  genericRSEstimator->Compute();
318 
319  // Update the Output Parameters
320  this->SetOutputProjectionRef(projectionRef);
321  this->SetOutputOrigin(genericRSEstimator->GetOutputOrigin());
322  this->SetOutputSpacing(genericRSEstimator->GetOutputSpacing());
323  this->SetOutputSize(genericRSEstimator->GetOutputSize());
324  this->UpdateTransform();
325 }
326 
327 template <class TInputImage, class TOutputImage>
328 void GenericRSResampleImageFilter<TInputImage, TOutputImage>::PrintSelf(std::ostream& os, itk::Indent indent) const
329 {
330  Superclass::PrintSelf(os, indent);
331  os << indent << "EstimateInputRpcModel:" << (m_EstimateInputRpcModel ? "On" : "Off") << std::endl;
332  os << indent << "EstimateOutputRpcModel:" << (m_EstimateOutputRpcModel ? "On" : "Off") << std::endl;
333  os << indent << "RpcEstimationUpdated:" << (m_RpcEstimationUpdated ? "True" : "False") << std::endl;
334  os << indent << "OutputOrigin: " << m_Resampler->GetOutputOrigin() << std::endl;
335  os << indent << "OutputSpacing: " << m_Resampler->GetOutputSpacing() << std::endl;
336  os << indent << "OutputStartIndex: " << m_Resampler->GetOutputStartIndex() << std::endl;
337  os << indent << "OutputSize: " << m_Resampler->GetOutputSize() << std::endl;
338  os << indent << "GenericRSTransform: " << std::endl;
339  m_Transform->Print(os, indent.GetNextIndent());
340 }
341 }
342 #endif
otb::GenericRSResampleImageFilter::SetOutputParametersFromImage
void SetOutputParametersFromImage(const ImageBaseType *image)
Definition: otbGenericRSResampleImageFilter.hxx:195
otb::PhysicalToRPCSensorModelImageFilter
This filter estimates a RPC sensor models from a physical model.
Definition: otbPhysicalToRPCSensorModelImageFilter.h:62
otb::GenericRSResampleImageFilter::GenericRSTransformPointerType
GenericRSTransformType::Pointer GenericRSTransformPointerType
Definition: otbGenericRSResampleImageFilter.h:100
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:32
otbImageToGenericRSOutputParameters.h
otb::ImageToGenericRSOutputParameters
This class is a helper class to estimate the output parameters of an image after projection in a targ...
Definition: otbImageToGenericRSOutputParameters.h:57
otb::SpatialReference::UTMFromGeoPoint
static void UTMFromGeoPoint(double lon, double lat, unsigned int &zone, hemisphere &hem)
otb::SpatialReference
This class is a wrapper around OGRSpatialReference.
Definition: otbSpatialReference.h:79
otb::GenericRSResampleImageFilter::IndexType
ResamplerType::IndexType IndexType
Definition: otbGenericRSResampleImageFilter.h:84
otbSpatialReference.h
otb::StreamingResampleImageFilter
This class is a composite filter resampling an input image by setting a transform....
Definition: otbStreamingResampleImageFilter.h:59
otb::GenericRSTransform
This is the class to handle generic remote sensing transform.
Definition: otbGenericRSTransform.h:57
otb::GenericRSResampleImageFilter::SetDisplacementFilterNumberOfThreads
void SetDisplacementFilterNumberOfThreads(unsigned int nbThread)
Definition: otbGenericRSResampleImageFilter.h:277
otb::SpatialReference::ToWkt
std::string ToWkt() const
otb::GenericRSResampleImageFilter::SetOutputParametersFromMap
void SetOutputParametersFromMap(const std::string &map, const SpacingType &spacing)
Definition: otbGenericRSResampleImageFilter.hxx:231
otb::GenericRSResampleImageFilter::SpacingType
ResamplerType::SpacingType SpacingType
Definition: otbGenericRSResampleImageFilter.h:82
otbGenericRSResampleImageFilter.h
otb::GenericRSResampleImageFilter::EstimateOutputRpcModel
void EstimateOutputRpcModel()
Definition: otbGenericRSResampleImageFilter.hxx:112
otb::GenericRSResampleImageFilter::GenerateOutputInformation
void GenerateOutputInformation() override
Definition: otbGenericRSResampleImageFilter.hxx:80
otb::SpatialReference::FromUTM
static SpatialReference FromUTM(unsigned int zone, hemisphere hem)
otb::GenericRSResampleImageFilter::ImageBaseType
itk::ImageBase< OutputImageType::ImageDimension > ImageBaseType
Definition: otbGenericRSResampleImageFilter.h:102
otb::GenericRSResampleImageFilter::UpdateTransform
virtual void UpdateTransform()
Definition: otbGenericRSResampleImageFilter.hxx:140
otb::GenericRSResampleImageFilter::GenericRSResampleImageFilter
GenericRSResampleImageFilter()
Definition: otbGenericRSResampleImageFilter.hxx:42
otb::SpatialReference::FromWGS84
static SpatialReference FromWGS84()
otb::GenericRSResampleImageFilter::PrintSelf
void PrintSelf(std::ostream &os, itk::Indent indent) const override
Definition: otbGenericRSResampleImageFilter.hxx:328
otb::GenericRSResampleImageFilter::GenerateData
void GenerateData() override
Definition: otbGenericRSResampleImageFilter.hxx:61
otb::GenericRSResampleImageFilter::EstimateInputRpcModel
void EstimateInputRpcModel()
Definition: otbGenericRSResampleImageFilter.hxx:169
otb::GenericRSResampleImageFilter::PropagateRequestedRegion
void PropagateRequestedRegion(itk::DataObject *output) override
Definition: otbGenericRSResampleImageFilter.hxx:152
otb::MDGeom::ProjectionWKT
@ ProjectionWKT
otb::GenericRSResampleImageFilter::InputImageType
TInputImage InputImageType
Definition: otbGenericRSResampleImageFilter.h:69
otb::SpatialReference::hemisphere
hemisphere
Definition: otbSpatialReference.h:127