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