OTB  7.2.0
Orfeo Toolbox
otbImageToGenericRSOutputParameters.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2020 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 otbImageToGenericRSOutputParameters_hxx
22 #define otbImageToGenericRSOutputParameters_hxx
23 
24 
26 #include "itkMacro.h"
27 #include "itkContinuousIndex.h"
28 
29 namespace otb
30 {
31 
32 template <class TImage>
34 {
35  m_Transform = GenericRSTransformType::New();
36  m_ForceSpacing = false;
37  m_ForceSize = false;
38  m_EstimateIsotropicSpacing = false;
39 }
40 
44 template <class TImage>
46 {
47  // Do some checks : exceptions if Null or empty projectionRef and empty keywordlist
48  if (m_Input.IsNull())
49  itkExceptionMacro(<< "The input is null , please set a non null input image");
50 
51  if (m_Input->GetProjectionRef().empty() && m_Input->GetImageKeywordlist().GetSize() == 0)
52  itkExceptionMacro(<< "No information in the metadata, please set an image with non empty metadata");
53 
54  // First Call to UpdateTransform : Initialize with the input image
55  // information
56  this->UpdateTransform();
57 
58  // Estimate the Output image Extent
59  this->EstimateOutputImageExtent();
60 
61  // Estimate the Output Spacing
62  if (!m_ForceSpacing)
63  this->EstimateOutputSpacing();
64 
65  // Finally Estimate the Output Size
66  this->EstimateOutputSize();
67 
68  // Estimate the Output Origin
69  this->EstimateOutputOrigin();
70 }
71 
72 
73 template <class TImage>
75 {
76  m_Transform->SetOutputDictionary(this->GetInput()->GetMetaDataDictionary());
77  m_Transform->SetOutputProjectionRef(this->GetInput()->GetProjectionRef());
78  m_Transform->SetOutputKeywordList(this->GetInput()->GetImageKeywordlist());
79  m_Transform->InstantiateTransform();
80 }
81 
82 
87 template <class TImage>
89 {
90  // Get the inverse transform again : used later
91  GenericRSTransformPointerType invTransform = GenericRSTransformType::New();
92  m_Transform->GetInverse(invTransform);
94 
95  // Compute the 4 corners in the cartographic coordinate system
96  std::vector<itk::ContinuousIndex<double, 2>> vindex;
97  std::vector<PointType> voutput;
98 
99  itk::ContinuousIndex<double, 2> index1(m_Input->GetLargestPossibleRegion().GetIndex());
100  index1[0] += -0.5;
101  index1[1] += -0.5;
102  itk::ContinuousIndex<double, 2> index2(index1);
103  itk::ContinuousIndex<double, 2> index3(index1);
104  itk::ContinuousIndex<double, 2> index4(index1);
105 
106  // Image size
107  SizeType size = m_Input->GetLargestPossibleRegion().GetSize();
108 
109  // project the 4 corners
110  index2[0] += size[0];
111  index3[0] += size[0];
112  index3[1] += size[1];
113  index4[1] += size[1];
114 
115  vindex.push_back(index1);
116  vindex.push_back(index2);
117  vindex.push_back(index3);
118  vindex.push_back(index4);
119 
120  for (unsigned int i = 0; i < vindex.size(); ++i)
121  {
122  PointType physicalPoint;
123  m_Input->TransformContinuousIndexToPhysicalPoint(vindex[i], physicalPoint);
124  voutput.push_back(invTransform->TransformPoint(physicalPoint));
125  }
126 
127  // Compute the boundaries
128  double minX = voutput[0][0];
129  double maxX = voutput[0][0];
130  double minY = voutput[0][1];
131  double maxY = voutput[0][1];
132 
133  for (unsigned int i = 0; i < voutput.size(); ++i)
134  {
135  // Origins
136  if (minX > voutput[i][0])
137  minX = voutput[i][0];
138  if (minY > voutput[i][1])
139  minY = voutput[i][1];
140 
141  // Sizes
142  if (maxX < voutput[i][0])
143  maxX = voutput[i][0];
144 
145  if (maxY < voutput[i][1])
146  maxY = voutput[i][1];
147  }
148 
149  // Edit the output image extent type
150  m_OutputExtent.maxX = maxX;
151  m_OutputExtent.minX = minX;
152  m_OutputExtent.maxY = maxY;
153  m_OutputExtent.minY = minY;
154 }
155 
156 
161 template <class TImage>
163 {
164  // Set the output origin in carto
165  // projection
166  PointType origin;
167  origin[0] = m_OutputExtent.minX + 0.5 * this->GetOutputSpacing()[0];
168  origin[1] = m_OutputExtent.maxY + 0.5 * this->GetOutputSpacing()[1];
169  this->SetOutputOrigin(origin);
170 }
172 
177 template <class TImage>
179 {
180  // Compute the output size
181  double sizeCartoX = std::abs(m_OutputExtent.maxX - m_OutputExtent.minX);
182  double sizeCartoY = std::abs(m_OutputExtent.minY - m_OutputExtent.maxY);
184 
185  PointType o, oX, oY;
186  o[0] = m_OutputExtent.minX;
187  o[1] = m_OutputExtent.maxY;
188 
189  oX = o;
190  oY = o;
191 
192  oX[0] += sizeCartoX;
193  oY[1] -= sizeCartoY;
194 
195  // Transform back into the input image
196  PointType io = m_Transform->TransformPoint(o);
197  PointType ioX = m_Transform->TransformPoint(oX);
198  PointType ioY = m_Transform->TransformPoint(oY);
199 
200  // Transform to indices
201  IndexType ioIndex, ioXIndex, ioYIndex;
202  m_Input->TransformPhysicalPointToIndex(io, ioIndex);
203  m_Input->TransformPhysicalPointToIndex(ioX, ioXIndex);
204  m_Input->TransformPhysicalPointToIndex(ioY, ioYIndex);
205 
206  // Evaluate Ox and Oy length in number of pixels
207  double OxLength, OyLength;
208 
209  OxLength = std::sqrt(std::pow((double)ioIndex[0] - (double)ioXIndex[0], 2) + std::pow((double)ioIndex[1] - (double)ioXIndex[1], 2));
210 
211  OyLength = std::sqrt(std::pow((double)ioIndex[0] - (double)ioYIndex[0], 2) + std::pow((double)ioIndex[1] - (double)ioYIndex[1], 2));
212 
213  // Evaluate spacing
214  SpacingType outputSpacing;
215 
216 
217  if (m_EstimateIsotropicSpacing)
218  {
219  double isotropicSpacing = std::min(sizeCartoX / OxLength, sizeCartoY / OyLength);
220  outputSpacing[0] = isotropicSpacing;
221  outputSpacing[1] = -isotropicSpacing;
222  }
223  else
224  {
225  outputSpacing[0] = sizeCartoX / OxLength;
226  outputSpacing[1] = -sizeCartoY / OyLength;
227  }
228 
229  this->SetOutputSpacing(outputSpacing);
230 }
231 
236 template <class TImage>
238 {
239  // Compute the output size
240  double sizeCartoX = std::abs(m_OutputExtent.maxX - m_OutputExtent.minX);
241  double sizeCartoY = std::abs(m_OutputExtent.minY - m_OutputExtent.maxY);
243 
244  // Evaluate output size
245  SizeType outputSize;
246  outputSize[0] = static_cast<unsigned int>(std::floor(std::abs(sizeCartoX / this->GetOutputSpacing()[0])));
247  outputSize[1] = static_cast<unsigned int>(std::floor(std::abs(sizeCartoY / this->GetOutputSpacing()[1])));
248 
249  // if ForceSizeTo used don't update the output size with the value
250  // computed : the value is computed to update the spacing knowing
251  // the forced size and the computed one.
252  if (!m_ForceSize)
253  this->SetOutputSize(outputSize);
254  else
255  {
256  // Compute the spacing knowing the size
257  SpacingType outputSpacing;
258  outputSpacing[0] = this->GetOutputSpacing()[0] * outputSize[0] / this->GetOutputSize()[0];
259  outputSpacing[1] = this->GetOutputSpacing()[1] * outputSize[1] / this->GetOutputSize()[1];
260  this->SetOutputSpacing(outputSpacing);
261  }
262  this->UpdateTransform();
263 }
264 }
265 
266 #endif
GenericRSTransformType::Pointer GenericRSTransformPointerType
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.