DGtal  1.4.beta
testConvexHull2D.cpp
Go to the documentation of this file.
1 
31 #include <iostream>
32 #include "DGtal/base/Common.h"
33 #include "DGtal/base/Circulator.h"
34 
35 #include "DGtal/kernel/PointVector.h"
36 #include "DGtal/geometry/tools/Hull2DHelpers.h"
37 #include "DGtal/geometry/tools/PolarPointComparatorBy2x2DetComputer.h"
38 #include "DGtal/geometry/tools/determinant/InHalfPlaneBySimple3x3Matrix.h"
39 #include "DGtal/geometry/tools/determinant/PredicateFromOrientationFunctor2.h"
40 #include "DGtal/io/boards/Board2D.h"
42 
43 using namespace std;
44 using namespace DGtal;
45 
46 
47 
48 
50 // Functions for testing the functions devoted to convex hull computations.
52 
65 template <typename ForwardIterator>
66 bool circularlyEqual(const ForwardIterator& first1, const ForwardIterator& last1,
67  const ForwardIterator& first2, const ForwardIterator& last2 )
68 {
69  ASSERT( first2 != last2 );
70 
71  //find a element of the first range equal to the first element of the second range
72  ForwardIterator start1 = find( first1, last1, *first2 );
73  if (start1 == last1)
74  return false;
75  else
76  {
77  bool areEqual = true; //true if the two ranges are equal up to circular shifts
78  //check whether the two ranges are equal or not
79  Circulator<ForwardIterator> c1(start1, first1, last1);
80  Circulator<ForwardIterator> cEnd1 = c1;
81  Circulator<ForwardIterator> c2(first2, first2, last2);
82  do
83  {
84  if (*c1 != *c2)
85  areEqual = false;
86  else
87  {
88  ++c1;
89  ++c2;
90  }
91  } while ( (c1 != cEnd1) && (areEqual) );
92  return areEqual;
93 
94  }
95 }
96 
97 
102 bool testConvexHull2D()
103 {
104  unsigned int nbok = 0;
105  unsigned int nb = 0;
106 
107  typedef PointVector<2,DGtal::int32_t> Point;
108 
109  trace.beginBlock ( "One simple test..." );
110 
111  vector<Point> data, g, res;
112  //data
113  data.push_back( Point(2,0) );
114  data.push_back( Point(4,0) );
115  data.push_back( Point(0,3) );
116  data.push_back( Point(0,-4) );
117  data.push_back( Point(3,4) );
118  data.push_back( Point(5,0) );
119  data.push_back( Point(4,3) );
120  data.push_back( Point(0,5) );
121  data.push_back( Point(-3,-4) );
122  data.push_back( Point(-5,0) );
123  data.push_back( Point(-4,-3) );
124  data.push_back( Point(0,-5) );
125  data.push_back( Point(3,-4) );
126  data.push_back( Point(4,-3) );
127  data.push_back( Point(-3,4) );
128  data.push_back( Point(-4,3) );
129  //ground truth
130  g.push_back( Point(5,0) );
131  g.push_back( Point(4,3) );
132  g.push_back( Point(3,4) );
133  g.push_back( Point(0,5) );
134  g.push_back( Point(-3,4) );
135  g.push_back( Point(-4,3) );
136  g.push_back( Point(-5,0) );
137  g.push_back( Point(-4,-3) );
138  g.push_back( Point(-3,-4) );
139  g.push_back( Point(0,-5) );
140  g.push_back( Point(3,-4) );
141  g.push_back( Point(4,-3) );
142 
143  //geometric predicate
144  typedef InHalfPlaneBySimple3x3Matrix<Point, DGtal::int64_t> Functor;
145  Functor functor;
146  typedef PredicateFromOrientationFunctor2<Functor> Predicate;
147  Predicate predicate( functor );
148 
149  //namespace
150  using namespace functions::Hull2D;
151 
152  //andrew algorithm
153  trace.info() << " andrew algorithm " << std::endl;
154  andrewConvexHullAlgorithm( data.begin(), data.end(), back_inserter( res ), predicate );
155 
156  copy(res.begin(), res.end(), ostream_iterator<Point>( cout, " " ) );
157  cout << endl;
158 
159  if ( (res.size() == g.size()) &&
160  (circularlyEqual(res.begin(), res.end(), g.begin(), g.end())) )
161  nbok++;
162  nb++;
163  trace.info() << "(" << nbok << "/" << nb << ") " << endl;
164 
165  //graham algorithm
166  res.clear();
167  trace.info() << " graham algorithm " << std::endl;
168  functors::PolarPointComparatorBy2x2DetComputer<Point> comparator;
169  grahamConvexHullAlgorithm( data.begin(), data.end(), back_inserter( res ), predicate, comparator );
170 
171  copy(res.begin(), res.end(), ostream_iterator<Point>( cout, " " ) );
172  cout << endl;
173 
174  if ( (res.size() == g.size()) &&
175  (circularlyEqual(res.begin(), res.end(), g.begin(), g.end())) )
176  nbok++;
177  nb++;
178  trace.info() << "(" << nbok << "/" << nb << ") " << endl;
179 
180  //melkman algorithm
181  res.clear();
182  trace.info() << " melkman algorithm " << std::endl;
183  sort( data.begin(), data.end() );
184  melkmanConvexHullAlgorithm( data.begin(), data.end(), back_inserter( res ), functor );
185 
186  copy(res.begin(), res.end(), ostream_iterator<Point>( cout, " " ) );
187  cout << endl;
188 
189  if ( (res.size() == g.size()) &&
190  (circularlyEqual(res.begin(), res.end(), g.begin(), g.end())) )
191  nbok++;
192  nb++;
193  trace.info() << "(" << nbok << "/" << nb << ") " << endl;
194  // melkman on line construction of convex hull:
195  trace.info() << "on line convex hull construction" << std::endl;
196  DGtal::MelkmanConvexHull<Point, Functor> ch;
197  for(vector<Point>::const_iterator it = data.begin(); it != data.end(); it++)
198  {
199  ch.add(*it);
200  }
201 
202 
203  unsigned int cvSize = 0;
204  for(DGtal::MelkmanConvexHull<Point, Functor>::ConstIterator it = ch.begin(); it != ch.end(); it++, cvSize++)
205  {
206  trace.info() << *it ;
207  };
208  trace.info() << std::endl;
209 
210  if(res.size() == cvSize && (cvSize == ch.size()))
211  nbok++;
212  nb++;
213 
214  trace.info() << "(" << nbok << "/" << nb << ") " << endl;
215  // test copy and [] operator on convex hull:
216  trace.info() << "test copy and [] operator on convex hull:" << std::endl;
217  DGtal::MelkmanConvexHull<Point, Functor> ch2 = ch;
218  unsigned int cvSize2 = 0;
219  for(DGtal::MelkmanConvexHull<Point, Functor>::ConstIterator it = ch.begin(); it != ch.end(); it++, cvSize2++)
220  {
221  trace.info() << *it ;
222  };
223  if(res.size() == cvSize2 && ch[0] == ch2[0])
224  nbok++;
225  nb++;
226  trace.info() << "(" << nbok << "/" << nb << ") " << endl;
227  trace.endBlock();
228 
229 
230  trace.beginBlock ( "Random Tests..." );
231  vector<Point> randomData, res1, res2;
232  const int numberOfPoints = 1000;
233  const int numberOfTries = 50;
234 
235  for (int i = 0; ( (i < numberOfTries)&&(nbok == nb) ); i++)
236  {
237  //new data
238  randomData.clear();
239  res1.clear();
240  res2.clear();
241  for (int j = 0; j < numberOfPoints; j++)
242  randomData.push_back( Point(rand()%256, rand()%256) );
243  //computation
244  andrewConvexHullAlgorithm( randomData.begin(), randomData.end(), back_inserter( res1 ), predicate );
245  grahamConvexHullAlgorithm( randomData.begin(), randomData.end(), back_inserter( res2 ), predicate, comparator );
246  //comparison
247  if ( (res1.size() == res2.size()) &&
248  (circularlyEqual(res1.begin(), res1.end(), res2.begin(), res2.end())) )
249  nbok++;
250  nb++;
251  trace.info() << "(" << nbok << "/" << nb << ") " << endl;
252  //another computation
253  res2.clear();
254  sort( randomData.begin(), randomData.end() );
255  melkmanConvexHullAlgorithm( randomData.begin(), randomData.end(), back_inserter( res2 ), functor );
256  //comparison
257  if ( (res1.size() == res2.size()) &&
258  (circularlyEqual(res1.begin(), res1.end(), res2.begin(), res2.end())) )
259  nbok++;
260  nb++;
261  trace.info() << "(" << nbok << "/" << nb << ") " << endl;
262  }
263 
264  trace.endBlock();
265 
266  return nbok == nb;
267 }
268 
269 
274 bool testConvexHullCompThickness()
275 {
276  unsigned int nbok = 0;
277  unsigned int nb = 0;
278  typedef PointVector<2,DGtal::int32_t> Point;
279  typedef InHalfPlaneBySimple3x3Matrix<Point, DGtal::int32_t> Functor;
280 
281  trace.beginBlock ( "One simple test..." );
282  DGtal::MelkmanConvexHull<Point, Functor> ch;
283  ch.add(Point(0,0));
284  ch.add(Point(11,1));
285  ch.add(Point(12,3));
286  ch.add(Point(8,3));
287  ch.add(Point(4,5));
288  ch.add(Point(2,6));
289  ch.add(Point(1,4));
290 
291  Point antipodalP, antipodalQ, antipodalS;
292  double thicknessE = DGtal::functions::Hull2D::computeHullThickness(ch.begin(), ch.end(),
293  DGtal::functions::Hull2D::EuclideanThickness,
294  antipodalP,
295  antipodalQ,
296  antipodalS);
297  double thicknessHV = DGtal::functions::Hull2D::computeHullThickness(ch.begin(), ch.end(),
298  DGtal::functions::Hull2D::HorizontalVerticalThickness,
299  antipodalP, antipodalQ, antipodalS);
300  double thicknessHVb = DGtal::functions::Hull2D::computeHullThickness(ch.begin(), ch.end(),
301  DGtal::functions::Hull2D::HorizontalVerticalThickness);
302 
303 
304  Board2D aBoard;
305  for(DGtal::MelkmanConvexHull<Point, Functor>::ConstIterator it = ch.begin(); it != ch.end(); it++){
306  if(it != ch.end()-1)
307  aBoard.drawLine((*it)[0], (*it)[1], (*(it+1))[0], (*(it+1))[1]);
308  else{
309  aBoard.drawLine((*it)[0], (*it)[1], (*(ch.begin()))[0], (*(ch.begin()))[1]);
310  }
311  aBoard << *it;
312  }
313 
314  aBoard.setPenColor(DGtal::Color::Red);
315  aBoard.drawCircle( antipodalS[0], antipodalS[1], 1.0) ;
316  aBoard.setPenColor(DGtal::Color::Blue);
317  aBoard.drawCircle(antipodalP[0], antipodalP[1], 1.0);
318  aBoard.drawCircle(antipodalQ[0], antipodalQ[1], 1.0);
319 
320  aBoard.drawLine(antipodalP[0], antipodalP[1], antipodalQ[0], antipodalQ[1]);
321  double awaitedThHV = DGtal::functions::Hull2D::getThicknessAntipodalPair(Point(0,0), Point(11,1), Point(2,6), DGtal::functions::Hull2D::HorizontalVerticalThickness );
322  double awaitedThE = DGtal::functions::Hull2D::getThicknessAntipodalPair(Point(0,0), Point(11,1), Point(2,6), DGtal::functions::Hull2D::EuclideanThickness );
323  trace.info() << "Thickness HV = " << thicknessHV << std::endl;
324  trace.info() << "Expected Thickness HV = " << awaitedThHV << std::endl;
325  trace.info() << "Thickness Euclidean = " << thicknessE << std::endl;
326  trace.info() << "Expected Euclidean Thickness = " << awaitedThE << std::endl;
327  aBoard.saveEPS("testConvexHull2D_Thickness.eps");
328 
329  // testing tickness after changing begin points.
330  std::vector<Z2i::RealPoint> hull {
331  {804.56832227024199, -68.471176393526704},
332  {804.96020257363512, -69.494933490400683},
333  {805.35232247974727, -70.519316530915930},
334  {825.15497274857830, -114.34711249674335},
335  {828.67425867587508, -121.69593677670120},
336  {829.05943325037651, -122.39546351563243},
337  {829.42436256177677, -122.73833140123513},
338  {827.82353937509288, -118.87280445059109},
339  {811.06291230576301, -82.124832029926566},
340  {806.83483124583609, -72.890457996792406},
341  {804.92531970702396, -68.803808853026638},
342  {804.55092650722997, -68.125792709291431},
343  };
344 
345  double th = DGtal::functions::Hull2D::computeHullThickness(hull.begin(), hull.end(),
346  DGtal::functions::Hull2D::HorizontalVerticalThickness);
347  std::rotate(hull.begin(), hull.begin() + 1, hull.end());
348  double th2 = DGtal::functions::Hull2D::computeHullThickness(hull.begin(), hull.end(),
349  DGtal::functions::Hull2D::HorizontalVerticalThickness);
350  std::rotate(hull.begin(), hull.begin() + 5, hull.end());
351  double th3 = DGtal::functions::Hull2D::computeHullThickness(hull.begin(), hull.end(),
352  DGtal::functions::Hull2D::HorizontalVerticalThickness);
353  trace.info() << "Thickness (before change init point) = " << th << std::endl;
354  trace.info() << "Thickness (after change init point) = " << th2 << std::endl;
355  trace.info() << "Thickness (after change init point) = " << th3 << std::endl;
356 
357  nbok += thicknessHV == awaitedThHV && thicknessE == awaitedThE &&
358  thicknessHVb == thicknessHV && abs(th - th2) < 0.000001 &&
359  th - 0.604414 < 0.000001 && abs(th - th3) < 0.000001 ;
360  nb++;
361  return nb==nbok;
362 }
363 
364 
365 
366 
368 // Standard services - public :
369 
370 int main( int argc, char** argv )
371 {
372  trace.beginBlock ( "Testing convex hull function" );
373  trace.info() << "Args:";
374  for ( int i = 0; i < argc; ++i )
375  trace.info() << " " << argv[ i ];
376  trace.info() << endl;
377 
378  bool res = testConvexHull2D() &&testConvexHullCompThickness();
379  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
380  trace.endBlock();
381  return res ? 0 : 1;
382 }
383 // //
DGtal::Board2D
Aim: This class specializes a 'Board' class so as to display DGtal objects more naturally (with <<)....
Definition: Board2D.h:71
DGtal::Circulator
Aim: Provides an adapter for classical iterators that can iterate through the underlying data structu...
Definition: Circulator.h:86
DGtal::Color::Red
static const Color Red
Definition: Color.h:416
DGtal::Color::Blue
static const Color Blue
Definition: Color.h:419
DGtal::InHalfPlaneBySimple3x3Matrix
Aim: Class that implements an orientation functor, ie. it provides a way to compute the orientation o...
Definition: InHalfPlaneBySimple3x3Matrix.h:92
DGtal::MelkmanConvexHull
Aim: This class implements the on-line algorithm of Melkman for the computation of the convex hull of...
Definition: MelkmanConvexHull.h:90
DGtal::MelkmanConvexHull::ConstIterator
std::deque< Point >::const_iterator ConstIterator
Definition: MelkmanConvexHull.h:123
DGtal::PointVector
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:593
DGtal::PredicateFromOrientationFunctor2
Aim: Small adapter to models of COrientationFunctor2. It is a model of concepts::CPointPredicate....
Definition: PredicateFromOrientationFunctor2.h:82
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::Trace::emphase
std::ostream & emphase()
DGtal::Trace::info
std::ostream & info()
DGtal::Trace::endBlock
double endBlock()
DGtal::functors::PolarPointComparatorBy2x2DetComputer
Aim: Class that implements a binary point predicate, which is able to compare the position of two giv...
Definition: PolarPointComparatorBy2x2DetComputer.h:80
LibBoard::Board::setPenColor
Board & setPenColor(const DGtal::Color &color)
Definition: Board.cpp:297
LibBoard::Board::drawLine
void drawLine(double x1, double y1, double x2, double y2, int depthValue=-1)
Definition: Board.cpp:367
LibBoard::Board::drawCircle
void drawCircle(double x, double y, double radius, int depthValue=-1)
Definition: Board.cpp:450
LibBoard::Board::saveEPS
void saveEPS(const char *filename, PageSize size=Board::BoundingBox, double margin=10.0) const
Definition: Board.cpp:804
DGtal::functions::Hull2D::getThicknessAntipodalPair
double getThicknessAntipodalPair(const TInputPoint &p, const TInputPoint &q, const TInputPoint &r, const ThicknessDefinition &def)
DGtal::functions::Hull2D::computeHullThickness
double computeHullThickness(const ForwardIterator &itb, const ForwardIterator &ite, const ThicknessDefinition &def)
Procedure to compute the convex hull thickness given from different definitions (Horizontal/vertical ...
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::trace
Trace trace
Definition: Common.h:153
Point
MyPointD Point
Definition: testClone2.cpp:383
ch
DGtal::MelkmanConvexHull< Point, Functor > ch
Definition: testConvexHull2DReverse.cpp:52
Functor
InHalfPlaneBySimple3x3Matrix< Point, double > Functor
Definition: testConvexHull2DReverse.cpp:49
thicknessE
const double thicknessE
Definition: testConvexHull2DReverse.cpp:77
thicknessHV
const double thicknessHV
Definition: testConvexHull2DReverse.cpp:73
circularlyEqual
bool circularlyEqual(const ForwardIterator &first1, const ForwardIterator &last1, const ForwardIterator &first2, const ForwardIterator &last2)
Definition: testConvexHull2D.cpp:66
main
int main(int argc, char **argv)
Definition: testConvexHull2D.cpp:370
testConvexHullCompThickness
bool testConvexHullCompThickness()
Definition: testConvexHull2D.cpp:274
testConvexHull2D
bool testConvexHull2D()
Definition: testConvexHull2D.cpp:102