OTB  6.7.0
Orfeo Toolbox
otbBandMathImageFilter.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 otbBandMathImageFilter_hxx
23 #define otbBandMathImageFilter_hxx
24 #include "otbBandMathImageFilter.h"
25 
26 #include "itkImageRegionIterator.h"
27 #include "itkNumericTraits.h"
28 #include "itkProgressReporter.h"
29 #include "otbMacro.h"
30 
31 
32 #include <iostream>
33 #include <string>
34 
35 namespace otb
36 {
37 
39 template <class TImage>
42 {
43  //This number will be incremented each time an image
44  //is added over the one minimumrequired
45  this->SetNumberOfRequiredInputs( 1 );
46  this->InPlaceOff();
48 
49  m_UnderflowCount = 0;
50  m_OverflowCount = 0;
51  m_ThreadUnderflow.SetSize(1);
52  m_ThreadOverflow.SetSize(1);
53 }
54 
56 template <class TImage>
59 {
60 }
61 
62 template <class TImage>
64 ::PrintSelf(std::ostream& os, itk::Indent indent) const
65 {
66  Superclass::PrintSelf(os, indent);
67 
68  os << indent << "Expression: " << m_Expression << std::endl;
69  os << indent << "Computed values follow:" << std::endl;
70  os << indent << "UnderflowCount: " << m_UnderflowCount << std::endl;
71  os << indent << "OverflowCount: " << m_OverflowCount << std::endl;
72  os << indent << "itk::NumericTraits<PixelType>::NonpositiveMin() : "
74  os << indent << "itk::NumericTraits<PixelType>::max() : "
75  << itk::NumericTraits<PixelType>::max() << std::endl;
76 }
77 
78 template <class TImage>
81 {
82  this->SetInput(idx, const_cast<TImage *>( image ));
83  DataObjectPointerArraySizeType nbInput = this->GetNumberOfInputs();
84  m_VVarName.resize(nbInput+4);
85  std::ostringstream varName;
86  varName << "b" << nbInput;
87  m_VVarName[idx] = varName.str();
88 
89  m_VVarName[idx+1] = "idxX";
90  m_VVarName[idx+2] = "idxY";
91  m_VVarName[idx+3] = "idxPhyX";
92  m_VVarName[idx+4] = "idxPhyY";
93 }
94 
95 template <class TImage>
97 ::SetNthInput(DataObjectPointerArraySizeType idx, const ImageType * image, const std::string& varName)
98 {
99  this->SetInput(idx, const_cast<TImage *>( image ));
100  m_VVarName.resize(this->GetNumberOfInputs()+4);
101  m_VVarName[idx] = varName;
102 
103  m_VVarName[idx+1] = "idxX";
104  m_VVarName[idx+2] = "idxY";
105  m_VVarName[idx+3] = "idxPhyX";
106  m_VVarName[idx+4] = "idxPhyY";
107 }
108 
109 template <class TImage>
111 ::SetNthInputName(DataObjectPointerArraySizeType idx, const std::string& varName)
112 {
113  m_VVarName[idx] = varName;
114 }
115 
116 template <typename TImage>
119 {
120  return const_cast<TImage *>(this->GetInput(idx));
121 }
122 
123 template< typename TImage >
125 ::SetExpression(const std::string& expression)
126 {
127  if (m_Expression != expression)
128  m_Expression = expression;
129  this->Modified();
130 }
131 
132 template< typename TImage >
133 std::string BandMathImageFilter<TImage>
135 {
136  return m_Expression;
137 }
138 
139 template< typename TImage >
140 std::string BandMathImageFilter<TImage>
142 {
143  return m_VVarName.at(idx);
144 }
145 
146 template< typename TImage >
149 {
150  typename std::vector<ParserType::Pointer>::iterator itParser;
151  unsigned int nbThreads = this->GetNumberOfThreads();
152  unsigned int nbInputImages = this->GetNumberOfInputs();
153  unsigned int nbAccessIndex = 4; //to give access to image and physical index
154  unsigned int i, j;
155  unsigned int inputSize[2];
156  std::vector< std::string > tmpIdxVarNames;
157 
158  tmpIdxVarNames.resize(nbAccessIndex);
159 
160  tmpIdxVarNames.resize(nbAccessIndex);
161  tmpIdxVarNames[0] = "idxX";
162  tmpIdxVarNames[1] = "idxY";
163  tmpIdxVarNames[2] = "idxPhyX";
164  tmpIdxVarNames[3] = "idxPhyY";
165 
166  // Check if input image dimensions matches
167  inputSize[0] = this->GetNthInput(0)->GetLargestPossibleRegion().GetSize(0);
168  inputSize[1] = this->GetNthInput(0)->GetLargestPossibleRegion().GetSize(1);
169 
170  for(unsigned int p = 1; p < nbInputImages; p++)
171  {
172  if((inputSize[0] != this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(0))
173  || (inputSize[1] != this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(1)))
174  {
175  itkExceptionMacro(<< "Input images must have the same dimensions." << std::endl
176  << "band #1 is [" << inputSize[0] << ";" << inputSize[1] << "]" << std::endl
177  << "band #" << p+1 << " is ["
178  << this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(0) << ";"
179  << this->GetNthInput(p)->GetLargestPossibleRegion().GetSize(1) << "]");
180  }
181  }
182 
183  // Store images specs
184  m_Spacing = this->GetNthInput(0)->GetSignedSpacing();
185  m_Origin = this->GetNthInput(0)->GetOrigin();
186 
187  // Allocate and initialize the thread temporaries
188  m_ThreadUnderflow.SetSize(nbThreads);
189  m_ThreadUnderflow.Fill(0);
190  m_ThreadOverflow.SetSize(nbThreads);
191  m_ThreadOverflow.Fill(0);
192  m_VParser.resize(nbThreads);
193  m_AImage.resize(nbThreads);
194  m_NbVar = nbInputImages+nbAccessIndex;
195  m_VVarName.resize(m_NbVar);
196 
197  for(itParser = m_VParser.begin(); itParser < m_VParser.end(); itParser++)
198  {
199  *itParser = ParserType::New();
200  }
201 
202  for(i = 0; i < nbThreads; ++i)
203  {
204  m_AImage[i].resize(m_NbVar);
205  m_VParser[i]->SetExpr(m_Expression);
206 
207  for(j=0; j < nbInputImages; ++j)
208  {
209  m_VParser[i]->DefineVar(m_VVarName[j], &(m_AImage[i][j]));
210  }
211 
212  for(j=nbInputImages; j < nbInputImages+nbAccessIndex; ++j)
213  {
214  m_VVarName[j] = tmpIdxVarNames[j-nbInputImages];
215  m_VParser[i]->DefineVar(m_VVarName[j], &(m_AImage[i][j]));
216  }
217  }
218 }
219 
220 template< typename TImage >
223 {
224  unsigned int nbThreads = this->GetNumberOfThreads();
225  unsigned int i;
226 
227  m_UnderflowCount = 0;
228  m_OverflowCount = 0;
229 
230  // Accumulate counts for each thread
231  for(i = 0; i < nbThreads; ++i)
232  {
233  m_UnderflowCount += m_ThreadUnderflow[i];
234  m_OverflowCount += m_ThreadOverflow[i];
235  }
236 
237  if((m_UnderflowCount != 0) || (m_OverflowCount!=0))
238  otbWarningMacro(<< std::endl
239  << "The Following Parsed Expression : "
240  << this->GetExpression() << std::endl
241  << "Generated " << m_UnderflowCount << " Underflow(s) "
242  << "And " << m_OverflowCount << " Overflow(s) " << std::endl
243  << "The Parsed Expression, The Inputs And The Output "
244  << "Type May Be Incompatible !");
245 }
246 
247 template< typename TImage >
249 ::ThreadedGenerateData(const ImageRegionType& outputRegionForThread,
250  itk::ThreadIdType threadId)
251 {
252  double value;
253  unsigned int j;
254  unsigned int nbInputImages = this->GetNumberOfInputs();
255 
256  typedef itk::ImageRegionConstIterator<TImage> ImageRegionConstIteratorType;
257 
258  assert(nbInputImages);
259  std::vector< ImageRegionConstIteratorType > Vit(nbInputImages);
260 
261  for(j=0; j < nbInputImages; ++j)
262  {
263  Vit[j] = ImageRegionConstIteratorType (this->GetNthInput(j), outputRegionForThread);
264  }
265 
266  itk::ImageRegionIterator<TImage> ot (this->GetOutput(), outputRegionForThread);
267 
268  // support progress methods/callbacks
269  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
270 
271  std::vector<double> & threadImage = m_AImage[threadId];
272  ParserType::Pointer const& threadParser = m_VParser[threadId];
273  long & threadUnderflow = m_ThreadUnderflow[threadId];
274  long & threadOverflow = m_ThreadOverflow[threadId];
275  ImageRegionConstIteratorType & firstImageRegion = Vit.front(); // alias for better perfs
276  while(!firstImageRegion.IsAtEnd())
277  {
278  for(j=0; j < nbInputImages; ++j)
279  {
280  threadImage[j] = static_cast<double>(Vit[j].Get());
281  }
282 
283  // Image Indexes
284  for(j=0; j < 2; ++j)
285  {
286  threadImage[nbInputImages+j] = static_cast<double>(firstImageRegion.GetIndex()[j]);
287  }
288  for(j=0; j < 2; ++j)
289  {
290  threadImage[nbInputImages+2+j] = static_cast<double>(m_Origin[j])
291  + static_cast<double>(firstImageRegion.GetIndex()[j]) * static_cast<double>(m_Spacing[j]);
292  }
293 
294  try
295  {
296 
297  value = threadParser->Eval();
298  }
299  catch(itk::ExceptionObject& err)
300  {
301  itkExceptionMacro(<< err);
302  }
303 
304  // Case value is equal to -inf or inferior to the minimum value
305  // allowed by the pixelType cast
306  if (value < double(itk::NumericTraits<PixelType>::NonpositiveMin()))
307  {
309  threadUnderflow++;
310  }
311  // Case value is equal to inf or superior to the maximum value
312  // allowed by the pixelType cast
313  else if (value > double(itk::NumericTraits<PixelType>::max()))
314  {
316  threadOverflow++;
317  }
318  else
319  {
320  ot.Set(static_cast<PixelType>(value));
321  }
322 
323  for(j=0; j < nbInputImages; ++j)
324  {
325  ++Vit[j];
326  }
327 
328  ++ot;
329 
330  progress.CompletedPixel();
331 
332  }
333 }
334 
335 }// end namespace otb
336 
337 #endif
void SetNthInput(DataObjectPointerArraySizeType idx, const ImageType *image)
Superclass::DataObjectPointerArraySizeType DataObjectPointerArraySizeType
void SetExpression(const std::string &expression)
void Set(const PixelType &value) const
void PrintSelf(std::ostream &os, itk::Indent indent) const override
std::string GetExpression() const
static ITK_CONSTEXPR_FUNC T max(const T &)
void ThreadedGenerateData(const ImageRegionType &outputRegionForThread, itk::ThreadIdType threadId) override
std::string GetNthInputName(DataObjectPointerArraySizeType idx) const
#define otbWarningMacro(x)
Definition: otbMacro.h:67
void SetNthInputName(DataObjectPointerArraySizeType idx, const std::string &expression)
unsigned int ThreadIdType
static ITK_CONSTEXPR_FUNC T NonpositiveMin()
ImageType::RegionType ImageRegionType
ImageType * GetNthInput(DataObjectPointerArraySizeType idx)