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