OTB  7.2.0
Orfeo Toolbox
otbLineSegmentDetector.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2020 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 otbLineSegmentDetector_hxx
22 #define otbLineSegmentDetector_hxx
23 
24 #include "otbLineSegmentDetector.h"
25 #include "itkImageRegionIterator.h"
26 #include "itkMinimumMaximumImageCalculator.h"
27 #include "itkImageConstIterator.h"
28 #include "itkNeighborhoodIterator.h"
29 #include "otbPolygon.h"
30 #include "itkCastImageFilter.h"
31 
32 #include "otbRectangle.h"
33 #include "otbRemoteSensingRegion.h"
34 
35 #include "otbMath.h"
36 
37 #include "itkMatrix.h"
38 #include "itkSymmetricEigenAnalysis.h"
39 
40 
41 extern "C" double dlngam_(double* x);
42 extern "C" double dbetai_(double* x, double* a, double* b);
43 
44 namespace otb
45 {
46 
47 template <class TInputImage, class TPrecision>
49 {
50  this->SetNumberOfRequiredInputs(1);
51  this->SetNumberOfRequiredOutputs(1);
52 
53  m_DirectionsAllowed = 1. / 8.;
54  m_Prec = CONST_PI * m_DirectionsAllowed;
55  m_Threshold = 5.2;
56 
58  m_GradientFilter = GradientFilterType::New();
59  m_MagnitudeFilter = MagnitudeFilterType::New();
60  m_OrientationFilter = OrientationFilterType::New();
62 
64  m_UsedPointImage = LabelImageType::New();
65 }
66 
67 template <class TInputImage, class TPrecision>
68 void LineSegmentDetector<TInputImage, TPrecision>::SetInput(const InputImageType* input)
69 {
70  this->Superclass::SetNthInput(0, const_cast<InputImageType*>(input));
71 }
72 
73 template <class TInputImage, class TPrecision>
75 {
76  if (this->GetNumberOfInputs() < 1)
77  {
78  return nullptr;
79  }
80 
81  return static_cast<const InputImageType*>(this->Superclass::GetInput(0));
82 }
83 
84 template <class TInputImage, class TPrecision>
86 {
87  // call the superclass' implementation of this method
88  Superclass::GenerateInputRequestedRegion();
89 
90  // get pointers to the inputs
91  typename InputImageType::Pointer input = const_cast<InputImageType*>(this->GetInput());
92 
93  if (!input)
94  {
95  return;
96  }
97 
98  // The input is necessarily the largest possible region.
99  // For a streamed implementation, use the StreamingLineSegmentDetector filter
100  input->SetRequestedRegionToLargestPossibleRegion();
101 }
102 
103 template <class TInputImage, class TPrecision>
105 {
106  if (this->GetInput()->GetRequestedRegion() != this->GetInput()->GetLargestPossibleRegion())
107  {
108  itkExceptionMacro(<< "Not streamed filter. ERROR : requested region is not the largest possible region.");
109  }
110 
112  m_UsedPointImage->SetRegions(this->GetInput()->GetLargestPossibleRegion());
113  m_UsedPointImage->Allocate();
114  m_UsedPointImage->FillBuffer(0);
116 
118  typedef itk::CastImageFilter<InputImageType, OutputImageType> castFilerType;
119  typename castFilerType::Pointer castFilter = castFilerType::New();
120  castFilter->SetInput(this->GetInput());
122 
124  m_GradientFilter->SetInput(castFilter->GetOutput());
125  m_GradientFilter->SetSigma(0.6);
126  m_MagnitudeFilter->SetInput(m_GradientFilter->GetOutput());
127  m_OrientationFilter->SetInput(m_GradientFilter->GetOutput());
129 
130  m_MagnitudeFilter->Update();
131  m_OrientationFilter->Update();
132 
134  CoordinateHistogramType CoordinateHistogram;
135  CoordinateHistogram = this->SortImageByModulusValue(m_MagnitudeFilter->GetOutput());
136 
138  this->LineSegmentDetection(CoordinateHistogram);
139 
141  this->ComputeRectangles();
142 }
143 
144 template <class TInputImage, class TPrecision>
146 LineSegmentDetector<TInputImage, TPrecision>::SortImageByModulusValue(MagnitudeImagePointerType modulusImage)
147 {
148  RegionType largestRegion = this->GetInput()->GetLargestPossibleRegion();
149 
150  // Compute the minimum region size
151  double logNT = 5. * std::log10(static_cast<double>(largestRegion.GetNumberOfPixels())) / 2.;
152  double log1_p = std::log10(m_DirectionsAllowed);
153  double rapport = logNT / log1_p;
154  m_MinimumRegionSize = static_cast<unsigned int>(-rapport);
155 
156  // Computing the min & max of the image
157  typedef itk::MinimumMaximumImageCalculator<OutputImageType> MinMaxCalculatorFilter;
158  typename MinMaxCalculatorFilter::Pointer minmaxCalculator = MinMaxCalculatorFilter::New();
159 
160  minmaxCalculator->SetImage(modulusImage);
161  minmaxCalculator->ComputeMinimum();
162  OutputPixelType min = minmaxCalculator->GetMinimum();
163  minmaxCalculator->ComputeMaximum();
164  OutputPixelType max = minmaxCalculator->GetMaximum();
165 
167  m_Threshold = m_Threshold * ((max - min) / 255.); // threshold normalized with min & max of the values
168 
170  unsigned int NbBin = 1024;
171  double lengthBin = static_cast<double>((max - min)) / static_cast<double>(NbBin - 1);
172 
175  // New region : without boundaries
176  RegionType region;
177  SizeType size = modulusImage->GetRequestedRegion().GetSize();
178  InputIndexType id = modulusImage->GetRequestedRegion().GetIndex();
179 
180  // Don't take in carre the boundary of the image.
181  // Special cases for streamed call
182  if (modulusImage->GetRequestedRegion().GetIndex()[0] == 0)
183  {
184  id[0]++;
185  size[0]--;
186  if (modulusImage->GetRequestedRegion().GetSize()[0] + modulusImage->GetRequestedRegion().GetIndex()[0] == largestRegion.GetSize(0))
187  size[0]--;
188  }
189  else if (modulusImage->GetRequestedRegion().GetSize()[0] + modulusImage->GetRequestedRegion().GetIndex()[0] == largestRegion.GetSize(0))
190  {
191  size[0]--;
192  }
193 
194  if (modulusImage->GetRequestedRegion().GetIndex()[1] == 0)
195  {
196  id[1]++;
197  size[1]--;
198  if (modulusImage->GetRequestedRegion().GetSize()[1] + modulusImage->GetRequestedRegion().GetIndex()[1] == largestRegion.GetSize(1))
199  size[1]--;
200  }
201  else if (modulusImage->GetRequestedRegion().GetSize()[1] + modulusImage->GetRequestedRegion().GetIndex()[1] == largestRegion.GetSize(1))
202  {
203  size[1]--;
204  }
205 
206  region.SetIndex(id);
207  region.SetSize(size);
208 
209  itk::ImageRegionIterator<OutputImageType> it(modulusImage, region);
210 
211  it.GoToBegin();
212  while (!it.IsAtEnd())
213  {
214  unsigned int bin = static_cast<unsigned int>(static_cast<double>(it.Value() - min) / lengthBin);
215 
216  // Highlights bug 498
217  assert(bin < NbBin);
218 
219  if (it.Value() - m_Threshold > 1e-10)
220  tempHisto[NbBin - bin - 1].push_back(it.GetIndex());
221  else
222  SetPixelToUsed(it.GetIndex());
223 
224  ++it;
225  }
226 
227  return tempHisto;
228 }
229 
230 /**************************************************************************************************************/
231 
235 template <class TInputImage, class TPrecision>
236 void LineSegmentDetector<TInputImage, TPrecision>::LineSegmentDetection(CoordinateHistogramType& CoordinateHistogram)
237 {
238 
240  CoordinateHistogramIteratorType ItCoordinateList = CoordinateHistogram.begin();
241 
242  while (ItCoordinateList != CoordinateHistogram.end())
243  {
244  typename IndexVectorType::iterator ItIndexVector = (*ItCoordinateList).begin();
245  while (ItIndexVector != (*ItCoordinateList).end())
246  {
247  InputIndexType index = *(ItIndexVector);
248 
250  if (this->IsNotUsed(index))
251  {
252  IndexVectorType region;
253  double regionAngle = 0.;
254 
255  bool fail = GrowRegion(index, region, regionAngle);
256 
257  if (!fail)
258  {
259  // region -> rectangle
260  RectangleType rectangle = Region2Rectangle(region, regionAngle);
261 
262  // ImprovRectangle(&rectangle)
263  double nfa = ImproveRectangle(rectangle);
264 
265  double density =
266  (double)region.size() /
267  (std::sqrt((rectangle[2] - rectangle[0]) * (rectangle[2] - rectangle[0]) + (rectangle[3] - rectangle[1]) * (rectangle[3] - rectangle[1])) *
268  rectangle[4]);
269  if (density < 0.7)
270  {
271  nfa = -1;
272  // std::cout << "Density = " << density << std::endl;
273  }
274 
275  // if NFA(ImprovRect(rec)) > 0.
276  if (nfa > 0.)
277  {
278  m_RectangleList.push_back(rectangle);
279  }
280  else
281  {
282  SetRegionToNotIni(region);
283  }
284  }
285  }
286  ++ItIndexVector;
287  }
288  ++ItCoordinateList;
289  }
290 }
291 
292 /**************************************************************************************************************/
293 
298 template <class TInputImage, class TPrecision>
300 {
301  // Output
302  this->GetOutput(0)->SetMetaDataDictionary(this->GetInput()->GetMetaDataDictionary());
303  // Retrieving root node
304  typename DataNodeType::Pointer root = this->GetOutput(0)->GetDataTree()->GetRoot()->Get();
305  // Create the document node
306  typename DataNodeType::Pointer document = DataNodeType::New();
307  document->SetNodeType(otb::DOCUMENT);
308  // Adding the layer to the data tree
309  this->GetOutput(0)->GetDataTree()->Add(document, root);
310  // Create the folder node
311  typename DataNodeType::Pointer folder = DataNodeType::New();
312  folder->SetNodeType(otb::FOLDER);
313  // Adding the layer to the data tree
314  this->GetOutput(0)->GetDataTree()->Add(folder, document);
315  this->GetOutput(0)->SetProjectionRef(this->GetInput()->GetProjectionRef());
317 
318  SpacingType spacing = this->GetInput()->GetSignedSpacing();
319  OriginType origin = this->GetInput()->GetOrigin();
320 
322  RectangleListTypeIterator itRec = m_RectangleList.begin();
323  while (itRec != m_RectangleList.end())
324  {
325  VertexType start, end;
327 
328  start[0] = origin[0] + static_cast<TPrecision>((*itRec)[0]) * spacing[0];
329  start[1] = origin[1] + static_cast<TPrecision>((*itRec)[1]) * spacing[1];
330 
331  end[0] = origin[0] + static_cast<TPrecision>((*itRec)[2]) * spacing[0];
332  end[1] = origin[1] + static_cast<TPrecision>((*itRec)[3]) * spacing[1];
333 
334  typename DataNodeType::Pointer CurrentGeometry = DataNodeType::New();
335  CurrentGeometry->SetNodeId("FEATURE_LINE");
336  CurrentGeometry->SetNodeType(otb::FEATURE_LINE);
337  typename LineType::Pointer line = LineType::New();
338  CurrentGeometry->SetLine(line);
339  this->GetOutput(0)->GetDataTree()->Add(CurrentGeometry, folder);
340  CurrentGeometry->GetLine()->AddVertex(start);
341  CurrentGeometry->GetLine()->AddVertex(end);
342 
343  ++itRec;
344  }
345 
346  return EXIT_SUCCESS;
347 }
348 
349 /**************************************************************************************************************/
350 
354 template <class TInputImage, class TPrecision>
355 void LineSegmentDetector<TInputImage, TPrecision>::CopyRectangle(RectangleType& rDst, RectangleType& rSrc) const
356 {
357  RectangleIteratorType itSrc = rSrc.begin();
358 
359  while (itSrc != rSrc.end())
360  {
361  rDst.push_back(*(itSrc));
362  ++itSrc;
363  }
364 }
365 
366 /**************************************************************************************************************/
367 
372 template <class TInputImage, class TPrecision>
374 {
375  int n = 0;
376  double nfa_new;
377  double delta = 0.5;
378  double delta_2 = delta / 2.0;
379  RectangleType r;
380 
381  double nfa_rect = this->ComputeRectNFA(rec);
382 
383  if (nfa_rect > 0.)
384  return nfa_rect;
385 
386  /*Try to improve the precision of the oriented */
387  CopyRectangle(r, rec);
388  for (n = 0; n < 5; ++n)
389  {
390  r[7] /= 2.0;
391  r[6] = CONST_PI * r[7]; // prec = rec[6]
392  nfa_new = this->ComputeRectNFA(r);
393  if (nfa_new > nfa_rect)
394  {
395  nfa_rect = nfa_new;
396  CopyRectangle(rec, r);
397  }
398  }
399 
400  if (nfa_rect > 0.)
401  return nfa_rect;
402 
403  /*Try to improve the width of the rectangle*/
404  CopyRectangle(r, rec);
405  for (n = 0; n < 5; ++n)
406  {
407  r[4] -= delta; // r[4] is stored as the width
408  nfa_new = this->ComputeRectNFA(r);
409  if (nfa_new > nfa_rect)
410  {
411  nfa_rect = nfa_new;
412  CopyRectangle(rec, r);
413  }
414  }
415  if (nfa_rect > 0.)
416  return nfa_rect;
417 
418  /*Try to improve the extremity of the segments*/
419  CopyRectangle(r, rec);
420  for (n = 0; n < 5; ++n)
421  {
422  if ((r[4] - delta) >= 0.5)
423  {
424  r[0] += -std::sin(r[5]) * delta_2;
425  r[1] += std::cos(r[5]) * delta_2;
426  r[2] += -std::sin(r[5]) * delta_2;
427  r[3] += std::cos(r[5]) * delta_2;
428  r[4] -= delta;
429 
430  nfa_new = this->ComputeRectNFA(r);
431  if (nfa_new > nfa_rect)
432  {
433  nfa_rect = nfa_new;
434  CopyRectangle(rec, r);
435  }
436  }
437  }
438  if (nfa_rect > 0.)
439  return nfa_rect;
440 
441  CopyRectangle(r, rec);
442  for (n = 0; n < 5; ++n)
443  {
444  if ((r[4] - delta) >= 0.5)
445  {
446  r[0] -= -std::sin(r[5]) * delta_2;
447  r[1] -= std::cos(r[5]) * delta_2;
448  r[2] -= -std::sin(r[5]) * delta_2;
449  r[3] -= std::cos(r[5]) * delta_2;
450  r[4] -= delta;
451 
452  nfa_new = this->ComputeRectNFA(r);
453  if (nfa_new > nfa_rect)
454  {
455  nfa_rect = nfa_new;
456  CopyRectangle(rec, r);
457  }
458  }
459  }
460  if (nfa_rect > 0.)
461  return nfa_rect;
462 
463  /*Try to improve the precision again */
464  CopyRectangle(r, rec);
465  for (n = 0; n < 5; ++n)
466  {
467  r[7] /= 2.0;
468  r[6] = CONST_PI * r[7]; // prec = rec[]
469  nfa_new = this->ComputeRectNFA(r);
470  if (nfa_new > nfa_rect)
471  {
472  nfa_rect = nfa_new;
473  CopyRectangle(rec, r);
474  }
475  }
476 
477  return nfa_rect;
478 }
479 
480 /**************************************************************************************************************/
481 
485 template <class TInputImage, class TPrecision>
486 bool LineSegmentDetector<TInputImage, TPrecision>::IsNotUsed(InputIndexType& index) const
487 {
488  bool isNotUsed = false;
489 
490  typedef itk::ImageConstIterator<LabelImageType> ImageIteratorType;
491  RegionType region = m_UsedPointImage->GetLargestPossibleRegion();
492  InputIndexType indexRef = region.GetIndex();
493  ImageIteratorType itLabel(m_UsedPointImage, region);
494  itLabel.GoToBegin();
495 
496  if (m_UsedPointImage->GetLargestPossibleRegion().IsInside(index))
497  {
498  itLabel.SetIndex(index);
499  if (itLabel.Get() == 0)
500  isNotUsed = true;
501  }
502  else
503  {
504  itkExceptionMacro(<< "Can't access to index " << index << ", outside the image largest region (" << indexRef << ", " << region.GetSize() << ")");
505  }
506 
507  return isNotUsed;
508 }
509 /**************************************************************************************************************/
510 
514 template <class TInputImage, class TPrecision>
515 bool LineSegmentDetector<TInputImage, TPrecision>::IsUsed(InputIndexType& index) const
516 {
517  bool isUsed = false;
518 
519  typedef itk::ImageConstIterator<LabelImageType> ImageIteratorType;
520  RegionType region = m_UsedPointImage->GetLargestPossibleRegion();
521  InputIndexType indexRef = region.GetIndex();
522  ImageIteratorType itLabel(m_UsedPointImage, region);
523  itLabel.GoToBegin();
524 
525  if (m_UsedPointImage->GetLargestPossibleRegion().IsInside(index))
526  {
527  itLabel.SetIndex(index);
528  if (itLabel.Get() == 255)
529  isUsed = true;
530  }
531  else
532  {
533  itkExceptionMacro(<< "Can't access to index " << index << ", outside the image largest region (" << indexRef << ", " << region.GetSize() << ")");
534  }
535 
536  return isUsed;
537 }
538 
539 /**************************************************************************************************************/
540 
544 template <class TInputImage, class TPrecision>
545 bool LineSegmentDetector<TInputImage, TPrecision>::IsNotIni(InputIndexType& index) const
546 {
547  bool isNotIni = false;
548 
549  typedef itk::ImageConstIterator<LabelImageType> ImageIteratorType;
550  RegionType region = m_UsedPointImage->GetLargestPossibleRegion();
551  InputIndexType indexRef = region.GetIndex();
552  ImageIteratorType itLabel(m_UsedPointImage, region);
553  itLabel.GoToBegin();
554 
555  if (m_UsedPointImage->GetLargestPossibleRegion().IsInside(index))
556  {
557  itLabel.SetIndex(index);
558  if (itLabel.Get() == 127)
559  isNotIni = true;
560  }
561  else
562  {
563  itkExceptionMacro(<< "Can't access to index " << index << ", outside the image largest region (" << indexRef << ", " << region.GetSize() << ")");
564  }
565 
566  return isNotIni;
567 }
568 
569 /**************************************************************************************************************/
570 
574 template <class TInputImage, class TPrecision>
576 {
577  typedef itk::NeighborhoodIterator<LabelImageType> NeighborhoodLabelIteratorType;
578  typename NeighborhoodLabelIteratorType::SizeType radiusLabel;
579  radiusLabel.Fill(0);
580  NeighborhoodLabelIteratorType itLabel(radiusLabel, m_UsedPointImage, m_UsedPointImage->GetRequestedRegion());
581  itLabel.SetLocation(index);
582  itLabel.SetCenterPixel(255); // 255 : Set the point status to : Used Point
583 }
585 
586 /**************************************************************************************************************/
587 
591 template <class TInputImage, class TPrecision>
593 {
594  typedef itk::NeighborhoodIterator<LabelImageType> NeighborhoodLabelIteratorType;
595  typename NeighborhoodLabelIteratorType::SizeType radiusLabel;
596  radiusLabel.Fill(0);
597  NeighborhoodLabelIteratorType itLabel(radiusLabel, m_UsedPointImage, m_UsedPointImage->GetRequestedRegion());
598  itLabel.SetLocation(index);
599  itLabel.SetCenterPixel(127); // 127 : Set the point status to : Not Ini Point
600 }
602 
603 /**************************************************************************************************************/
604 
608 template <class TInputImage, class TPrecision>
610 {
611  IndexVectorIteratorType it = region.begin();
612  while (it != region.end())
613  {
614  this->SetPixelToNotIni(*it);
615  it++;
616  }
617 }
619 
620 /**************************************************************************************************************/
621 
626 template <class TInputImage, class TPrecision>
627 bool LineSegmentDetector<TInputImage, TPrecision>::GrowRegion(InputIndexType index, IndexVectorType& region, double& regionAngle)
628 {
629 
631  this->SetPixelToUsed(index);
632 
634  typedef itk::ConstNeighborhoodIterator<OutputImageType> NeighborhoodIteratorType;
635  typename NeighborhoodIteratorType::SizeType radius;
636  radius.Fill(1);
637  NeighborhoodIteratorType itNeigh(radius, m_MagnitudeFilter->GetOutput(), m_MagnitudeFilter->GetOutput()->GetRequestedRegion());
638  NeighborhoodIteratorType itNeighDir(radius, m_OrientationFilter->GetOutput(), m_OrientationFilter->GetOutput()->GetRequestedRegion());
640 
642  unsigned int neighSize = itNeigh.GetSize()[0] * itNeigh.GetSize()[1];
643 
645  region.push_back(index);
646  double sumX = 0.;
647  double sumY = 0.;
648 
652  for (unsigned int cpt = 0; cpt < region.size(); ++cpt)
653  {
654  itNeigh.SetLocation(region[cpt]);
655  itNeighDir.SetLocation(region[cpt]);
656  sumX += std::cos(*(itNeighDir.GetCenterValue()));
657  sumY += std::sin(*(itNeighDir.GetCenterValue()));
658  regionAngle = std::atan2(sumY, sumX);
660 
661  unsigned int s = 0;
662  while (s < neighSize)
663  {
664  InputIndexType NeighIndex = itNeigh.GetIndex(s);
665  double angleComp = itNeighDir.GetPixel(s);
666 
668  {
669  if ((this->IsNotUsed(NeighIndex) || this->IsNotIni(NeighIndex)) && this->IsAligned(angleComp, regionAngle, m_Prec))
670  {
671  this->SetPixelToUsed(NeighIndex);
672  region.push_back(NeighIndex);
673  }
674  }
675  ++s;
676  }
677 
680  unsigned int nbPixels = this->GetInput()->GetLargestPossibleRegion().GetNumberOfPixels();
681  if (region.size() > m_MinimumRegionSize && region.size() < nbPixels / 4)
682  {
683  return EXIT_SUCCESS;
684  }
685  else
686  {
687  return EXIT_FAILURE;
688  }
689 }
690 
691 /**************************************************************************************************************/
692 
696 template <class TInputImage, class TPrecision>
697 bool LineSegmentDetector<TInputImage, TPrecision>::IsAligned(double Angle, double regionAngle, double prec) const
698 {
699  double diff = Angle - regionAngle;
700 
701  if (diff < 0.0)
702  diff = -diff;
703  if (diff > 1.5 * CONST_PI)
704  {
705  diff -= CONST_2PI;
706  if (diff < 0.0)
707  diff = -diff;
708  }
709 
710  return diff < prec;
711 }
712 
713 /**************************************************************************************************************/
714 
718 template <class TInputImage, class TPrecision>
720  double regionAngle)
721 {
722 
724  double weight = 0., sumWeight = 0.;
725  double x = 0., y = 0.;
726  double l_min = 0., l_max = 0., l = 0., w = 0., w_min = 0., w_max = 0.;
727 
729  typedef itk::ConstNeighborhoodIterator<OutputImageType> NeighborhoodIteratorType;
730  typename NeighborhoodIteratorType::SizeType radius;
731  radius.Fill(0);
732  NeighborhoodIteratorType itNeigh(radius, m_MagnitudeFilter->GetOutput(), m_MagnitudeFilter->GetOutput()->GetRequestedRegion());
734 
736  IndexVectorIteratorType it = region.begin();
737  while (it != region.end())
738  {
739  itNeigh.SetLocation(*it);
740  weight = *itNeigh.GetCenterValue();
741  x += static_cast<double>((*it)[0]) * weight;
742  y += static_cast<double>((*it)[1]) * weight;
743  sumWeight += weight;
744  ++it;
745  }
746  if (sumWeight < 1e-10)
747  {
748  x = 0.;
749  y = 0.;
750  }
751  else
752  {
753  x /= sumWeight;
754  y /= sumWeight;
755  }
757 
759  double theta = this->ComputeRegionOrientation(region, x, y, regionAngle);
760 
761  /* Length & Width of the rectangle **/
762  typedef std::vector<MagnitudePixelType> MagnitudeVector;
763 
764  RegionType largestRegion = this->GetInput()->GetLargestPossibleRegion();
765  unsigned int Diagonal =
766  static_cast<unsigned int>(vnl_math_hypot(static_cast<double>(largestRegion.GetSize(1)), static_cast<double>(largestRegion.GetSize(0))) + 2);
767 
768  MagnitudeVector sum_l(2 * Diagonal, itk::NumericTraits<MagnitudePixelType>::Zero);
769  MagnitudeVector sum_w(2 * Diagonal, itk::NumericTraits<MagnitudePixelType>::Zero);
770 
771  double dx = std::cos(theta);
772  double dy = std::sin(theta);
773 
774  it = region.begin();
775  while (it != region.end())
776  {
777  itNeigh.SetLocation(*it);
778  weight = *itNeigh.GetCenterValue();
779  l = (static_cast<double>((*it)[0]) - x) * dx + (static_cast<double>((*it)[1]) - y) * dy;
780  w = -(static_cast<double>((*it)[0]) - x) * dy + (static_cast<double>((*it)[1]) - y) * dx;
781 
782  if (l < l_min)
783  l_min = l;
784  if (l > l_max)
785  l_max = l;
786  if (w < w_min)
787  w_min = w;
788  if (w > w_max)
789  w_max = w;
790 
791  sum_l[static_cast<int>(std::floor(l) + 0.5) + Diagonal] += static_cast<MagnitudePixelType>(weight);
792  sum_w[static_cast<int>(std::floor(w) + 0.5) + Diagonal] += static_cast<MagnitudePixelType>(weight);
793 
794  ++it;
795  }
796 
798  double sum_th = 0.01 * sumWeight; /* weight threshold for selecting region */
799  double s = 0.;
800  int i = 0;
801 
802  for (s = 0.0, i = static_cast<int>(l_min); s < sum_th && i <= static_cast<int>(l_max); ++i)
803  s += sum_l[Diagonal + i];
804 
805  double lb = (static_cast<double>(i - 1) - 0.5);
806 
807  for (s = 0.0, i = static_cast<int>(l_max); s < sum_th && i >= static_cast<int>(l_min); i--)
808  s += sum_l[Diagonal + i];
809  double lf = (static_cast<double>(i + 1) + 0.5);
810 
811  for (s = 0.0, i = static_cast<int>(w_min); s < sum_th && i <= static_cast<int>(w_max); ++i)
812  s += sum_w[Diagonal + i];
813 
814  double wr = (static_cast<double>(i - 1) - 0.5);
815 
816  for (s = 0.0, i = static_cast<int>(w_max); s < sum_th && i >= static_cast<int>(w_min); i--)
817  s += sum_w[Diagonal + i];
818 
819  double wl = (static_cast<double>(i + 1) + 0.5);
820 
831  RectangleType rec(8, 0.); // Definition of a
832  // rectangle : 8 components
833 
834  if (std::abs(wl - wr) -
835  std::sqrt(static_cast<double>(largestRegion.GetSize(0) * largestRegion.GetSize(0) + largestRegion.GetSize(1) * largestRegion.GetSize(1))) <
836  1e-10)
837  {
838  rec[0] = (x + lb * dx > 0) ? x + lb * dx : 0.;
839  rec[1] = (y + lb * dy > 0) ? y + lb * dy : 0.;
840  rec[2] = (x + lf * dx > 0) ? x + lf * dx : 0.;
841  rec[3] = (y + lf * dy > 0) ? y + lf * dy : 0;
842  rec[4] = wl - wr;
843  rec[5] = theta;
844  rec[6] = m_Prec;
845  rec[7] = m_DirectionsAllowed;
846 
847  if (rec[4] - 1. < 1e-10)
848  rec[4] = 1.;
849  }
850  return rec;
851 }
852 
853 /**************************************************************************************************************/
854 
859 template <class TInputImage, class TPrecision>
860 double LineSegmentDetector<TInputImage, TPrecision>::ComputeRegionOrientation(IndexVectorType region, double x, double y, double angleRegion) const
861 {
862 
863  double Ixx = 0.0;
864  double Iyy = 0.0;
865  double Ixy = 0.0;
866  double theta = 0.;
867  double weight = 0., sum = 0.;
868 
870  typedef itk::ConstNeighborhoodIterator<OutputImageType> NeighborhoodIteratorType;
871  typename NeighborhoodIteratorType::SizeType radius;
872  radius.Fill(0);
873  NeighborhoodIteratorType itNeigh(radius, m_MagnitudeFilter->GetOutput(), m_MagnitudeFilter->GetOutput()->GetRequestedRegion());
875 
877  IndexVectorIteratorType it = region.begin();
878  while (it != region.end())
879  {
880  itNeigh.SetLocation(*it);
881  weight = *itNeigh.GetCenterValue();
882  double Iyy2 = static_cast<double>((*it)[0]) - x;
883  double Ixx2 = static_cast<double>((*it)[1]) - y;
885 
886  Ixx += Ixx2 * Ixx2 * weight;
887  Iyy += Iyy2 * Iyy2 * weight;
888  Ixy -= Ixx2 * Iyy2 * weight;
889  sum += weight;
890  ++it;
891  }
892 
894  typedef itk::Matrix<double, 2, 2> MatrixType;
895  typedef std::vector<double> MatrixEigenType;
896  MatrixType Inertie, eigenVector;
897  MatrixEigenType eigenMatrix(2, 0.);
898  Inertie[1][1] = Ixx;
899  Inertie[0][0] = Iyy;
900  Inertie[0][1] = Ixy;
901  Inertie[1][0] = Ixy;
902 
903  typedef itk::SymmetricEigenAnalysis<MatrixType, MatrixEigenType> EigenAnalysisType;
904  EigenAnalysisType eigenFilter(2);
905  eigenFilter.ComputeEigenValuesAndVectors(Inertie, eigenMatrix, eigenVector);
906  theta = std::atan2(eigenVector[1][1], -eigenVector[1][0]);
907 
908  /* the previous procedure don't cares orientations,
909  so it could be wrong by 180 degrees.
910  here is corrected if necessary */
911 
912  if (this->angle_diff(theta, angleRegion) > m_Prec)
913  theta += CONST_PI;
914 
915  return theta;
916 }
917 
918 /**************************************************************************************************************/
919 
923 template <class TInputImage, class TPrecision>
924 double LineSegmentDetector<TInputImage, TPrecision>::angle_diff(double a, double b) const
925 {
926  a -= b;
927  while (a <= -CONST_PI)
928  a += CONST_2PI;
929  while (a > CONST_PI)
930  a -= CONST_2PI;
931  if (a < 0.0)
932  a = -a;
933  return a;
934 }
936 
937 /**************************************************************************************************************/
938 
942 template <class TInputImage, class TPrecision>
943 double LineSegmentDetector<TInputImage, TPrecision>::ComputeRectNFA(const RectangleType& rec) const
944 {
945  int NbAligned = 0;
946  double nfa_val = 0.;
947 
954  typedef otb::Rectangle<double> OTBRectangleType;
955  OTBRectangleType::Pointer rectangle = OTBRectangleType::New();
956 
958  for (int i = 0; i < 2; ++i)
959  {
960  typename OTBRectangleType::VertexType vertex;
961  vertex[0] = rec[2 * i];
962  vertex[1] = rec[2 * i + 1];
963  rectangle->AddVertex(vertex);
964  }
965  rectangle->SetWidth(rec[4]);
966  rectangle->SetOrientation(rec[5]);
968 
970  OutputImageDirRegionType region = m_OrientationFilter->GetOutput()->GetLargestPossibleRegion();
971  region.Crop(rectangle->GetBoundingRegion());
973 
974 
975  itk::ImageRegionIterator<OutputImageDirType> it(m_OrientationFilter->GetOutput(), region);
976  it.GoToBegin();
977 
978  int pts = 0;
979 
980  while (!it.IsAtEnd())
981  {
982  if (rectangle->IsInside(it.GetIndex()) && m_OrientationFilter->GetOutput()->GetBufferedRegion().IsInside(it.GetIndex()))
983  {
984  ++pts;
985 
986  if (this->IsAligned(it.Get(), rec[5] /*theta*/, rec[6] /*Prec*/))
987  NbAligned++;
988  }
989  ++it;
990  }
991 
993  RegionType largestRegion = const_cast<Self*>(this)->GetInput()->GetLargestPossibleRegion();
994  double logNT = 5. * std::log10(static_cast<double>(largestRegion.GetNumberOfPixels())) / 2.;
996 
997  nfa_val = NFA(pts, NbAligned, m_DirectionsAllowed, logNT);
998 
999  return nfa_val;
1000 }
1001 
1002 /**************************************************************************************************************/
1003 /*
1004  compute logarithm of binomial NFA
1005  NFA = NT.b(n, k, p)
1006  log10 NFA = log10( NFA )
1007 
1008  n, k, p - binomial parameters.
1009  logNT - logarithm of Number of Tests
1010  */
1011 template <class TInputImage, class TPrecision>
1012 double LineSegmentDetector<TInputImage, TPrecision>::NFA(int n, int k, double p, double logNT) const
1013 {
1014  double val;
1015  double n_d = static_cast<double>(n);
1016  double k_d = static_cast<double>(k);
1017 
1018  if (k == 0)
1019  return -logNT;
1020 
1021  // double x = p;
1022  double a = k_d;
1023  double b = n_d - k_d + 1.0;
1024  val = -logNT - log10(dbetai_(&p, &a, &b));
1025 
1026  if (vnl_math_isinf(val)) /* approximate by the first term of the tail */
1027  {
1028  double x1 = n_d + 1.0;
1029  double x2 = k_d + 1.0;
1030  double x3 = n_d - k_d + 1.0;
1031 
1032  val = -logNT - (dlngam_(&x1) - dlngam_(&x2) - dlngam_(&x3)) / CONST_LN10 - k_d * log10(p) - (n_d - k_d) * log10(1.0 - p);
1033  }
1034  return val;
1035 }
1036 
1040 template <class TInputImage, class TPrecision>
1041 void LineSegmentDetector<TInputImage, TPrecision>::PrintSelf(std::ostream& os, itk::Indent indent) const
1042 {
1043  Superclass::PrintSelf(os, indent);
1044 }
1045 
1046 } // end namespace otb
1047 
1048 #endif
void PrintSelf(std::ostream &os, itk::Indent indent) const override
itk::Size< Monteverdi_DIMENSION > SizeType
Definition: mvdTypes.h:137
virtual double angle_diff(double a, double b) const
virtual void SetRegionToNotIni(IndexVectorType region)
virtual void LineSegmentDetection(CoordinateHistogramType &CoordinateHistogram)
This class represent a Rectangle.
Definition: otbRectangle.h:39
constexpr double CONST_PI
Definition: otbMath.h:49
constexpr double CONST_2PI
Definition: otbMath.h:55
virtual void CopyRectangle(RectangleType &rDst, RectangleType &rSrc) const
virtual double ComputeRegionOrientation(IndexVectorType region, double x, double y, double angleRegion) const
virtual double ImproveRectangle(RectangleType &rectangle) const
constexpr double CONST_LN10
Definition: otbMath.h:48
virtual CoordinateHistogramType SortImageByModulusValue(MagnitudeImagePointerType modulusImage)
virtual bool IsAligned(double Angle, double regionAngle, double prec) const
virtual bool GrowRegion(InputIndexType index, IndexVectorType &region, double &regionAngle)
virtual double NFA(int n, int k, double p, double logNT) const
virtual int ComputeRectangles()
virtual void SetInput(const InputImageType *input)
void GenerateInputRequestedRegion() override
virtual void SetPixelToNotIni(InputIndexType index)
double dlngam_(double *x)
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
virtual bool IsNotIni(InputIndexType &index) const
virtual bool IsNotUsed(InputIndexType &index) const
virtual RectangleType Region2Rectangle(IndexVectorType region, double regionAngle)
std::vector< IndexVectorType > CoordinateHistogramType
void GenerateData() override
double dbetai_(double *x, double *a, double *b)
VectorImageType::SpacingType SpacingType
Definition: mvdTypes.h:175
std::vector< double > RectangleType
virtual double ComputeRectNFA(const RectangleType &rec) const
virtual const InputImageType * GetInput(void)
virtual void SetPixelToUsed(InputIndexType index)
virtual bool IsUsed(InputIndexType &index) const