OTB  6.7.0
Orfeo Toolbox
otbStreamingWarpImageFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 1999-2011 Insight Software Consortium
3  * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
4  *
5  * This file is part of Orfeo Toolbox
6  *
7  * https://www.orfeo-toolbox.org/
8  *
9  * Licensed under the Apache License, Version 2.0 (the "License");
10  * you may not use this file except in compliance with the License.
11  * You may obtain a copy of the License at
12  *
13  * http://www.apache.org/licenses/LICENSE-2.0
14  *
15  * Unless required by applicable law or agreed to in writing, software
16  * distributed under the License is distributed on an "AS IS" BASIS,
17  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
18  * See the License for the specific language governing permissions and
19  * limitations under the License.
20  */
21 
22 #ifndef otbStreamingWarpImageFilter_hxx
23 #define otbStreamingWarpImageFilter_hxx
24 
28 #include "itkMetaDataObject.h"
29 #include "otbMetaDataKey.h"
30 
31 namespace otb
32 {
33 
34 template<class TInputImage, class TOutputImage, class TDisplacementField>
37  {
38  // Fill the default maximum displacement
39  m_MaximumDisplacement.Fill(1);
40  m_OutputSignedSpacing = this->Superclass::GetOutputSpacing();
41  }
42 
43 
44 template<class TInputImage, class TOutputImage, class TDisplacementField>
45 void
47 ::SetOutputSpacing( const SpacingType outputSpacing )
48 {
49  m_OutputSignedSpacing = outputSpacing;
50  SpacingType spacing = outputSpacing;
51  typename TInputImage::DirectionType direction =
52  this->GetOutput()->GetDirection();
53  for( unsigned int i = 0 ; i < TInputImage::ImageDimension ; ++i )
54  {
55  if ( spacing[i] < 0 )
56  {
57  if ( direction[i][i] > 0 )
58  {
59  for( unsigned int j = 0 ; j < TInputImage::ImageDimension ; ++j )
60  {
61  direction[j][i] = - direction[j][i];
62  }
63  }
64  spacing[i] = - spacing[i];
65  }
66  }
67  this->Superclass::SetOutputSpacing( spacing );
68  this->Superclass::SetOutputDirection( direction );
69  this->Modified();
70 }
71 
72 template<class TInputImage, class TOutputImage, class TDisplacementField>
73 void
75 ::SetOutputSpacing( const double *values )
76 {
77  SpacingType s;
78  for(unsigned int i = 0; i < TInputImage::ImageDimension; ++i)
79  {
80  s[i] = static_cast< typename SpacingType::ValueType >(values[i]);
81  }
82  this->SetOutputSpacing(s);
83 }
84 
85 template<class TInputImage, class TOutputImage, class TDisplacementField>
86 void
89  {
90  Superclass::GenerateInputRequestedRegion();
91 
92  // Get the input and displacement field pointers
93  InputImageType * inputPtr
94  = const_cast<InputImageType *>(this->GetInput());
95  DisplacementFieldType * displacementPtr
96  = const_cast<DisplacementFieldType*>(this->GetDisplacementField());
97  OutputImageType * outputPtr
98  = this->GetOutput();
99  // Check if the input and the displacement field exist
100  if (!inputPtr || !displacementPtr || !outputPtr)
101  {
102  return;
103  }
104 
105  // Here we are breaking traditional pipeline steps because we need to access the displacement field data
106  // so as to compute the input image requested region
107 
108  // 1) First, evaluate the displacement field requested region corresponding to the output requested region
109  // (Here we suppose that the displacement field and the output image are in the same geometry/map projection)
110  typename OutputImageType::RegionType outputRequestedRegion = outputPtr->GetRequestedRegion();
111  typename OutputImageType::IndexType outIndexStart = outputRequestedRegion.GetIndex();
112  typename OutputImageType::IndexType outIndexEnd;
113  for(unsigned int dim = 0; dim<OutputImageType::ImageDimension; ++dim)
114  outIndexEnd[dim]= outIndexStart[dim] + outputRequestedRegion.GetSize()[dim]-1;
115  typename OutputImageType::PointType outPointStart, outPointEnd;
116  outputPtr->TransformIndexToPhysicalPoint(outIndexStart, outPointStart);
117  outputPtr->TransformIndexToPhysicalPoint(outIndexEnd, outPointEnd);
118 
119  typename DisplacementFieldType::IndexType defIndexStart, defIndexEnd;
120  displacementPtr->TransformPhysicalPointToIndex(outPointStart, defIndexStart);
121  displacementPtr->TransformPhysicalPointToIndex(outPointEnd, defIndexEnd);
122 
123  typename DisplacementFieldType::SizeType defRequestedSize;
124  typename DisplacementFieldType::IndexType defRequestedIndex;
125 
126  for(unsigned int dim = 0; dim<OutputImageType::ImageDimension; ++dim)
127  {
128  defRequestedIndex[dim] = std::min(defIndexStart[dim], defIndexEnd[dim]);
129  defRequestedSize[dim] = std::max(defIndexStart[dim], defIndexEnd[dim]) - defRequestedIndex[dim] + 1;
130  }
131 
132  // Finally, build the displacement field requested region
133  typename DisplacementFieldType::RegionType displacementRequestedRegion;
134  displacementRequestedRegion.SetIndex(defRequestedIndex);
135  displacementRequestedRegion.SetSize(defRequestedSize);
136 
137  // Avoid extrapolation
138  displacementRequestedRegion.PadByRadius(1);
139 
140  // crop the input requested region at the input's largest possible region
141  if (displacementRequestedRegion.Crop(displacementPtr->GetLargestPossibleRegion()))
142  {
143  displacementPtr->SetRequestedRegion(displacementRequestedRegion);
144  }
145  else
146  {
147  // Couldn't crop the region (requested region is outside the largest
148  // possible region). Throw an exception.
149 
150  // store what we tried to request (prior to trying to crop)
151  displacementPtr->SetRequestedRegion(displacementRequestedRegion);
152 
153  // build an exception
154  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
155  e.SetLocation(ITK_LOCATION);
156  e.SetDescription("Requested region is (at least partially) outside the largest possible region of the displacement field.");
157  e.SetDataObject(inputPtr);
158  throw e;
159  }
160 
161  // 2) If we are still there, we have a correct displacement field requested region.
162  // This next step breaks pipeline rule but we need to do it to compute the input requested region,
163  // since it depends on displacement value.
164 
165  // Trigger pipeline update on the displacement field
166 
167  displacementPtr->PropagateRequestedRegion();
168  displacementPtr->UpdateOutputData();
169 
170  // Walk the loaded displacement field to derive maximum and minimum
171  itk::ImageRegionIteratorWithIndex<DisplacementFieldType> defIt(displacementPtr, displacementRequestedRegion);
172  defIt.GoToBegin();
173 
174  typename InputImageType::PointType currentPoint;
175  typename InputImageType::PointType inputStartPoint, inputEndPoint;
176 
177  // Initialise start and end points
178  displacementPtr->TransformIndexToPhysicalPoint(defIt.GetIndex(), currentPoint);
179  for(unsigned int dim = 0; dim<DisplacementFieldType::ImageDimension; ++dim)
180  {
181  currentPoint[dim]+=defIt.Get()[dim];
182  }
183  inputStartPoint = currentPoint;
184  inputEndPoint = currentPoint;
185 
186  ++defIt;
187 
188  // 3) Now Walk the field and compute the physical bounding box
189  while(!defIt.IsAtEnd())
190  {
191  displacementPtr->TransformIndexToPhysicalPoint(defIt.GetIndex(), currentPoint);
192  for(unsigned int dim = 0; dim<DisplacementFieldType::ImageDimension; ++dim)
193  {
194  currentPoint[dim]+=defIt.Get()[dim];
195  if(inputStartPoint[dim] > currentPoint[dim])
196  inputStartPoint[dim] = currentPoint[dim];
197  if(inputEndPoint[dim] < currentPoint[dim])
198  inputEndPoint[dim] = currentPoint[dim];
199  }
200  ++defIt;
201  }
202 
203  // Convert physical bounding box to requested region
204  typename InputImageType::IndexType inputStartIndex, inputEndIndex;
205  inputPtr->TransformPhysicalPointToIndex(inputStartPoint, inputStartIndex);
206  inputPtr->TransformPhysicalPointToIndex(inputEndPoint, inputEndIndex);
207 
208  typename InputImageType::SizeType inputFinalSize;
209  typename InputImageType::IndexType inputFinalIndex;
210 
211  for(unsigned int dim = 0; dim<DisplacementFieldType::ImageDimension; ++dim)
212  {
213  inputFinalIndex[dim] = std::min(inputStartIndex[dim], inputEndIndex[dim]);
214  inputFinalSize[dim] = std::max(inputStartIndex[dim], inputEndIndex[dim])-inputFinalIndex[dim]+1;
215  }
216 
217  typename InputImageType::RegionType inputRequestedRegion;
218  inputRequestedRegion.SetIndex(inputFinalIndex);
219  inputRequestedRegion.SetSize(inputFinalSize);
220 
221  // Compute the padding due to the interpolator
222  unsigned int interpolatorRadius =
224 
225  // pad the input requested region by the operator radius
226  inputRequestedRegion.PadByRadius(interpolatorRadius);
227 
228  // crop the input requested region at the input's largest possible region
229  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
230  {
231  inputPtr->SetRequestedRegion(inputRequestedRegion);
232  }
233  else
234  {
235 
236  inputFinalSize.Fill(0);
237  inputRequestedRegion.SetSize(inputFinalSize);
238  inputFinalIndex.Fill(0);
239  inputRequestedRegion.SetIndex(inputFinalIndex);
240 
241  // store what we tried to request (prior to trying to crop)
242  inputPtr->SetRequestedRegion(inputRequestedRegion);
243 
244 // // build an exception
245 // itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
246 // e.SetLocation(ITK_LOCATION);
247 // e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
248 // e.SetDataObject(inputPtr);
249 // throw e;
250  }
251  }
252 
253 
254 template<class TInputImage, class TOutputImage, class TDisplacementField>
255 void
258 {
259  Superclass::GenerateOutputInformation();
260 
261  // Set the NoData flag to the edge padding value
262  itk::MetaDataDictionary& dict = this->GetOutput()->GetMetaDataDictionary();
263  std::vector<bool> noDataValueAvailable;
264  bool ret = itk::ExposeMetaData<std::vector<bool> >(dict,MetaDataKey::NoDataValueAvailable,noDataValueAvailable);
265  if (!ret)
266  {
267  noDataValueAvailable.resize(this->GetOutput()->GetNumberOfComponentsPerPixel(),false);
268  }
269  std::vector<double> noDataValue;
270  ret = itk::ExposeMetaData<std::vector<double> >(dict,MetaDataKey::NoDataValue,noDataValue);
271  if (!ret)
272  {
273  noDataValue.resize(this->GetOutput()->GetNumberOfComponentsPerPixel(),0.0);
274  }
275  PixelType edgePadding = this->GetEdgePaddingValue();
276  if (itk::NumericTraits<PixelType>::GetLength(edgePadding) != this->GetOutput()->GetNumberOfComponentsPerPixel())
277  {
278  itk::NumericTraits<PixelType>::SetLength(edgePadding,this->GetOutput()->GetNumberOfComponentsPerPixel());
279  this->SetEdgePaddingValue(edgePadding);
280  }
281  for (unsigned int i=0; i<noDataValueAvailable.size() ; ++i)
282  {
283  if (!noDataValueAvailable[i])
284  {
285  noDataValueAvailable[i] = true;
286  noDataValue[i] = itk::DefaultConvertPixelTraits<PixelType>::GetNthComponent(i,edgePadding);
287  }
288  }
289  itk::EncapsulateMetaData<std::vector<bool> >(dict,MetaDataKey::NoDataValueAvailable,noDataValueAvailable);
290  itk::EncapsulateMetaData<std::vector<double> >(dict,MetaDataKey::NoDataValue,noDataValue);
291 }
292 
293 template<class TInputImage, class TOutputImage, class TDisplacementField>
294 void
297  const OutputImageRegionType& outputRegionForThread,
298  itk::ThreadIdType threadId )
299  {
300  // the superclass itk::WarpImageFilter is doing the actual warping
301  Superclass::ThreadedGenerateData(outputRegionForThread,threadId);
302 
303  // second pass on the thread region to mask pixels outside the displacement grid
304  const PixelType paddingValue = this->GetEdgePaddingValue();
305  OutputImagePointerType outputPtr = this->GetOutput();
306 
307  // ITK 4.13 fix const correctness of GetDisplacementField.
308  // Related commit in ITK: https://github.com/InsightSoftwareConsortium/ITK/commit/0070848b91baf69f04893bc3ce85bcf110c3c63a
309 
310  // DisplacementFieldPointerType fieldPtr = this->GetDisplacementField();
311  const DisplacementFieldType * fieldPtr = this->GetDisplacementField();
312 
313 
314  DisplacementFieldRegionType defRegion = fieldPtr->GetLargestPossibleRegion();
315 
317  outputPtr, outputRegionForThread );
318  IndexType currentIndex;
319  PointType currentPoint;
321 
322  while(!outputIt.IsAtEnd())
323  {
324  // get the output image index
325  currentIndex = outputIt.GetIndex();
326  outputPtr->TransformIndexToPhysicalPoint(currentIndex,currentPoint);
327  fieldPtr->TransformPhysicalPointToContinuousIndex(currentPoint,contiIndex);
328 
329  for (unsigned int dim = 0; dim<DisplacementFieldType::ImageDimension; ++dim)
330  {
331  if (contiIndex[dim] < static_cast<double>(defRegion.GetIndex(dim)) ||
332  contiIndex[dim] > static_cast<double>(defRegion.GetIndex(dim)+defRegion.GetSize(dim)-1))
333  {
334  outputIt.Set(paddingValue);
335  break;
336  }
337  }
338  ++outputIt;
339  }
340  }
341 
342 template<class TInputImage, class TOutputImage, class TDisplacementField>
343 void
345 ::PrintSelf(std::ostream& os, itk::Indent indent) const
346  {
347  Superclass::PrintSelf(os, indent);
348  os << indent << "Maximum displacement: " << m_MaximumDisplacement << std::endl;
349  }
350 
351 } // end namespace otb
352 #endif
itk::Size< Monteverdi_DIMENSION > SizeType
Definition: mvdTypes.h:137
static unsigned int CalculateNeededRadiusForInterpolator(const InterpolationType *interpolator)
void SetDataObject(DataObject *dobj)
PixelType Get(void) const
virtual void SetDescription(const std::string &s)
virtual void SetLocation(const std::string &s)
OTBOSSIMAdapters_EXPORT char const * NoDataValue
DisplacementFieldType::RegionType DisplacementFieldRegionType
void Set(const PixelType &value) const
OutputImageType::SpacingType SpacingType
const IndexType & GetIndex() const
OutputImageType::Pointer OutputImagePointerType
static ComponentType GetNthComponent(int c, const PixelType &pixel)
void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, itk::ThreadIdType threadId) override
itk::Index< Monteverdi_DIMENSION > IndexType
Definition: mvdTypes.h:133
OTBOSSIMAdapters_EXPORT char const * NoDataValueAvailable
TInputImage InputImageType
static void SetLength(T &m, const unsigned int s)
unsigned int ThreadIdType
OutputImageType::PixelType PixelType
OutputImageType::RegionType OutputImageRegionType
OutputImageType::IndexType IndexType
void SetOutputSpacing(const SpacingType OutputSpacing) override
TOutputImage OutputImageType
void PrintSelf(std::ostream &os, itk::Indent indent) const override
VectorImageType::PointType PointType
Definition: mvdTypes.h:189
OutputImageType::PointType PointType