OTB  6.7.0
Orfeo Toolbox
otbImageToGenericRSOutputParameters.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2019 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 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>
45 void
48 {
49  // Do some checks : exceptions if Null or empty projectionRef and empty keywordlist
50  if(m_Input.IsNull())
51  itkExceptionMacro(<<"The input is null , please set a non null input image");
52 
53  if(m_Input->GetProjectionRef().empty() && m_Input->GetImageKeywordlist().GetSize() == 0)
54  itkExceptionMacro(<<"No information in the metadata, please set an image with non empty metadata");
55 
56  // First Call to UpdateTransform : Initialize with the input image
57  // information
58  this->UpdateTransform();
59 
60  // Estimate the Output image Extent
61  this->EstimateOutputImageExtent();
62 
63  // Estimate the Output Spacing
64  if(!m_ForceSpacing)
65  this->EstimateOutputSpacing();
66 
67  // Finally Estimate the Output Size
68  this->EstimateOutputSize();
69 
70  // Estimate the Output Origin
71  this->EstimateOutputOrigin();
72 }
73 
74 
75 template<class TImage>
76 void
79 {
80  m_Transform->SetOutputDictionary(this->GetInput()->GetMetaDataDictionary());
81  m_Transform->SetOutputProjectionRef(this->GetInput()->GetProjectionRef());
82  m_Transform->SetOutputKeywordList(this->GetInput()->GetImageKeywordlist());
83  m_Transform->InstantiateTransform();
84 }
85 
86 
91 template <class TImage>
92 void
95 {
96  // Get the inverse transform again : used later
97  GenericRSTransformPointerType invTransform = GenericRSTransformType::New();
98  m_Transform->GetInverse(invTransform);
100 
101  // Compute the 4 corners in the cartographic coordinate system
102  std::vector< itk::ContinuousIndex<double,2> > vindex;
103  std::vector<PointType> voutput;
104 
105  itk::ContinuousIndex<double,2> index1(m_Input->GetLargestPossibleRegion().GetIndex());
106  index1[0] += -0.5;
107  index1[1] += -0.5;
108  itk::ContinuousIndex<double,2> index2(index1);
109  itk::ContinuousIndex<double,2> index3(index1);
110  itk::ContinuousIndex<double,2> index4(index1);
111 
112  // Image size
113  SizeType size = m_Input->GetLargestPossibleRegion().GetSize();
114 
115  // project the 4 corners
116  index2[0] += size[0];
117  index3[0] += size[0];
118  index3[1] += size[1];
119  index4[1] += size[1];
120 
121  vindex.push_back(index1);
122  vindex.push_back(index2);
123  vindex.push_back(index3);
124  vindex.push_back(index4);
125 
126  for (unsigned int i = 0; i < vindex.size(); ++i)
127  {
128  PointType physicalPoint;
129  m_Input->TransformContinuousIndexToPhysicalPoint(vindex[i], physicalPoint);
130  voutput.push_back(invTransform->TransformPoint(physicalPoint));
131  }
132 
133  // Compute the boundaries
134  double minX = voutput[0][0];
135  double maxX = voutput[0][0];
136  double minY = voutput[0][1];
137  double maxY = voutput[0][1];
138 
139  for (unsigned int i = 0; i < voutput.size(); ++i)
140  {
141  // Origins
142  if (minX > voutput[i][0])
143  minX = voutput[i][0];
144  if (minY > voutput[i][1])
145  minY = voutput[i][1];
146 
147  // Sizes
148  if (maxX < voutput[i][0])
149  maxX = voutput[i][0];
150 
151  if (maxY < voutput[i][1])
152  maxY = voutput[i][1];
153  }
154 
155  // Edit the output image extent type
156  m_OutputExtent.maxX = maxX;
157  m_OutputExtent.minX = minX;
158  m_OutputExtent.maxY = maxY;
159  m_OutputExtent.minY = minY;
160 }
161 
162 
167 template <class TImage>
168 void
171 {
172  // Set the output origin in carto
173  // projection
174  PointType origin;
175  origin[0] = m_OutputExtent.minX + 0.5 * this->GetOutputSpacing()[0];
176  origin[1] = m_OutputExtent.maxY + 0.5 * this->GetOutputSpacing()[1];
177  this->SetOutputOrigin(origin);
178 }
180 
185 template <class TImage>
186 void
189 {
190  // Compute the output size
191  double sizeCartoX = std::abs(m_OutputExtent.maxX - m_OutputExtent.minX);
192  double sizeCartoY = std::abs(m_OutputExtent.minY - m_OutputExtent.maxY);
194 
195  PointType o, oX, oY;
196  o[0] = m_OutputExtent.minX;
197  o[1] = m_OutputExtent.maxY;
198 
199  oX = o;
200  oY = o;
201 
202  oX[0] += sizeCartoX;
203  oY[1] -= sizeCartoY;
204 
205  // Transform back into the input image
206  PointType io = m_Transform->TransformPoint(o);
207  PointType ioX = m_Transform->TransformPoint(oX);
208  PointType ioY = m_Transform->TransformPoint(oY);
209 
210  // Transform to indices
211  IndexType ioIndex, ioXIndex, ioYIndex;
212  m_Input->TransformPhysicalPointToIndex(io, ioIndex);
213  m_Input->TransformPhysicalPointToIndex(ioX, ioXIndex);
214  m_Input->TransformPhysicalPointToIndex(ioY, ioYIndex);
215 
216  // Evaluate Ox and Oy length in number of pixels
217  double OxLength, OyLength;
218 
219  OxLength = std::sqrt(std::pow((double) ioIndex[0] - (double) ioXIndex[0], 2)
220  + std::pow((double) ioIndex[1] - (double) ioXIndex[1], 2));
221 
222  OyLength = std::sqrt(std::pow((double) ioIndex[0] - (double) ioYIndex[0], 2)
223  + std::pow((double) ioIndex[1] - (double) ioYIndex[1], 2));
224 
225  // Evaluate spacing
226  SpacingType outputSpacing;
227 
228 
229  if(m_EstimateIsotropicSpacing)
230  {
231  double isotropicSpacing = std::min(sizeCartoX / OxLength, sizeCartoY / OyLength);
232  outputSpacing[0] = isotropicSpacing;
233  outputSpacing[1] = -isotropicSpacing;
234  }
235  else
236  {
237  outputSpacing[0] = sizeCartoX / OxLength;
238  outputSpacing[1] = -sizeCartoY / OyLength;
239  }
240 
241  this->SetOutputSpacing(outputSpacing);
242 }
243 
248 template <class TImage>
249 void
252 {
253  // Compute the output size
254  double sizeCartoX = std::abs(m_OutputExtent.maxX - m_OutputExtent.minX);
255  double sizeCartoY = std::abs(m_OutputExtent.minY - m_OutputExtent.maxY);
257 
258  // Evaluate output size
259  SizeType outputSize;
260  outputSize[0] = static_cast<unsigned int>(std::floor(std::abs(sizeCartoX / this->GetOutputSpacing()[0])));
261  outputSize[1] = static_cast<unsigned int>(std::floor(std::abs(sizeCartoY / this->GetOutputSpacing()[1])));
262 
263  // if ForceSizeTo used don't update the output size with the value
264  // computed : the value is computed to update the spacing knowing
265  // the forced size and the computed one.
266  if(!m_ForceSize)
267  this->SetOutputSize(outputSize);
268  else
269  {
270  // Compute the spacing knowing the size
271  SpacingType outputSpacing;
272  outputSpacing[0] = this->GetOutputSpacing()[0] * outputSize[0]/this->GetOutputSize()[0];
273  outputSpacing[1] = this->GetOutputSpacing()[1] * outputSize[1]/this->GetOutputSize()[1];
274  this->SetOutputSpacing(outputSpacing);
275  }
276  this->UpdateTransform();
277 }
278 
279 
280 }
281 
282 #endif