DGtal 1.4.0
Loading...
Searching...
No Matches
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
42namespace 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