OTB  6.7.0
Orfeo Toolbox
otbSFSTexturesFunctor.h
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 otbSFSTexturesFunctor_h
22 #define otbSFSTexturesFunctor_h
23 
24 #include "otbMath.h"
25 #include "itkNumericTraits.h"
26 #include <vector>
27 
28 namespace otb
29 {
56 namespace Functor
57 {
58 template<class TIter, class TOutputValue>
60 {
61 public:
63  {
64  m_SpatialThreshold = 100;
67  m_Alpha = 1;
68  this->SetNumberOfDirections(20); // set the step too
69  m_SelectedTextures = std::vector<bool>(6, 1);
70  }
71  virtual ~SFSTexturesFunctor() {}
72 
73  typedef typename TIter::InternalPixelType InternalPixelType;
74  typedef typename TIter::SizeType SizeType;
75  typedef typename TIter::IndexType IndexType;
76  typedef typename TIter::OffsetType OffsetType;
77  typedef TOutputValue OutputValueType;
78  typedef std::vector<OutputValueType> OutputType;
79 
80  void SetSpatialThreshold(unsigned int thresh){ m_SpatialThreshold = thresh; }
83  void SetAlpha(double alpha){ m_Alpha = alpha; }
84  void SetNumberOfDirections(unsigned int D)
85  {
87  m_DirectionStep = CONST_PI / static_cast<double>(D);
88  }
89  void SetDirectionStep(double step){ m_DirectionStep = step; }
90  void SetSelectedTextures(std::vector<bool> vect)
91  {
92  m_SelectedTextures.clear();
93  m_SelectedTextures = vect;
94  }
95  void SetTextureStatus(unsigned int id, bool isSelected){ m_SelectedTextures[id] = isSelected; }
96 
97  // Define some getters using itkGetMacro
98  itkGetMacro(SpatialThreshold,unsigned int);
99  itkGetMacro(SpectralThreshold,InternalPixelType);
100  itkGetMacro(RatioMaxConsiderationNumber,unsigned int);
101  itkGetMacro(Alpha,double);
102  itkGetMacro(NumberOfDirections,unsigned int);
103 
104  std::vector<bool> GetTexturesStatus(){ return m_SelectedTextures; }
105 
106  inline OutputType operator ()(const TIter& it)
107  {
109  double width = itk::NumericTraits<double>::max();
110  double SpatialThresholdDouble = static_cast<double>(m_SpatialThreshold);
111  double NumberOfDirectionsDouble = static_cast<double>(m_NumberOfDirections);
112  double dist = 0.;
113  double angle = 0.;
114  double sdiVal = 0.;
115  double sumWMean = 0.;
116  double sum = 0.;
117  std::vector<double> di(m_NumberOfDirections, 0.);
118  std::vector<double> minSorted(m_RatioMaxConsiderationNumber, width);
119  std::vector<double> maxSorted(m_RatioMaxConsiderationNumber, length);
120  std::vector<double> sti(m_NumberOfDirections, 0.);
121  std::vector<unsigned int> lengthLine(m_NumberOfDirections, 0);
122 
123  std::vector<double>::iterator itVector;
124  OutputType out(6, 0);
125 
126  OffsetType off;
127  off.Fill(0);
128 
129  for (unsigned int d = 0; d < m_NumberOfDirections; d++)
130  {
131  // Current angle direction
132  angle = m_DirectionStep * static_cast<double>(d);
133 
134  // last offset in the direction respecting spatial threshold
135  off[0] = static_cast<int>(std::floor(SpatialThresholdDouble * std::cos(angle) + 0.5));
136  off[1] = static_cast<int>(std::floor(SpatialThresholdDouble * std::sin(angle) + 0.5));
137  // last indices in the direction respecting spectral threshold
138  OffsetType offEnd = this->FindLastOffset(it, off);
139 
140  // Check the opposite side
141  off[0] *= -1.0;
142  off[1] *= -1.0;
143  OffsetType offStart = this->FindLastOffset(it, off);
144 
145  // computes distance = dist between the 2 segment point.
146  dist = std::sqrt(std::pow(static_cast<double>(offEnd[0]-offStart[0]), 2) + std::pow(static_cast<double>(offEnd[1]-offStart[1]), 2));
147 
148  // for length computation
149  if (m_SelectedTextures[0] == true) if (dist > length) length = dist;
150 
151  // for width computation
152  if (m_SelectedTextures[1] == true) if (dist < width) width = dist;
153 
154  // for PSI computation
155  if (m_SelectedTextures[2] == true || m_SelectedTextures[5] == true) sum += dist;
156 
157  // for w-mean computation
158  if (m_SelectedTextures[3] == true) sdiVal = this->ComputeSDi(it, offEnd);
159 
160  // for Ratio computation
161  if (m_SelectedTextures[4] == true)
162  {
163  bool doo = false;
164  itVector = maxSorted.begin();
165  while (itVector != maxSorted.end() && doo == false)
166  {
167  if (dist > (*itVector))
168  {
169  itVector = maxSorted.insert(itVector, dist);
170  maxSorted.pop_back();
171  doo = true;
172  }
173  ++itVector;
174  }
175  doo = false;
176  itVector = minSorted.begin();
177  while (itVector != minSorted.end() && doo == false)
178  {
179  if (dist < (*itVector))
180  {
181  itVector = minSorted.insert(itVector, dist);
182  minSorted.pop_back();
183  doo = true;
184  }
185  ++itVector;
186  }
187  }
188 
189  di[d] = dist;
190  if (m_SelectedTextures[3] == true)
191  {
192  lengthLine[d] = static_cast<unsigned int>(dist); //static_cast<unsigned int>( std::sqrt(std::pow(static_cast<double>(offEnd[0]), 2) + std::pow(static_cast<double>(offEnd[1]), 2)) );
193  sti[d] = sdiVal;
194  if (sdiVal != 0.) sumWMean += (m_Alpha * (dist - 1) * dist /*lengthLine[n]*di[n]*/) / sdiVal;
195  }
196  }
197 
199  // length
200  if (m_SelectedTextures[0] == true) out[0] = static_cast<OutputValueType>(length);
201  // width
202  if (m_SelectedTextures[1] == true) out[1] = static_cast<OutputValueType>(width);
203  // PSI
204  if (m_SelectedTextures[2] == true) out[2] = static_cast<OutputValueType>(sum / NumberOfDirectionsDouble);
205  // w-mean
206  if (m_SelectedTextures[3] == true) out[3] = static_cast<OutputValueType>(sumWMean / NumberOfDirectionsDouble);
207  // ratio
208  if (m_SelectedTextures[4] == true)
209  {
210  double sumMin = 0;
211  double sumMax = 0;
212  for (unsigned int t = 0; t < m_RatioMaxConsiderationNumber; t++)
213  {
214  sumMin += minSorted[t];
215  sumMax += maxSorted[t];
216  }
217  if (sumMax != 0.) out[4] = static_cast<OutputValueType>(std::atan(sumMin / sumMax));
218  else if (sumMax == 0. && sumMin == 0.) out[4] = static_cast<OutputValueType>(1.);
219  }
220  // SD
221  if (m_SelectedTextures[5] == true)
222  {
223  double sumPSI = 0;
224  for (unsigned int n = 0; n < di.size(); ++n)
225  sumPSI += std::pow(di[n] - sumWMean / NumberOfDirectionsDouble, 2);
226  out[5] = static_cast<OutputValueType>(std::sqrt(sumPSI) / (NumberOfDirectionsDouble - 1.));
227  }
228 
229  return out;
230 
231  }
232 
237  OffsetType FindLastOffset(const TIter& it, const OffsetType& stopOffset)
238  {
239  bool res = true;
240  int signX = this->ComputeStep(stopOffset[0]);
241  int signY = this->ComputeStep(stopOffset[1]);
243 
244  OffsetType currentOff;
245  currentOff.Fill(0);
246  currentOff[0] = signX;
247 
248  double slop = 0.;
249  if (stopOffset[0] != 0) slop = static_cast<double>(stopOffset[1] / static_cast<double>(stopOffset[0]));
250 
251  bool isInside = true;
252  while (isInside == true && res == true)
253  {
254  this->ComputePointLine(currentOff, slop, signY, stopOffset[0]);
255 
256  if (std::abs(it.GetPixel(currentOff) - it.GetCenterPixel()) > m_SpectralThreshold)
257  {
258  res = false;
259  }
260  else currentOff[0] += signX;
261 
262  isInside = this->CheckIsInside(signX, signY, currentOff, stopOffset);
263  }
264 
265  return currentOff;
266  }
267 
269  double ComputeSDi(const TIter& it, const OffsetType& stopOffset)
270  {
271  bool canGo = true;
272  unsigned int nbElt = 0;
273  double SDi = 0.;
274  double mean = 0.;
275  double slop = 0.;
276  if (stopOffset[0] != 0) slop = static_cast<double>(stopOffset[1] / static_cast<double>(stopOffset[0]));
278 
279  int signX = this->ComputeStep(stopOffset[0]);
280  int signY = this->ComputeStep(stopOffset[1]);
281 
282  OffsetType currentOff;
283  currentOff.Fill(0);
284  currentOff[0] = signX;
285 
286  bool isInside = true;
287  // First compute mean
288  while (isInside == true && canGo == true)
289  {
290  this->ComputePointLine(currentOff, slop, signY, stopOffset[0]);
291 
292  mean += static_cast<double>(it.GetPixel(currentOff));
293  nbElt++;
294 
295  if (std::abs(it.GetPixel(currentOff) - it.GetCenterPixel()) >= m_SpectralThreshold) canGo = false;
296  else currentOff[0] += signX;
297 
298  isInside = this->CheckIsInside(signX, signY, currentOff, stopOffset);
299  }
300 
301  mean /= static_cast<double>(nbElt);
302  currentOff[0] = signX;
303  currentOff[1] = 0;
304  isInside = true;
305 
306  while (isInside == true && canGo == true)
307  {
308  this->ComputePointLine(currentOff, slop, signY, stopOffset[0]);
309 
310  SDi += std::pow((static_cast<double>(it.GetPixel(currentOff)) - mean), 2);
311  if (std::abs(it.GetPixel(currentOff) - it.GetCenterPixel()) >= m_SpectralThreshold) canGo = false;
312  else currentOff[0] += signX;
313 
314  isInside = this->CheckIsInside(signX, signY, currentOff, stopOffset);
315 
316  }
317  return std::sqrt(SDi);
318  }
319 
321  bool CheckIsInside(const int& signX, const int& signY, const OffsetType& currentOff, const OffsetType& stopOffset)
322  {
323  bool isInside = true;
324  if (signX * currentOff[0] >= signX * stopOffset[0] && stopOffset[0] != 0) isInside = false;
325  else if (signY * currentOff[1] >= signY * stopOffset[1] && stopOffset[1] != 0) isInside = false;
327 
328  return isInside;
329  }
330 
335  void ComputePointLine(OffsetType& currentOff, const double& slop, const int& signY, const int& stopOffsetX)
336  {
337  if (stopOffsetX != 0) currentOff[1] = static_cast<int>(std::floor(slop * static_cast<double>(currentOff[0]) + 0.5));
338  else currentOff[1] += signY;
339  }
341 
345  int ComputeStep(const int& stopOffset)
346  {
347  int sign = 1;
348  if (stopOffset < 0) sign = -1;
350 
351  return sign;
352  }
353 
354 protected:
357 
359  unsigned int m_SpatialThreshold;
360 
363 
365  double m_Alpha;
366 
368  unsigned int m_NumberOfDirections;
369 
372 
382  std::vector<bool> m_SelectedTextures;
383 };
384 
385 } // end namespace functor
386 } // end namespace otb
387 
388 #endif
itk::Size< Monteverdi_DIMENSION > SizeType
Definition: mvdTypes.h:137
bool CheckIsInside(const int &signX, const int &signY, const OffsetType &currentOff, const OffsetType &stopOffset)
constexpr double CONST_PI
Definition: otbMath.h:48
double ComputeSDi(const TIter &it, const OffsetType &stopOffset)
OutputType operator()(const TIter &it)
static ITK_CONSTEXPR_FUNC T max(const T &)
itk::Index< Monteverdi_DIMENSION > IndexType
Definition: mvdTypes.h:133
void SetSelectedTextures(std::vector< bool > vect)
OffsetType FindLastOffset(const TIter &it, const OffsetType &stopOffset)
std::vector< OutputValueType > OutputType
void SetTextureStatus(unsigned int id, bool isSelected)
static ITK_CONSTEXPR_FUNC T NonpositiveMin()
void ComputePointLine(OffsetType &currentOff, const double &slop, const int &signY, const int &stopOffsetX)
void SetSpectralThreshold(InternalPixelType thresh)
TIter::InternalPixelType InternalPixelType
void SetRatioMaxConsiderationNumber(unsigned int value)
void SetNumberOfDirections(unsigned int D)
int ComputeStep(const int &stopOffset)
void SetSpatialThreshold(unsigned int thresh)