Orfeo Toolbox  3.16
otbVectorizationPathListFilter.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12 
13  This software is distributed WITHOUT ANY WARRANTY; without even
14  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15  PURPOSE. See the above copyright notices for more information.
16 
17 =========================================================================*/
18 #ifndef __otbVectorizationPathListFilter_txx
19 #define __otbVectorizationPathListFilter_txx
20 
21 #include "otbVectorizationPathListFilter.h"
22 #include "otbMacro.h"
23 #include "otbMath.h"
24 
25 namespace otb
26 {
30 template <class TInputModulus, class TInputDirection, class TOutputPath>
31 VectorizationPathListFilter<TInputModulus, TInputDirection, TOutputPath>
32 ::VectorizationPathListFilter()
33 {
34  this->SetNumberOfRequiredInputs(2);
35  this->SetNumberOfInputs(2);
36  m_AmplitudeThreshold = 1.0;
37 }
38 
39 template <class TInputModulus, class TInputDirection, class TOutputPath>
40 void
41 VectorizationPathListFilter<TInputModulus, TInputDirection, TOutputPath>
42 ::SetInput(InputModulusType * inputModulus)
43 {
44  this->itk::ProcessObject::SetNthInput(0, const_cast<InputModulusType *>(inputModulus));
45 }
46 
47 template <class TInputModulus, class TInputDirection, class TOutputPath>
48 const typename VectorizationPathListFilter<TInputModulus, TInputDirection, TOutputPath>
49 ::InputModulusType *
50 VectorizationPathListFilter<TInputModulus, TInputDirection, TOutputPath>
51 ::GetInput(void)
52 {
53  if (this->GetNumberOfInputs() < 1)
54  {
55  return 0;
56  }
57  return static_cast<const TInputModulus*>(this->itk::ProcessObject::GetInput(0));
58 }
59 
60 template <class TInputModulus, class TInputDirection, class TOutputPath>
61 void
62 VectorizationPathListFilter<TInputModulus, TInputDirection, TOutputPath>
63 ::SetInputDirection(InputDirectionType * inputDirection)
64 {
65  this->itk::ProcessObject::SetNthInput(1, const_cast<InputDirectionType *>(inputDirection));
66 }
67 
68 template <class TInputModulus, class TInputDirection, class TOutputPath>
69 const typename VectorizationPathListFilter<TInputModulus, TInputDirection, TOutputPath>
70 ::InputDirectionType *
71 VectorizationPathListFilter<TInputModulus, TInputDirection, TOutputPath>
72 ::GetInputDirection(void)
73 {
74  if (this->GetNumberOfInputs() < 2)
75  {
76  return 0;
77  }
78  return static_cast<const TInputDirection *>(this->itk::ProcessObject::GetInput(1));
79 }
80 
84 template <class TInputModulus, class TInputDirection, class TOutputPath>
85 void
86 VectorizationPathListFilter<TInputModulus, TInputDirection, TOutputPath>
87 ::GenerateData(void)
88 {
89  InputModulusConstPointerType modPtr = this->GetInput();
90  InputDirectionConstPointerType dirPtr = this->GetInputDirection();
91  OutputPathListPointerType outPtr = this->GetOutput();
92 
93  typedef typename OffsetVectorType::iterator OffsetIteratorType;
94 
95  RadiusType radius;
96  radius.Fill(2);
97  OffsetVectorType offsetVector;
98 
99  // Creation of the flag image
100  FlagImagePointerType flagImage = FlagImageType::New();
101  flagImage->SetRegions(modPtr->GetLargestPossibleRegion());
102  flagImage->Allocate();
103  flagImage->FillBuffer(false);
104 
105  // Iterators instantiation
106  ModRegionIteratorType modIt(modPtr, modPtr->GetLargestPossibleRegion());
107  DirRegionIteratorType dirIt(dirPtr, dirPtr->GetLargestPossibleRegion());
108  FlagRegionIteratorType flagIt(flagImage, flagImage->GetLargestPossibleRegion());
109 
110  for (modIt.GoToBegin(), dirIt.GoToBegin(), flagIt.GoToBegin();
111  (!modIt.IsAtEnd()) && (!dirIt.IsAtEnd()) && (!flagIt.IsAtEnd());
112  ++modIt, ++dirIt, ++flagIt)
113  {
114  if ((modIt.Get() > m_AmplitudeThreshold) && (!flagIt.Get()))
115  {
116  //this is a beginning, to follow in two directions
117  OutputPathPointerType pathTempDirect = OutputPathType::New();
118  OutputPathPointerType pathTempReverse = OutputPathType::New();
119  OutputPathPointerType path = OutputPathType::New();
120 
121  bool flagFinish;
122  int flagReverse = 0;
123  double totalAmplitude = 0;
124 
125  ModNeighborhoodIteratorType nModIt(radius, modPtr, modPtr->GetLargestPossibleRegion());
126  DirNeighborhoodIteratorType nDirIt(radius, dirPtr, dirPtr->GetLargestPossibleRegion());
127  FlagNeighborhoodIteratorType nFlagIt(radius, flagImage, flagImage->GetLargestPossibleRegion());
128 
129  for (flagReverse = 0; flagReverse < 2; ++flagReverse)
130  {
131  nModIt.SetLocation(modIt.GetIndex());
132  nDirIt.SetLocation(dirIt.GetIndex());
133  nFlagIt.SetLocation(flagIt.GetIndex());
134  // temporary point
135  PointType point;
136  VertexType vertex;
137  modPtr->TransformIndexToPhysicalPoint(nModIt.GetIndex(), point);
138  modPtr->TransformPhysicalPointToContinuousIndex(point, vertex);
139  if (flagReverse == 0)
140  {
141  flagIt.Set(true);
142 
143  // otbMsgDebugMacro(<<"Adding new vertex: "<<vertex);
144 
145  pathTempDirect->AddVertex(vertex);
146  }
147  flagFinish = false;
148  while (!flagFinish)
149  {
150  offsetVector = GetThreeNeighborOffsetFromDirection(nDirIt.GetCenterPixel(), flagReverse);
151  OffsetIteratorType vecIt = offsetVector.begin();
152  bool flagFound = false;
153  while (vecIt != offsetVector.end() && !flagFound)
154  {
155  flagFound = nModIt.GetPixel(*vecIt) > 0
156  && !nFlagIt.GetPixel(*vecIt);
157  ++vecIt;
158  }
159  if (flagFound)
160  {
161  point.Fill(0);
162  PointType tmpPoint;
163  totalAmplitude = 0;
164  for (vecIt = offsetVector.begin(); vecIt != offsetVector.end(); ++vecIt)
165  {
166  double currentAmplitude = nModIt.GetPixel(*vecIt);
167  modPtr->TransformIndexToPhysicalPoint(nModIt.GetIndex(*vecIt), tmpPoint);
168  point[0] += currentAmplitude * tmpPoint[0];
169  point[1] += currentAmplitude * tmpPoint[1];
170  totalAmplitude += currentAmplitude;
171  }
172  point[0] = point[0] / totalAmplitude;
173  point[1] = point[1] / totalAmplitude;
174  modPtr->TransformPhysicalPointToContinuousIndex(point, vertex);
175  if (flagReverse == 0)
176  {
177 // otbMsgDevMacro(<<"Adding new vertex (direct): "<<vertex);
178 
179  pathTempDirect->AddVertex(vertex);
180  }
181  else
182  {
183 
184 // otbMsgDevMacro(<<"Adding new vertex (reverse): "<<vertex);
185 
186  pathTempReverse->AddVertex(vertex);
187  }
188  // flag the pixel use
189  nFlagIt.SetCenterPixel(true);
190  //update the neighbor iterators so they are centered on the nearest pixel to the barycenter
191  IndexType newIndex;
192  if (modPtr->TransformPhysicalPointToIndex(point, newIndex))
193  {
194 // otbMsgDevMacro(<<"Moving to new center: " << newIndex);
195  nModIt.SetLocation(newIndex);
196  nDirIt.SetLocation(newIndex);
197  nFlagIt.SetLocation(newIndex);
198 
199  if (nModIt.GetCenterPixel() == 0)
200  {
201  //we need to check that in case the barycenter is out...
202  flagFinish = true;
203  }
204  if (nFlagIt.GetCenterPixel())
205  {
206  //we don't want to go back to the same pixels
207  flagFinish = true;
208  }
209  }
210  else
211  {
212  //new point outside image
213  flagFinish = true;
214  }
215  }
216  else
217  {
218  flagFinish = true;
219  }
220  }
221  }
222  VertexListPointerType vertexDirect = pathTempDirect->GetVertexList();
223  VertexListPointerType vertexReverse = pathTempReverse->GetVertexList();
224 
225  unsigned int numberVertex = 0;
226 
227  VertexIteratorType vertexReverseIt = vertexReverse->End();
228  if (vertexReverseIt != vertexReverse->Begin())
229  {
230  --vertexReverseIt;
231  while (vertexReverseIt != vertexReverse->Begin())
232  {
233  path->AddVertex(vertexReverseIt.Value());
234  ++numberVertex;
235  --vertexReverseIt;
236  }
237  path->AddVertex(vertexReverseIt.Value());
238  }
239 
240  VertexIteratorType vertexDirectIt = vertexDirect->Begin();
241  while (vertexDirectIt != vertexDirect->End())
242  {
243  path->AddVertex(vertexDirectIt.Value());
244  ++vertexDirectIt;
245  ++numberVertex;
246  }
247 
248  // otbMsgDebugMacro(<<"Path number of vertices: "<<numberVertex);
249 
250  if (numberVertex > 3)
251  {
252  outPtr->PushBack(path);
253  }
254  }
255  }
256 }
257 
264 template <class TInputModulus, class TInputDirection, class TOutputPath>
265 typename VectorizationPathListFilter<TInputModulus, TInputDirection, TOutputPath>
266 ::OffsetVectorType
267 VectorizationPathListFilter<TInputModulus, TInputDirection, TOutputPath>
268 ::GetEightNeighborOffsetFromDirection(double direction, unsigned int flagReverse)
269 {
270  int neighborhoodNumber = 0;
271  OffsetVectorType offset;
272  offset.reserve(8);
273  if (direction > 0)
274  {
275  //find the direction in terms of 0, 1, 2, 3
276  neighborhoodNumber = (int) (direction / (CONST_PI_4) -1);
277  }
278  else
279  {
280  neighborhoodNumber = (int) ((direction + CONST_PI) / (CONST_PI_4) -1);
281  neighborhoodNumber = (neighborhoodNumber + 4);
282  //if the direction was <0 need to convert to 4, 5, 6, 7
283  }
284  if (flagReverse)
285  {
286  //if the reverse flag is activated we need to look on the other side
287  neighborhoodNumber = (neighborhoodNumber + 4) % 8;
288  }
289  OffsetType tmpOffset;
290  switch (neighborhoodNumber)
291  {
292  case 0:
293  tmpOffset[0] = 1;
294  tmpOffset[1] = 0;
295  offset.push_back(tmpOffset);
296  tmpOffset[0] = 1;
297  tmpOffset[1] = 1;
298  offset.push_back(tmpOffset);
299  tmpOffset[0] = 0;
300  tmpOffset[1] = 1;
301  offset.push_back(tmpOffset);
302 
303  tmpOffset[0] = 2;
304  tmpOffset[1] = 0;
305  offset.push_back(tmpOffset);
306  tmpOffset[0] = 2;
307  tmpOffset[1] = 1;
308  offset.push_back(tmpOffset);
309  tmpOffset[0] = 2;
310  tmpOffset[1] = 2;
311  offset.push_back(tmpOffset);
312  tmpOffset[0] = 1;
313  tmpOffset[1] = 2;
314  offset.push_back(tmpOffset);
315  tmpOffset[0] = 0;
316  tmpOffset[1] = 2;
317  offset.push_back(tmpOffset);
318 
319  break;
320 
321  case 1:
322  tmpOffset[0] = 1;
323  tmpOffset[1] = 1;
324  offset.push_back(tmpOffset);
325  tmpOffset[0] = 0;
326  tmpOffset[1] = 1;
327  offset.push_back(tmpOffset);
328  tmpOffset[0] = -1;
329  tmpOffset[1] = 1;
330  offset.push_back(tmpOffset);
331 
332  tmpOffset[0] = 2;
333  tmpOffset[1] = 2;
334  offset.push_back(tmpOffset);
335  tmpOffset[0] = 1;
336  tmpOffset[1] = 2;
337  offset.push_back(tmpOffset);
338  tmpOffset[0] = 0;
339  tmpOffset[1] = 2;
340  offset.push_back(tmpOffset);
341  tmpOffset[0] = -1;
342  tmpOffset[1] = 2;
343  offset.push_back(tmpOffset);
344  tmpOffset[0] = -2;
345  tmpOffset[1] = 2;
346  offset.push_back(tmpOffset);
347  break;
348 
349  case 2:
350  tmpOffset[0] = 0;
351  tmpOffset[1] = 1;
352  offset.push_back(tmpOffset);
353  tmpOffset[0] = -1;
354  tmpOffset[1] = 1;
355  offset.push_back(tmpOffset);
356  tmpOffset[0] = -1;
357  tmpOffset[1] = 0;
358  offset.push_back(tmpOffset);
359 
360  tmpOffset[0] = 0;
361  tmpOffset[1] = 2;
362  offset.push_back(tmpOffset);
363  tmpOffset[0] = -1;
364  tmpOffset[1] = 2;
365  offset.push_back(tmpOffset);
366  tmpOffset[0] = -2;
367  tmpOffset[1] = 2;
368  offset.push_back(tmpOffset);
369  tmpOffset[0] = -2;
370  tmpOffset[1] = 1;
371  offset.push_back(tmpOffset);
372  tmpOffset[0] = -2;
373  tmpOffset[1] = 0;
374  offset.push_back(tmpOffset);
375  break;
376 
377  case 3:
378  tmpOffset[0] = -1;
379  tmpOffset[1] = 1;
380  offset.push_back(tmpOffset);
381  tmpOffset[0] = -1;
382  tmpOffset[1] = 0;
383  offset.push_back(tmpOffset);
384  tmpOffset[0] = -1;
385  tmpOffset[1] = -1;
386  offset.push_back(tmpOffset);
387 
388  tmpOffset[0] = -2;
389  tmpOffset[1] = 2;
390  offset.push_back(tmpOffset);
391  tmpOffset[0] = -2;
392  tmpOffset[1] = 1;
393  offset.push_back(tmpOffset);
394  tmpOffset[0] = -2;
395  tmpOffset[1] = 0;
396  offset.push_back(tmpOffset);
397  tmpOffset[0] = -2;
398  tmpOffset[1] = -1;
399  offset.push_back(tmpOffset);
400  tmpOffset[0] = -2;
401  tmpOffset[1] = -2;
402  offset.push_back(tmpOffset);
403  break;
404 
405  case 4:
406  tmpOffset[0] = -1;
407  tmpOffset[1] = 0;
408  offset.push_back(tmpOffset);
409  tmpOffset[0] = -1;
410  tmpOffset[1] = -1;
411  offset.push_back(tmpOffset);
412  tmpOffset[0] = 0;
413  tmpOffset[1] = -1;
414  offset.push_back(tmpOffset);
415 
416  tmpOffset[0] = -2;
417  tmpOffset[1] = 0;
418  offset.push_back(tmpOffset);
419  tmpOffset[0] = -2;
420  tmpOffset[1] = -1;
421  offset.push_back(tmpOffset);
422  tmpOffset[0] = -2;
423  tmpOffset[1] = -2;
424  offset.push_back(tmpOffset);
425  tmpOffset[0] = -1;
426  tmpOffset[1] = -2;
427  offset.push_back(tmpOffset);
428  tmpOffset[0] = 0;
429  tmpOffset[1] = -2;
430  offset.push_back(tmpOffset);
431  break;
432 
433  case 5:
434  tmpOffset[0] = -1;
435  tmpOffset[1] = -1;
436  offset.push_back(tmpOffset);
437  tmpOffset[0] = 0;
438  tmpOffset[1] = -1;
439  offset.push_back(tmpOffset);
440  tmpOffset[0] = 1;
441  tmpOffset[1] = -1;
442  offset.push_back(tmpOffset);
443 
444  tmpOffset[0] = -2;
445  tmpOffset[1] = -2;
446  offset.push_back(tmpOffset);
447  tmpOffset[0] = -1;
448  tmpOffset[1] = -2;
449  offset.push_back(tmpOffset);
450  tmpOffset[0] = 0;
451  tmpOffset[1] = -2;
452  offset.push_back(tmpOffset);
453  tmpOffset[0] = 1;
454  tmpOffset[1] = -2;
455  offset.push_back(tmpOffset);
456  tmpOffset[0] = 2;
457  tmpOffset[1] = -2;
458  offset.push_back(tmpOffset);
459  break;
460 
461  case 6:
462  tmpOffset[0] = 0;
463  tmpOffset[1] = -1;
464  offset.push_back(tmpOffset);
465  tmpOffset[0] = 1;
466  tmpOffset[1] = -1;
467  offset.push_back(tmpOffset);
468  tmpOffset[0] = 1;
469  tmpOffset[1] = 0;
470  offset.push_back(tmpOffset);
471 
472  tmpOffset[0] = 0;
473  tmpOffset[1] = -2;
474  offset.push_back(tmpOffset);
475  tmpOffset[0] = 1;
476  tmpOffset[1] = -2;
477  offset.push_back(tmpOffset);
478  tmpOffset[0] = 2;
479  tmpOffset[1] = -2;
480  offset.push_back(tmpOffset);
481  tmpOffset[0] = 2;
482  tmpOffset[1] = -1;
483  offset.push_back(tmpOffset);
484  tmpOffset[0] = 2;
485  tmpOffset[1] = 0;
486  offset.push_back(tmpOffset);
487  break;
488 
489  case 7:
490  tmpOffset[0] = 1;
491  tmpOffset[1] = -1;
492  offset.push_back(tmpOffset);
493  tmpOffset[0] = 1;
494  tmpOffset[1] = 0;
495  offset.push_back(tmpOffset);
496  tmpOffset[0] = 1;
497  tmpOffset[1] = 1;
498  offset.push_back(tmpOffset);
499 
500  tmpOffset[0] = 2;
501  tmpOffset[1] = -2;
502  offset.push_back(tmpOffset);
503  tmpOffset[0] = 2;
504  tmpOffset[1] = -1;
505  offset.push_back(tmpOffset);
506  tmpOffset[0] = 2;
507  tmpOffset[1] = 0;
508  offset.push_back(tmpOffset);
509  tmpOffset[0] = 2;
510  tmpOffset[1] = 1;
511  offset.push_back(tmpOffset);
512  tmpOffset[0] = 2;
513  tmpOffset[1] = 2;
514  offset.push_back(tmpOffset);
515  break;
516  }
517  return offset;
518 }
519 
526 template <class TInputModulus, class TInputDirection, class TOutputPath>
527 typename VectorizationPathListFilter<TInputModulus, TInputDirection, TOutputPath>
528 ::OffsetVectorType
529 VectorizationPathListFilter<TInputModulus, TInputDirection, TOutputPath>
530 ::GetThreeNeighborOffsetFromDirection(double direction, unsigned int flagReverse)
531 {
532  int neighborhoodNumber = 0;
533  OffsetVectorType offset;
534  offset.reserve(3);
535  if (direction > 0)
536  {
537  //find the direction in terms of 0, 1, 2, 3
538  neighborhoodNumber = (int) (direction / (CONST_PI_4) -1);
539  }
540  else
541  {
542  neighborhoodNumber = (int) ((direction + CONST_PI) / (CONST_PI_4) -1);
543  neighborhoodNumber = (neighborhoodNumber + 4);
544  //if the direction was <0 need to convert to 4, 5, 6, 7
545  }
546  if (flagReverse)
547  {
548  //if the reverse flag is activated we need to look on the other side
549  neighborhoodNumber = (neighborhoodNumber + 4) % 8;
550  }
551  OffsetType tmpOffset;
552 // otbMsgDevMacro(<<"Direction: " << neighborhoodNumber)
553  switch (neighborhoodNumber)
554  {
555  case 0:
556  tmpOffset[0] = 1;
557  tmpOffset[1] = 0;
558  offset.push_back(tmpOffset);
559  tmpOffset[0] = 1;
560  tmpOffset[1] = 1;
561  offset.push_back(tmpOffset);
562  tmpOffset[0] = 0;
563  tmpOffset[1] = 1;
564  offset.push_back(tmpOffset);
565 
566  break;
567 
568  case 1:
569  tmpOffset[0] = 1;
570  tmpOffset[1] = 1;
571  offset.push_back(tmpOffset);
572  tmpOffset[0] = 0;
573  tmpOffset[1] = 1;
574  offset.push_back(tmpOffset);
575  tmpOffset[0] = -1;
576  tmpOffset[1] = 1;
577  offset.push_back(tmpOffset);
578 
579  break;
580 
581  case 2:
582  tmpOffset[0] = 0;
583  tmpOffset[1] = 1;
584  offset.push_back(tmpOffset);
585  tmpOffset[0] = -1;
586  tmpOffset[1] = 1;
587  offset.push_back(tmpOffset);
588  tmpOffset[0] = -1;
589  tmpOffset[1] = 0;
590  offset.push_back(tmpOffset);
591 
592  break;
593 
594  case 3:
595  tmpOffset[0] = -1;
596  tmpOffset[1] = 1;
597  offset.push_back(tmpOffset);
598  tmpOffset[0] = -1;
599  tmpOffset[1] = 0;
600  offset.push_back(tmpOffset);
601  tmpOffset[0] = -1;
602  tmpOffset[1] = -1;
603  offset.push_back(tmpOffset);
604 
605  break;
606 
607  case 4:
608  tmpOffset[0] = -1;
609  tmpOffset[1] = 0;
610  offset.push_back(tmpOffset);
611  tmpOffset[0] = -1;
612  tmpOffset[1] = -1;
613  offset.push_back(tmpOffset);
614  tmpOffset[0] = 0;
615  tmpOffset[1] = -1;
616  offset.push_back(tmpOffset);
617 
618  break;
619 
620  case 5:
621  tmpOffset[0] = -1;
622  tmpOffset[1] = -1;
623  offset.push_back(tmpOffset);
624  tmpOffset[0] = 0;
625  tmpOffset[1] = -1;
626  offset.push_back(tmpOffset);
627  tmpOffset[0] = 1;
628  tmpOffset[1] = -1;
629  offset.push_back(tmpOffset);
630 
631  break;
632 
633  case 6:
634  tmpOffset[0] = 0;
635  tmpOffset[1] = -1;
636  offset.push_back(tmpOffset);
637  tmpOffset[0] = 1;
638  tmpOffset[1] = -1;
639  offset.push_back(tmpOffset);
640  tmpOffset[0] = 1;
641  tmpOffset[1] = 0;
642  offset.push_back(tmpOffset);
643 
644  break;
645 
646  case 7:
647  tmpOffset[0] = 1;
648  tmpOffset[1] = -1;
649  offset.push_back(tmpOffset);
650  tmpOffset[0] = 1;
651  tmpOffset[1] = 0;
652  offset.push_back(tmpOffset);
653  tmpOffset[0] = 1;
654  tmpOffset[1] = 1;
655  offset.push_back(tmpOffset);
656 
657  break;
658  }
659  return offset;
660 }
661 
665 template <class TInputModulus, class TInputDirection, class TOutputPath>
666 void
667 VectorizationPathListFilter<TInputModulus, TInputDirection, TOutputPath>
668 ::PrintSelf(std::ostream& os, itk::Indent indent) const
669 {
670  Superclass::PrintSelf(os, indent);
671 }
672 
673 } // End namespace otb
674 #endif
Generated at Sun Feb 3 2013 00:54:52 for Orfeo Toolbox with doxygen 1.8.1.1