DGtal 1.3.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
34#include "DGtal/kernel/PointVector.h"
35//////////////////////////////////////////////////////////////////////////////
36
37///////////////////////////////////////////////////////////////////////////////
38// Implementation of inline functions //
39
40namespace 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 {
339 struct WrapIt{
340 WrapIt(const ForwardIterator &b, const ForwardIterator &e): myB(b), myE(e),
341 mySize(std::distance(b, e)){}
342 TInputPoint operator[](unsigned int i){
343 unsigned int j = i%(mySize);
344 return *(myB+j);
345 }
346 ForwardIterator myB;
347 ForwardIterator myE;
348 unsigned int mySize;
349 };
350 WrapIt aConvexHull(itb, ite);
351 double resThickness = std::numeric_limits<double>::max();
352 unsigned int i = 0;
353 unsigned int j = 1;
354 unsigned int size = aConvexHull.mySize;
355 if(size<3)
356 {
357 if(size>0)
358 {
359 antipodalEdgeP = aConvexHull[0];
360 antipodalEdgeQ = aConvexHull[size==1?0:1];
361 antipodalVertexR = aConvexHull[size==1?0:1];
362 }
363 return 0.0;
364 }
365 while(getAngle(aConvexHull[i], aConvexHull[i+1],
366 aConvexHull[j], aConvexHull[j+1]) < M_PI - angleTolerance ){
367 j++;
368 }
369 double th = getThicknessAntipodalPair(aConvexHull[i], aConvexHull[i+1], aConvexHull[j], def);
370 if(th<resThickness){
371 resThickness = th;
372 antipodalVertexR = aConvexHull[j];
373 antipodalEdgeP = aConvexHull[i];
374 antipodalEdgeQ = aConvexHull[i+1];
375 }
376 i++;
377 while(i < size){
378 if(getAngle(aConvexHull[i], aConvexHull[i+1],
379 aConvexHull[j], aConvexHull[j+1])< M_PI - angleTolerance){
380 j++;
381 }else{
382 th = getThicknessAntipodalPair(aConvexHull[i], aConvexHull[i+1],
383 aConvexHull[j], def);
384 if(th<resThickness){
385 resThickness = th;
386 antipodalVertexR = aConvexHull[j];
387 antipodalEdgeP = aConvexHull[i];
388 antipodalEdgeQ = aConvexHull[i+1];
389 }
390 if(isCoLinearOpp(aConvexHull[i], aConvexHull[i+1],
391 aConvexHull[j], aConvexHull[j+1]) ){
392
393 th = getThicknessAntipodalPair(aConvexHull[i], aConvexHull[i+1],
394 aConvexHull[j+1], def);
395 if(th<resThickness){
396 resThickness = th;
397 antipodalVertexR = aConvexHull[j+1];
398 antipodalEdgeP = aConvexHull[i];
399 antipodalEdgeQ = aConvexHull[i+1];
400 }
401 }
402 i++;
403 }
404 }
405 return resThickness;
406 }
407
408
409 template<typename TInputPoint>
410 inline
411 double getAngle(const TInputPoint &a, const TInputPoint &b, const TInputPoint &c, const TInputPoint &d)
412 {
413 TInputPoint v1 (b-a);
414 TInputPoint v2 (d-c);
415 double det = (v1[0]*v2[1] - v1[1]*v2[0])/ (v1.norm()*v2.norm());
416 double dot = (v1[0]*v2[0] + v1[1]*v2[1])/ (v1.norm()*v2.norm());
417 double ang = atan2(det, dot);
418 return ( ang < -angleTolerance) ? 2*M_PI+ang : ang;
419 }
420
421
422 template<typename TInputPoint>
423 inline
424 bool isCoLinearOpp(const TInputPoint& a, const TInputPoint& b,
425 const TInputPoint& c,const TInputPoint& d)
426 {
427 double ang = getAngle(a, b, c, d);
428 return ang < M_PI + angleTolerance &&
429 ang > M_PI - angleTolerance ;
430 }
431
432
433
434 template<typename TInputPoint>
435 inline
436 double getThicknessAntipodalPair(const TInputPoint &p, const TInputPoint &q,
437 const TInputPoint &r, const ThicknessDefinition &def){
438 bool isInside;
439 if(def == HorizontalVerticalThickness){
440 double dH = computeHProjDistance(p, q, r, isInside);
441 double dV = computeVProjDistance(p, q, r, isInside);
442 return dH > dV ? dV : dH;
443 }else{
444 return computeEuclideanDistance(p, q, r, isInside);
445 }
446 }
447
448 template< typename TInputPoint>
449 inline
450 double
451 computeHProjDistance(const TInputPoint &a, const TInputPoint &b, const TInputPoint &c,
452 bool &isInside )
453 {
454 if(a[1]==b[1])
455 {
456 return std::numeric_limits<double>::max();
457 }
458 else
459 {
460 auto k = -(a[0]-b[0])*c[1]-(b[0]*a[1]-a[0]*b[1]);
461 isInside = ( a[1] <= c[1] && c[1] <= b[1] )
462 || ( b[1] <= c[1] && c[1] <= a[1] );
463 return std::abs((k/static_cast<double>(b[1]-a[1])) - c[0]);
464 }
465 }
466
467 template< typename TInputPoint>
468 inline
469 double
470 computeVProjDistance(const TInputPoint &a, const TInputPoint &b,
471 const TInputPoint &c, bool &isInside )
472 {
473 if(a[0]==b[0])
474 {
475 return std::numeric_limits<double>::max();
476 }
477 else
478 {
479 auto k = -(b[1]-a[1])*c[0]-(b[0]*a[1]-a[0]*b[1]);
480 isInside = ( a[0] <= c[0] && c[0] <= b[0] )
481 || ( b[0] <= c[0] && c[0] <= a[0] );
482 return std::abs((k/static_cast<double>(a[0]-b[0])) - c[1]);
483 }
484 }
485
486
487 template< typename TInputPoint>
488 inline
489 double
490 computeEuclideanDistance(const TInputPoint &a, const TInputPoint &b, const TInputPoint &c,
491 bool &isInside )
492 {
493 auto dy = b[1]-a[1];
494 auto dx = a[0]-b[0];
495 auto dc = b[0]*a[1]-a[0]*b[1];
496 auto const pos = (b-a).dot(c-a);
497 isInside = (0 <= pos) && (pos <= (b-a).dot(b-a));
498 return std::abs(dy*c[0]+dx*c[1]+dc)/(std::sqrt(dy*dy+dx*dx));
499 }
500
501
502
503 } // namespace convexHull2D
504 } // namespace functions
505} // namespace DGtal
506
507// //
508///////////////////////////////////////////////////////////////////////////////
509
510
Generated on Sun Nov 27 2022 15:15:39 for DGtal by doxygen 1.9.5