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