DGtal  1.2.0
Hull2DHelpers.ih
1 /**
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU Lesser General Public License as
4  * published by the Free Software Foundation, either version 3 of the
5  * License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  *
15  **/
16 
17 /**
18  * @file Hull2DHelpers.ih
19  * @author Tristan Roussillon (\c tristan.roussillon@liris.cnrs.fr )
20  * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
21  *
22  * @date 2013/12/02
23  *
24  * Implementation of inline methods defined in Hull2DHelpers.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cstdlib>
32 #include <cmath>
33 
34 #include "DGtal/kernel/PointVector.h"
35 //////////////////////////////////////////////////////////////////////////////
36 
37 ///////////////////////////////////////////////////////////////////////////////
38 // Implementation of inline functions //
39 
40 namespace DGtal
41 {
42  namespace functions
43  {
44  namespace Hull2D
45  {
46 
47  //----------------------------------------------------------------------------
48  template <typename Stack, typename Point, typename Predicate>
49  inline
50  void updateHullWithStack(Stack& aStack, const Point& aNewPoint, const Predicate& aPredicate)
51  {
52  BOOST_CONCEPT_ASSERT(( concepts::CStack<Stack> ));
53 
54  Point Q = aStack.top();
55  aStack.pop();
56  if (aStack.size() != 0)
57  {
58  Point P = aStack.top();
59  while ( ( !aPredicate(P,Q,aNewPoint) )&&(aStack.size() != 0) )
60  {
61  //remove Q
62  Q = P;
63  aStack.pop();
64  if (aStack.size() != 0)
65  P = aStack.top();
66  }
67  //add Q
68  aStack.push(Q);
69  }
70  }
71 
72  //----------------------------------------------------------------------------
73  template <typename Stack, typename Point, typename Predicate>
74  inline
75  void updateHullWithAdaptedStack(Stack aStack, const Point& aNewPoint, const Predicate& aPredicate)
76  {
77  updateHullWithStack( aStack, aNewPoint, aPredicate );
78  }
79 
80  //----------------------------------------------------------------------------
81  template <typename Stack, typename ForwardIterator, typename Predicate>
82  inline
83  void buildHullWithStack(Stack& aStack,
84  const ForwardIterator& itb, const ForwardIterator& ite,
85  const Predicate& aPredicate)
86  {
87  BOOST_CONCEPT_ASSERT(( concepts::CStack<Stack> ));
88  BOOST_CONCEPT_ASSERT(( boost_concepts::ForwardTraversalConcept<ForwardIterator> ));
89  BOOST_CONCEPT_ASSERT(( boost_concepts::ReadableIteratorConcept<ForwardIterator> ));
90 
91  //for all points
92  for(ForwardIterator it = itb; it != ite; ++it)
93  {
94  if(aStack.size() < 2)
95  {
96  aStack.push( *it );
97  }
98  else
99  {
100  //we update the hull so that the predicate returns 'true'
101  //for each sets of three consecutive points
102  updateHullWithStack(aStack, *it, aPredicate);
103  //add new point
104  aStack.push( *it );
105  }
106  }//end for all points
107  }
108 
109  //----------------------------------------------------------------------------
110  template <typename Stack, typename ForwardIterator, typename Predicate>
111  inline
112  void buildHullWithAdaptedStack(Stack aStack,
113  const ForwardIterator& itb, const ForwardIterator& ite,
114  const Predicate& aPredicate)
115  {
116  buildHullWithStack( aStack, itb, ite, aPredicate );
117  }
118 
119  //----------------------------------------------------------------------------
120  template <typename ForwardIterator, typename OutputIterator, typename Predicate>
121  inline
122  void openGrahamScan(const ForwardIterator& itb, const ForwardIterator& ite,
123  OutputIterator res, const Predicate& aPredicate)
124  {
125  BOOST_CONCEPT_ASSERT(( boost_concepts::ForwardTraversalConcept<ForwardIterator> ));
126  BOOST_CONCEPT_ASSERT(( boost_concepts::ReadableIteratorConcept<ForwardIterator> ));
127  typedef typename IteratorCirculatorTraits<ForwardIterator>::Value Point;
128  BOOST_CONCEPT_ASSERT(( boost_concepts::IncrementableIteratorConcept<OutputIterator> ));
129  BOOST_CONCEPT_ASSERT(( boost_concepts::WritableIteratorConcept<OutputIterator,Point> ));
130 
131  //container
132  std::deque<Point> container;
133 
134  //hull computation
135  buildHullWithAdaptedStack( backStack(container), itb, ite, aPredicate );
136 
137  //Copy
138  std::copy(container.begin(), container.end(), res);
139  }
140 
141  //----------------------------------------------------------------------------
142  template <typename ForwardIterator, typename OutputIterator, typename Predicate>
143  inline
144  void closedGrahamScanFromVertex(const ForwardIterator& itb, const ForwardIterator& ite,
145  OutputIterator res, const Predicate& aPredicate)
146  {
147  BOOST_CONCEPT_ASSERT(( boost_concepts::ForwardTraversalConcept<ForwardIterator> ));
148  BOOST_CONCEPT_ASSERT(( boost_concepts::ReadableIteratorConcept<ForwardIterator> ));
149  typedef typename IteratorCirculatorTraits<ForwardIterator>::Value Point;
150  BOOST_CONCEPT_ASSERT(( boost_concepts::IncrementableIteratorConcept<OutputIterator> ));
151  BOOST_CONCEPT_ASSERT(( boost_concepts::WritableIteratorConcept<OutputIterator,Point> ));
152 
153  //container
154  std::deque<Point> container;
155 
156  //hull computation
157  buildHullWithAdaptedStack( backStack(container), itb, ite, aPredicate);
158 
159  //we update the hull to take into account the starting point
160  if ( container.size() > 3 )
161  updateHullWithAdaptedStack( backStack(container), container.front(), aPredicate );
162 
163  std::copy(container.begin(), container.end(), res);
164  }
165 
166  //----------------------------------------------------------------------------
167  template <typename ForwardIterator, typename OutputIterator, typename Predicate>
168  inline
169  void closedGrahamScanFromAnyPoint(const ForwardIterator& itb, const ForwardIterator& ite,
170  OutputIterator res, const Predicate& aPredicate)
171  {
172  BOOST_CONCEPT_ASSERT(( boost_concepts::ForwardTraversalConcept<ForwardIterator> ));
173  BOOST_CONCEPT_ASSERT(( boost_concepts::ReadableIteratorConcept<ForwardIterator> ));
174  typedef typename IteratorCirculatorTraits<ForwardIterator>::Value Point;
175  BOOST_CONCEPT_ASSERT(( boost_concepts::IncrementableIteratorConcept<OutputIterator> ));
176  BOOST_CONCEPT_ASSERT(( boost_concepts::WritableIteratorConcept<OutputIterator,Point> ));
177 
178  //container
179  std::deque<Point> container;
180 
181  //hull computation
182  buildHullWithAdaptedStack( backStack(container), itb, ite, aPredicate);
183 
184  //we update the hull to take into account the starting point
185  if ( container.size() > 3 )
186  {
187  updateHullWithAdaptedStack( backStack(container), container.front(), aPredicate );
188  //we move forward the starting point
189  container.push_back( container.front() );
190  container.pop_front();
191  }
192 
193  //while the last two points and the first one are not correctly oriented, we update the hull
194  while ( (container.size() > 3)&&
195  !aPredicate( *boost::next(container.rbegin()), container.back(), container.front() ) )
196  {
197  updateHullWithAdaptedStack( backStack(container), container.front(), aPredicate );
198  //we move forward the starting point
199  container.push_back( container.front() );
200  container.pop_front();
201  }
202 
203  std::copy(container.begin(), container.end(), res);
204  }
205 
206 
207  //----------------------------------------------------------------------------
208  template <typename ForwardIterator,
209  typename OutputIterator,
210  typename Predicate,
211  typename PolarComparator >
212  inline
213  void grahamConvexHullAlgorithm(const ForwardIterator& itb, const ForwardIterator& ite,
214  OutputIterator res,
215  const Predicate& aPredicate,
216  PolarComparator& aPolarComparator)
217  {
218  BOOST_CONCEPT_ASSERT(( boost_concepts::ForwardTraversalConcept<ForwardIterator> ));
219  BOOST_CONCEPT_ASSERT(( boost_concepts::ReadableIteratorConcept<ForwardIterator> ));
220  typedef typename IteratorCirculatorTraits<ForwardIterator>::Value Point;
221  BOOST_CONCEPT_ASSERT(( boost_concepts::IncrementableIteratorConcept<OutputIterator> ));
222  BOOST_CONCEPT_ASSERT(( boost_concepts::WritableIteratorConcept<OutputIterator,Point> ));
223  BOOST_CONCEPT_ASSERT(( concepts::CPolarPointComparator2D<PolarComparator> ));
224 
225  if ( itb != ite )
226  {
227  //container
228  std::vector<Point> container;
229  std::copy( itb, ite, std::back_inserter( container ) );
230 
231  //find an extremal point
232  //NB: we choose the point of greatest x-coordinate
233  //so that the sort step (by a polar comparator)
234  //returns a weakly externally visible polygon
235  typename std::vector<Point>::iterator itMax
236  = std::max_element( container.begin(), container.end() );
237 
238  //sort around this point
239  aPolarComparator.setPole( *itMax );
240  std::sort( container.begin(), container.end(), aPolarComparator );
241 
242  //scan from an extremal point
243  closedGrahamScanFromVertex( container.begin(), container.end(), res, aPredicate );
244  }
245  }
246 
247  //----------------------------------------------------------------------------
248  template <typename ForwardIterator,
249  typename OutputIterator,
250  typename Predicate>
251  inline
252  void grahamConvexHullAlgorithm(const ForwardIterator& itb, const ForwardIterator& ite,
253  OutputIterator res,
254  const Predicate& aPredicate)
255  {
256  BOOST_CONCEPT_ASSERT(( boost_concepts::ForwardTraversalConcept<ForwardIterator> ));
257  BOOST_CONCEPT_ASSERT(( boost_concepts::ReadableIteratorConcept<ForwardIterator> ));
258  typedef typename IteratorCirculatorTraits<ForwardIterator>::Value Point;
259 
260  //define a default comparator
261  typedef typename Point::Coordinate Integer;
262  typedef AvnaimEtAl2x2DetSignComputer<Integer> DetComputer;
263  typedef functors::PolarPointComparatorBy2x2DetComputer<Point, DetComputer> Comparator;
264  Comparator comparator;
265 
266  //call overloaded version
267  grahamConvexHullAlgorithm( itb, ite, res, aPredicate, comparator );
268  }
269 
270  //----------------------------------------------------------------------------
271  template <typename ForwardIterator,
272  typename OutputIterator,
273  typename Predicate >
274  inline
275  void andrewConvexHullAlgorithm(const ForwardIterator& itb, const ForwardIterator& ite,
276  OutputIterator res,
277  const Predicate& aPredicate )
278  {
279  BOOST_CONCEPT_ASSERT(( boost_concepts::ForwardTraversalConcept<ForwardIterator> ));
280  BOOST_CONCEPT_ASSERT(( boost_concepts::ReadableIteratorConcept<ForwardIterator> ));
281  typedef typename IteratorCirculatorTraits<ForwardIterator>::Value Point;
282  BOOST_CONCEPT_ASSERT(( boost_concepts::IncrementableIteratorConcept<OutputIterator> ));
283  BOOST_CONCEPT_ASSERT(( boost_concepts::WritableIteratorConcept<OutputIterator,Point> ));
284 
285  if ( itb != ite )
286  {
287  //containers
288  std::vector<Point> container;
289  std::copy( itb, ite, std::back_inserter( container ) );
290  std::vector<Point> upperHull, lowerHull;
291 
292  //sort according to the x-coordinate
293  std::sort( container.begin(), container.end() );
294 
295  //lower hull computation
296  openGrahamScan( container.begin(), container.end(), std::back_inserter(lowerHull), aPredicate );
297 
298  //upper hull computation
299  openGrahamScan( container.rbegin(), container.rend(), std::back_inserter(upperHull), aPredicate );
300 
301  //lower hull output
302  typename std::vector<Point>::iterator lowerHullStart = lowerHull.begin();
303  if ( lowerHull.front() == upperHull.back() )
304  lowerHullStart++;
305  std::copy( lowerHullStart, lowerHull.end(), res );
306 
307  //upper hull output
308  typename std::vector<Point>::iterator upperHullStart = upperHull.begin();
309  if ( lowerHull.back() == upperHull.front() )
310  upperHullStart++;
311  std::copy( upperHullStart, upperHull.end(), res );
312  }
313  }
314 
315 
316 
317  template <typename ForwardIterator>
318  inline
319  double computeHullThickness(const ForwardIterator &itb,
320  const ForwardIterator &ite,
321  const ThicknessDefinition &def)
322  {
323  typedef typename std::iterator_traits<ForwardIterator>::value_type TInputPoint;
324  TInputPoint p, q, r;
325  return computeHullThickness(itb, ite, def, p, q, r);
326  }
327 
328 
329  template <typename ForwardIterator,
330  typename TInputPoint>
331  inline
332  double computeHullThickness(const ForwardIterator &itb,
333  const ForwardIterator &ite,
334  const ThicknessDefinition &def,
335  TInputPoint& antipodalEdgeP,
336  TInputPoint& antipodalEdgeQ,
337  TInputPoint& antipodalVertexR){
338  struct WrapIt{
339  WrapIt(const ForwardIterator &b, const ForwardIterator &e): myB(b), myE(e),
340  mySize(std::distance(b, e)){}
341  TInputPoint operator[](unsigned int i){
342  unsigned int j = i%(mySize);
343  return *(myB+j);
344  }
345  ForwardIterator myB;
346  ForwardIterator myE;
347  unsigned int mySize;
348  };
349  WrapIt aConvexHull(itb, ite);
350  double resThickness = std::numeric_limits<double>::max();
351  unsigned int i = 0;
352  unsigned int j = 1;
353  unsigned int size = aConvexHull.mySize;
354  if(size<3)
355  {
356  if(size>0)
357  {
358  antipodalEdgeP = aConvexHull[0];
359  antipodalEdgeQ = aConvexHull[size==1?0:1];
360  antipodalVertexR = aConvexHull[size==1?0:1];
361  }
362  return 0.0;
363  }
364 
365  while(getAngle(aConvexHull[i], aConvexHull[i+1], aConvexHull[j], aConvexHull[j+1]) < M_PI ){
366  j++;
367  }
368  double th = getThicknessAntipodalPair(aConvexHull[i], aConvexHull[i+1], aConvexHull[j], def);
369  if(th<resThickness){
370  resThickness = th;
371  antipodalVertexR = aConvexHull[j];
372  antipodalEdgeP = aConvexHull[i];
373  antipodalEdgeQ = aConvexHull[i+1];
374  }
375  i++;
376  while(i < size){
377  if(getAngle(aConvexHull[i], aConvexHull[i+1],
378  aConvexHull[j], aConvexHull[j+1])<M_PI){
379  j++;
380  }else{
381  th = getThicknessAntipodalPair(aConvexHull[i], aConvexHull[i+1],
382  aConvexHull[j], def);
383  if(th<resThickness){
384  resThickness = th;
385  antipodalVertexR = aConvexHull[j];
386  antipodalEdgeP = aConvexHull[i];
387  antipodalEdgeQ = aConvexHull[i+1];
388  }
389  if(getAngle(aConvexHull[i], aConvexHull[i+1],
390  aConvexHull[j], aConvexHull[j+1])==M_PI){
391 
392  th = getThicknessAntipodalPair(aConvexHull[i], aConvexHull[i+1],
393  aConvexHull[j+1], def);
394  if(th<resThickness){
395  resThickness = th;
396  antipodalVertexR = aConvexHull[j+1];
397  antipodalEdgeP = aConvexHull[i];
398  antipodalEdgeQ = aConvexHull[i+1];
399  }
400  }
401  i++;
402  }
403  }
404  return resThickness;
405  }
406 
407 
408  template<typename TInputPoint>
409  inline
410  double getAngle(const TInputPoint &a, const TInputPoint &b, const TInputPoint &c, const TInputPoint &d){
411  double angle1 = atan2(b[1]-a[1], b[0]-a[0]);
412  double angle2 = atan2(d[1]-c[1], d[0]-c[0]);
413  double r =angle2-angle1;
414  return ( r < 0) ? 2*M_PI+r :r;
415  }
416 
417 
418  template<typename TInputPoint>
419  inline
420  double getThicknessAntipodalPair(const TInputPoint &p, const TInputPoint &q,
421  const TInputPoint &r, const ThicknessDefinition &def){
422  bool isInside;
423  if(def == HorizontalVerticalThickness){
424  double dH = computeHProjDistance(p, q, r, isInside);
425  double dV = computeVProjDistance(p, q, r, isInside);
426  return dH > dV ? dV : dH;
427  }else{
428  return computeEuclideanDistance(p, q, r, isInside);
429  }
430  }
431 
432  template< typename TInputPoint>
433  inline
434  double
435  computeHProjDistance(const TInputPoint &a, const TInputPoint &b, const TInputPoint &c,
436  bool &isInside )
437  {
438  if(a[1]==b[1])
439  {
440  return std::numeric_limits<double>::max();
441  }
442  else
443  {
444  auto k = -(a[0]-b[0])*c[1]-(b[0]*a[1]-a[0]*b[1]);
445  isInside = ( a[1] <= c[1] && c[1] <= b[1] )
446  || ( b[1] <= c[1] && c[1] <= a[1] );
447  return std::abs((k/static_cast<double>(b[1]-a[1])) - c[0]);
448  }
449  }
450 
451  template< typename TInputPoint>
452  inline
453  double
454  computeVProjDistance(const TInputPoint &a, const TInputPoint &b,
455  const TInputPoint &c, bool &isInside )
456  {
457  if(a[0]==b[0])
458  {
459  return std::numeric_limits<double>::max();
460  }
461  else
462  {
463  auto k = -(b[1]-a[1])*c[0]-(b[0]*a[1]-a[0]*b[1]);
464  isInside = ( a[0] <= c[0] && c[0] <= b[0] )
465  || ( b[0] <= c[0] && c[0] <= a[0] );
466  return std::abs((k/static_cast<double>(a[0]-b[0])) - c[1]);
467  }
468  }
469 
470 
471  template< typename TInputPoint>
472  inline
473  double
474  computeEuclideanDistance(const TInputPoint &a, const TInputPoint &b, const TInputPoint &c,
475  bool &isInside )
476  {
477  auto dy = b[1]-a[1];
478  auto dx = a[0]-b[0];
479  auto dc = b[0]*a[1]-a[0]*b[1];
480  auto const pos = (b-a).dot(c-a);
481  isInside = (0 <= pos) && (pos <= (b-a).dot(b-a));
482  return std::abs(dy*c[0]+dx*c[1]+dc)/(std::sqrt(dy*dy+dx*dx));
483  }
484 
485 
486 
487  } // namespace convexHull2D
488  } // namespace functions
489 } // namespace DGtal
490 
491 // //
492 ///////////////////////////////////////////////////////////////////////////////
493 
494 
Generated on Tue Jun 1 2021 17:05:20 for DGtal by doxygen 1.9.1