DGtal 2.1.0
Loading...
Searching...
No Matches
testAffineBasis.cpp
Go to the documentation of this file.
1
31#include <iostream>
32#include <vector>
33#include <algorithm>
34#include <random>
35#include "DGtal/base/Common.h"
36#include "DGtal/kernel/SpaceND.h"
37#include "DGtal/geometry/tools/AffineGeometry.h"
38#include "DGtal/geometry/tools/AffineBasis.h"
39#include "DGtalCatch.h"
41
42using namespace std;
43using namespace DGtal;
44
45std::random_device rd;
46std::mt19937 g(rd());
47
48template < typename Point >
49std::vector< Point >
50makeRandomVectors( int nb, int amplitude )
51{
52 std::uniform_int_distribution<int> U(-amplitude, amplitude);
53 std::vector< Point > P;
54 for ( auto n = 0; n < nb; ++n )
55 {
56 Point A;
57 for ( auto i = 0; i < Point::dimension; i++ )
58 A[ i ] = U( g );
59 P.push_back( A );
60 }
61 return P;
62}
63
64template < typename Point >
65std::vector< Point >
66makeRandomLatticePointsFromDirVectors( int nb, const vector< Point>& V )
67{
68 std::uniform_int_distribution<int> U(-10, 10);
69 vector< Point > P;
70 int n = V[0].size();
71 int m = V.size();
72 Point A;
73 for ( auto i = 0; i < n; i++ )
74 A[ i ] = U( g );
75 P.push_back( A );
76 for ( auto k = 0; k < nb; k++ )
77 {
78 Point B = A;
79 for ( auto i = 0; i < m; i++ )
80 {
81 int l = U( g );
82 B += l * V[ i ];
83 }
84 P.push_back( B );
85 }
86 std::shuffle( P.begin(), P.end(), g );
87 return P;
88}
89
90template <typename TPoint>
91bool sameDirection( const TPoint& a, const TPoint& b )
92{
93 return ( a[ 0 ] == b[ 0 ] ) ? ( a == b ) : ( a == -b );
94}
95
96template <typename TPoint>
97std::pair<TPoint, TPoint>
98boundingBox( const std::vector< TPoint >& v )
99{
100 TPoint l, u;
101 if ( ! v.empty() )
102 {
103 l = u = v[ 0 ];
104 for ( auto i = 1; i < v.size(); i++ )
105 {
106 l = l.inf( v[ i ] );
107 u = u.sup( v[ i ] );
108 }
109 }
110 return std::make_pair( l, u );
111}
112
113template <typename T, typename C>
114T
115determinant( const PointVector<2, T, C>& u,
116 const PointVector<2, T, C>& v )
117{
118 return u [ 0 ] * v [ 1 ] - u [ 1 ] * v [ 0 ];
119}
120
121template <typename T, typename C>
122T
123determinant( const PointVector<3, T, C>& u,
124 const PointVector<3, T, C>& v,
125 const PointVector<3, T, C>& w )
126{
127 typedef PointVector<2,T,C> P2;
128 return u[ 0 ] * determinant( P2( v[1], v[2] ), P2( w[1], w[2] ) )
129 - u[ 1 ] * determinant( P2( u[0], u[2] ), P2( w[0], w[2] ) )
130 + u[ 2 ] * determinant( P2( u[0], u[1] ), P2( v[0], v[1] ) );
131}
132
134// Functions for testing class AffineBasis in 2D.
136
137SCENARIO( "AffineBasis< Point2i > unit tests", "[affine_basis][2i]" )
138{
139 typedef SpaceND< 2, int > Space;
140 typedef Space::Point Point;
141 typedef AffineBasis< Point > Basis;
142 GIVEN( "Given B = (0,0) + { (8,2), (-4,-1), (-8,-2), (16,4), (200,50) } of affine dimension 1" ) {
143 Point o( 0, 0 );
144 std::vector<Point> X
145 = { Point(8,2), Point(-4,-1), Point(-8,-2), Point(16,4), Point(200,50) };
146 Basis B( o, X, Basis::Type::ECHELON_REDUCED );
147 THEN( "When reduced, it has dimension 1" ) {
148 CAPTURE( B.basis() );
149 REQUIRE( B.dimension() == 1 );
150 }
151 THEN( "When reduced, it is the vector (-4,-1) or (4,1)" ) {
152 Point b0 = B.basis()[ 0 ];
153 CAPTURE( b0 );
154 REQUIRE( ((b0 == Point(-4,-1)) || (b0 == Point(4,1))) );
155 }
156 }
157}
158
160// Functions for testing class AffineBasis in 3D.
162
163SCENARIO( "AffineBasis< Z3 > LLL tests", "[affine_basis][3d][LLL]" )
164{
165 typedef SpaceND< 3, int> Space;
166 typedef Space::Point Point;
167 typedef AffineGeometry< Point > Affine;
168 typedef AffineBasis< Point > Basis;
169
170 WHEN( "Using basis B = (1, 0, -2) (1, 0, -1)" ) {
171 std::vector< Point > b = { Point( 0, 0, 0), Point(1, 0, -2), Point(1, 0, -1) };
172 const auto [ o, B ] = Affine::affineBasis( b );
173 Point e = functions::computeIndependentVector( B );
174 Basis AB( b, Basis::Type::LLL_REDUCED );
175 const auto [ d, L, r ] = AB.decomposeVector( e );
176 CAPTURE( B );
177 CAPTURE( AB.basis() );
178 CAPTURE( d );
179 CAPTURE( L );
180 CAPTURE( r );
181 CAPTURE( e );
182 REQUIRE( r != Point::zero );
183 }
184}
185
186
188// Functions for testing class AffineBasis in 4D.
190
191SCENARIO( "AffineBasis< Point4i > unit tests", "[affine_basis][4i]" )
192{
193 typedef SpaceND< 4, int > Space;
194 typedef Space::Point Point;
195 typedef AffineBasis< Point > Basis;
196 GIVEN( "Given X a set of randomly generated points by adding linear combinations of 1 lattice vectors" ) {
197 std::vector< Point > V = { Point{ 3, 1, 0, 2 } };
198 auto X = makeRandomLatticePointsFromDirVectors( 20, V );
199 Basis B( X, Basis::Type::ECHELON_REDUCED );
200 THEN( "When reduced, it has dimension 1" ) {
201 CAPTURE( B.basis() );
202 REQUIRE( B.dimension() == 1 );
203 }
204 THEN( "When reduced, it is the vector V[0] or -V[0]" ) {
205 CAPTURE( V );
206 Point b0 = B.basis()[ 0 ];
207 CAPTURE( b0 );
208 REQUIRE( ((b0 == V[0]) || (b0 == -V[0])) );
209 }
210 }
211 GIVEN( "Given X a set of randomly generated points by adding linear combinations of 2 lattice vectors" ) {
212 std::vector< Point > V = { Point{ 3, 1, 0, 2 }, Point{ -2, -1, 2, 7 } };
213 auto X = makeRandomLatticePointsFromDirVectors( 20, V );
214 Basis B( X, Basis::Type::SHORTEST_ECHELON_REDUCED );
215 THEN( "When reduced, basis has dimension 2" ) {
216 CAPTURE( B.basis() );
217 REQUIRE( B.dimension() == 2 );
218 }
219 THEN( "When reduced, basis spans vectors V[i] or -V[i]" ) {
220 CAPTURE( V );
221 Point b0 = B.basis()[ 0 ];
222 Point b1 = B.basis()[ 1 ];
223 CAPTURE( b0 );
224 CAPTURE( b1 );
225 REQUIRE( B.isParallel( V[ 0 ] ) );
226 REQUIRE( B.isParallel( V[ 1 ] ) );
227 }
228 THEN( "every point of X has rational coordinates with no remainder" ) {
229 unsigned int nb_ok = 0;
230 for ( auto p : X )
231 {
232 const auto [d, lambda, rem ] = B.decompose( p );
233 nb_ok += ( rem == Point::zero ) ? : 1;
234 if ( rem != Point::zero )
235 std::cout << "p=" << p << " d=" << d
236 << " lambda=" << lambda << " rem=" << rem << "\n";
237 }
238 REQUIRE( nb_ok == X.size() );
239 }
240 THEN( "every lattice point can be written as a linear combination" ) {
241 auto Y = makeRandomVectors<Point>( 20, 10 );
242 unsigned int nb_ok = 0;
243 for ( auto y : Y )
244 {
245 const auto p = y + B.first;
246 const auto [d, lambda, rem ] = B.decompose( p );
247 auto q = B.recompose( d, lambda, rem );
248 nb_ok += ( p == q ) ? : 1;
249 if ( p != q )
250 std::cout << "p=" << p << " d=" << d
251 << " lambda=" << lambda << " rem=" << rem
252 << " q=" << q << "\n";
253 }
254 REQUIRE( nb_ok == Y.size() );
255 }
256 }
257}
258
259SCENARIO( "AffineBasis< Point4i > projection tests", "[affine_basis][4i][4d]" )
260{
261 typedef SpaceND< 4, int64_t > Space;
262 typedef Space::Point Point;
263 typedef SpaceND< 2, int64_t > Space2;
264 typedef Space2::Point PPoint;
265 typedef AffineBasis< Point > Basis;
266 typedef Space::RealPoint RealPoint;
267 typedef Space2::RealPoint PRealPoint;
268 typedef AffineBasis< RealPoint > RealBasis;
269
270 GIVEN( "Given X a set of randomly generated points by adding linear combinations of 2 lattice vectors, and Y the same set but with real coordinates" ) {
271 std::vector< Point > V = { Point{ 3, 4, 0, 2 }, Point{ -2, -1, 5, -7 } };
272 auto X = makeRandomLatticePointsFromDirVectors( 20, V );
273 Basis B( X, Basis::Type::ECHELON_REDUCED );
274 std::vector< RealPoint > Y( X.size() );
275 for ( auto i = 0; i < Y.size(); i++ )
276 Basis::transform( Y[ i ], X[ i ] );
277 RealBasis RB( Y, RealBasis::Type::ECHELON_REDUCED );
278 std::vector< PPoint > pX;
279 std::vector< PRealPoint > pY;
280 auto lcm = B .projectPoints( pX, X );
281 auto rlcm = RB.projectPoints( pY, Y );
282 auto rbox = boundingBox( pY );
283 double factor = ( rbox.second - rbox.first ).squaredNorm();
284 THEN( "When reduced, their affine bases has same dimension 2" ) {
285 CAPTURE( B.basis() );
286 CAPTURE( RB.basis() );
287 REQUIRE( B.dimension() == 2 );
288 REQUIRE( RB.dimension() == 2 );
289 }
290 THEN( "Their projections have the same geometry (i.e. orientations within points)" ) {
291 CAPTURE( lcm );
292 CAPTURE( rlcm );
293 CAPTURE( pX );
294 CAPTURE( pY );
295 REQUIRE( pX.size() == pY.size() );
296 // Computing arbitrary determinants between triplets of points
297 const std::size_t nb = 2000;
298 std::size_t nb_ok = 0;
299 const double eps = 1e-12 * factor;
300 for ( auto i = 0; i < nb; i++ )
301 {
302 const std::size_t j = rand() % pX.size();
303 const std::size_t k = rand() % pX.size();
304 const std::size_t l = rand() % pX.size();
305 const auto u = pX[ k ] - pX[ j ];
306 const auto v = pX[ l ] - pX[ j ];
307 const auto ru = pY[ k ] - pY[ j ];
308 const auto rv = pY[ l ] - pY[ j ];
309 const auto det = u [ 0 ] * v [ 1 ] - u [ 1 ] * v [ 0 ];
310 const auto rdet = ru[ 0 ] * rv[ 1 ] - ru[ 1 ] * rv[ 0 ];
311 if ( rdet > eps ) nb_ok += ( det > 0 ) ? 1 : 0;
312 else if ( rdet < -eps ) nb_ok += ( det < 0 ) ? 1 : 0;
313 else nb_ok += ( det == 0 ) ? 1 : 0;
314 if ( nb_ok != i+1 )
315 {
316 std::cout << "eps=" << eps << " factor=" << factor << "\n";
317 std::cout << "u =" << u << " v =" << v << " det =" << det << "\n";
318 std::cout << "ru=" << ru << " rv=" << rv << " rdet=" << rdet << "\n";
319 break;
320 }
321 }
322 REQUIRE( nb_ok == nb );
323 }
324 }
325}
326
327
329// Functions for testing class AffineBasis in 5D.
331
332SCENARIO( "AffineBasis< Point5i > unit tests", "[affine_basis][5i]" )
333{
334 typedef SpaceND< 5, BigInteger > Space;
335 typedef Space::Point Point;
336 typedef AffineGeometry< Point > Affine;
337 typedef AffineBasis< Point > Basis;
338
339 unsigned int nb = 0;
340 unsigned int nb_s_parallel_l = 0;
341 unsigned int nb_equal_big = 0;
342 unsigned int nb_equal_N = 0;
343 unsigned int nb_equal_Nr = 0;
344 for ( auto n = 0; n < 10; n++ )
345 {
346 std::vector< Point > X = makeRandomVectors<Point>( 5, 100 );
347 std::vector< Point > Y;
348 for ( auto i = 1; i < X.size(); i++ ) Y.push_back( X[ i ] - X[ 0 ] );
349 Basis B ( X[ 0 ], Y, Basis::Type::LLL_REDUCED );
350 Basis Br( X[ 0 ], Y, Basis::Type::ECHELON_REDUCED );
351 if ( functions::computeAffineDimension( X ) != 4 ) continue;
352 auto N = functions::computeOrthogonalVectorToBasis( B.basis() );
353 auto Nr = functions::computeOrthogonalVectorToBasis( Br.basis() );
354 auto N_big = Affine::template orthogonalVector<BigInteger>( B.basis() );
355 auto Nr_big = Affine::template orthogonalVector<BigInteger>( Br.basis() );
356 auto N_cast = N_big;
357 auto Nr_cast = Nr_big;
358 for ( auto i = 0; i < 5; i++ ) N_cast[ i ] = N[ i ];
359 for ( auto i = 0; i < 5; i++ ) Nr_cast[ i ] = Nr[ i ];
360 nb_equal_big += sameDirection( N_big, Nr_big ) ? 1 : 0;
361 nb_equal_N += sameDirection( N_cast, N_big ) ? 1 : 0;
362 nb_equal_Nr += sameDirection( Nr_cast, Nr_big ) ? 1 : 0;
363 nb += 1;
364 bool s_parallel_l = Br.isParallel( B );
365 nb_s_parallel_l += s_parallel_l ? 1 : 0;
366 if ( ! s_parallel_l )
367 {
368 std::cout << "dim(B)=" << B.dimension() << " dim(Br)=" << Br.dimension() << "\n";
369 std::cout << "* Br is not // to B:\n" << Br << "\n";
370 for ( auto v : B.basis() )
371 {
372 bool parallel = Br.isParallel( v );
373 std::cout << " v=" << v << " //=" << ( parallel ? "True" : "False" ) << "\n";
374 if ( ! parallel )
375 {
376 const auto [d, lambda, r] = B.decomposeVector( v );
377 std::cout << "d=" << d << " L=" << lambda << " r=" << r << "\n";
378 }
379 }
380 }
381 // std::cout << B << "\n" << Br << "\n";
382 // std::cout << "N =" << N << " Nr =" << Nr << "\n";
383 // std::cout << "N_big=" << N_big << " Nr_big=" << Nr_big << "\n";
384 }
385 THEN( "Normals with big integers are the same" ) {
386 REQUIRE( nb_equal_big == nb );
387 }
388 THEN( "There is less overflow in orthogonal vector computation using non reduced basis" ) {
389 REQUIRE( nb_equal_N >= nb_equal_Nr );
390 }
391 THEN( "The LLL-reduced basis and the echelon-reduced basis are parallel" ) {
392 REQUIRE( nb_s_parallel_l == nb );
393 }
394}
395
396
398// Functions for testing class AffineBasis in 10D.
400
401SCENARIO( "AffineBasis< Z10 > LLL tests", "[affine_basis][10d][LLL]" )
402{
403 typedef SpaceND< 10, int64_t> Space;
404 typedef Space::Point Point;
405 typedef AffineBasis< Point > Basis;
406
407
408 WHEN( "Using basis is unimodular" ) {
409 // when the matrix is unimodular, outputs canonic vectors.
410 std::vector< Point > B = {
411 Point{ -3, 10, 47, 61, -53, -126, 713, 601,-1476, 1569},
412 Point{ 2, -7, -33, -43, 37, 89, -502, -425, 1047,-1103},
413 Point{ -3, 11, 53, 69, -59, -142, 800, 677,-1663, 1764},
414 Point{ 1, -9, -48, -63, 52, 130, -727, -623, 1543,-1604},
415 Point{ -2, 9, 48, 63, -49, -124, 680, 583,-1409, 1533},
416 Point{ 5, -25, -118, -163, 113, 334,-1761,-1595, 4030,-3838},
417 Point{ -3, 17, 85, 118, -84, -245, 1297, 1173,-2974, 2824},
418 Point{ 5, -24, -119, -156, 126, 321,-1782,-1534, 3799,-3921},
419 Point{ 2, -10, -44, -65, 42, 137, -699, -659, 1713,-1489},
420 Point{ 1, -5, -23, -27, 33, 64, -405, -315, 784, -868}
421 };
422 Point o = Point::zero;
423 Basis S( o, B, Basis::Type::SHORTEST_ECHELON_REDUCED );
424 Basis L( o, B, Basis::Type::LLL_REDUCED );
425 THEN( "The LLL-reduced basis is canonic" ) {
426 const auto& V = L.basis();
427 unsigned int nb = 0;
428 unsigned int nbok = 0;
429 for ( auto i = 0; i < V.size(); i++ )
430 {
431 nbok += V[ i ].norm1() == 1 ? : 0;
432 nb++;
433 }
434 CAPTURE( L.basis() );
435 REQUIRE( nbok == nb );
436 }
437 THEN( "The echelon-reduced basis is triangular" ) {
438 const auto& V = S.basis();
439 unsigned int nb = 0;
440 unsigned int nbok = 0;
441 for ( auto i = 0; i < V.size(); i++ )
442 for ( auto j = 0; j < i; j++ )
443 {
444 nbok += ( V[ i ][ j ] == 0 ) ? 1 : 0;
445 nb++;
446 }
447 CAPTURE( S.basis() );
448 REQUIRE( nbok == nb );
449 }
450 THEN( "The LLL-reduced basis is parallel to the echelon-reduced basis" ) {
451 REQUIRE( S.isParallel( L ) );
452 }
453 }
454}
455
456
457SCENARIO( "AffineBasis< Point5i > projection tests", "[affine_basis][5i][5d]" )
458{
459 typedef SpaceND< 5, int64_t > Space;
460 typedef Space::Point Point;
461 typedef AffineBasis< Point > Basis;
462 typedef Space::RealPoint RealPoint;
463 typedef AffineBasis< RealPoint > RealBasis;
464
465 GIVEN( "Given X a set of randomly generated points by adding linear combinations of 2 lattice vectors, and Y the same set but with real coordinates" ) {
466 typedef SpaceND< 2, int64_t > Space2;
467 typedef Space2::Point PPoint;
468 typedef Space2::RealPoint PRealPoint;
469 std::vector< Point > V = { Point{ 3, 4, 0, 2, -5 }, Point{ -2, -1, 5, -7, 1 } };
470 auto X = makeRandomLatticePointsFromDirVectors( 20, V );
471 Basis B( X, Basis::Type::ECHELON_REDUCED );
472 std::vector< RealPoint > Y( X.size() );
473 for ( auto i = 0; i < Y.size(); i++ )
474 Basis::transform( Y[ i ], X[ i ] );
475 RealBasis RB( Y, RealBasis::Type::ECHELON_REDUCED, 0.99, 1e-10 );
476 std::vector< PPoint > pX;
477 std::vector< PRealPoint > pY;
478 auto lcm = B .projectPoints( pX, X );
479 auto rlcm = RB.projectPoints( pY, Y );
480 auto rbox = boundingBox( pY );
481 double factor = ( rbox.second - rbox.first ).squaredNorm();
482 THEN( "When reduced, their affine bases has same dimension 2" ) {
483 CAPTURE( B.basis() );
484 CAPTURE( RB.basis() );
485 REQUIRE( B.dimension() == 2 );
486 REQUIRE( RB.dimension() == 2 );
487 }
488 THEN( "Their projections have the same geometry (i.e. orientations within points)" ) {
489 CAPTURE( lcm );
490 CAPTURE( rlcm );
491 CAPTURE( pX );
492 CAPTURE( pY );
493 REQUIRE( pX.size() == pY.size() );
494 // Computing arbitrary determinants between triplets of points
495 const std::size_t nb = 200;
496 std::size_t nb_ok = 0;
497 const double eps = 1e-12 * factor;
498 for ( auto i = 0; i < nb; i++ )
499 {
500 const std::size_t j = rand() % pX.size();
501 const std::size_t k = rand() % pX.size();
502 const std::size_t l = rand() % pX.size();
503 const auto u = pX[ k ] - pX[ j ];
504 const auto v = pX[ l ] - pX[ j ];
505 const auto ru = pY[ k ] - pY[ j ];
506 const auto rv = pY[ l ] - pY[ j ];
507 const auto det = determinant( u, v );//u [ 0 ] * v [ 1 ] - u [ 1 ] * v [ 0 ];
508 const auto rdet = determinant(ru,rv );//ru[ 0 ] * rv[ 1 ] - ru[ 1 ] * rv[ 0 ];
509 if ( rdet > eps ) nb_ok += ( det > 0 ) ? 1 : 0;
510 else if ( rdet < -eps ) nb_ok += ( det < 0 ) ? 1 : 0;
511 else nb_ok += ( det == 0 ) ? 1 : 0;
512 if ( nb_ok != i+1 )
513 {
514 std::cout << "eps=" << eps << " factor=" << factor << "\n";
515 std::cout << "u =" << u << " v =" << v << " det =" << det << "\n";
516 std::cout << "ru=" << ru << " rv=" << rv << " rdet=" << rdet << "\n";
517 break;
518 }
519 }
520 REQUIRE( nb_ok == nb );
521 }
522 }
523
524 GIVEN( "Given X a set of randomly generated points by adding linear combinations of 3 lattice vectors, and Y the same set but with real coordinates" ) {
525 typedef SpaceND< 3, int64_t > Space3;
526 typedef Space3::Point PPoint;
527 typedef Space3::RealPoint PRealPoint;
528 std::vector< Point > V = { Point{ 3, 4, 0, 2, 5 },
529 Point{ -2, -1, 5, -7, 1 },
530 Point{ -5, 2, 11, 1, 4 } };
531 auto X = makeRandomLatticePointsFromDirVectors( 20, V );
532 Basis B( X, Basis::Type::ECHELON_REDUCED );
533 std::vector< RealPoint > Y( X.size() );
534 for ( auto i = 0; i < Y.size(); i++ )
535 Basis::transform( Y[ i ], X[ i ] );
536 RealBasis RB( Y, RealBasis::Type::ECHELON_REDUCED );
537 std::vector< PPoint > pX;
538 std::vector< PRealPoint > pY;
539 auto lcm = B .projectPoints( pX, X );
540 auto rlcm = RB.projectPoints( pY, Y );
541 auto rbox = boundingBox( pY );
542 double factor = ( rbox.second - rbox.first ).squaredNorm();
543 THEN( "When reduced, their affine bases has same dimension 2" ) {
544 CAPTURE( B.basis() );
545 CAPTURE( RB.basis() );
546 REQUIRE( B.dimension() == 3 );
547 REQUIRE( RB.dimension() == 3 );
548 }
549 THEN( "Their projections have null remainder" ) {
550 CAPTURE( B.basis() );
551 unsigned int nb_null_rem = 0;
552 for ( auto i = 0; i < X.size(); i++ )
553 {
554 const auto [ d, L, r ] = B.decomposeVector( X[i] - B.origin() );
555 nb_null_rem += ( r == Point::zero ) ? 1 : 0;
556 if ( r != Point::zero )
557 std::cout << (X[i] - B.origin()) << " => " << d << "/" << L << "/" << r << "\n";
558 }
559 REQUIRE( nb_null_rem == X.size() );
560 }
561 THEN( "Their projections have the same geometry (i.e. orientations within points)" ) {
562 CAPTURE( B.basis() );
563 CAPTURE( RB.basis() );
564 CAPTURE( lcm );
565 CAPTURE( rlcm );
566 CAPTURE( pX );
567 CAPTURE( pY );
568 REQUIRE( pX.size() == pY.size() );
569 // Computing arbitrary determinants between quadruplets of points
570 const std::size_t nb = 200;
571 std::size_t nb_ok = 0;
572 const double eps = 1e-12 * factor;
573 for ( auto i = 0; i < nb; i++ )
574 {
575 const std::size_t j = rand() % pX.size();
576 const std::size_t k = rand() % pX.size();
577 const std::size_t l = rand() % pX.size();
578 const std::size_t m = rand() % pX.size();
579 const auto u = pX[ k ] - pX[ j ];
580 const auto v = pX[ l ] - pX[ j ];
581 const auto w = pX[ m ] - pX[ j ];
582 const auto ru = pY[ k ] - pY[ j ];
583 const auto rv = pY[ l ] - pY[ j ];
584 const auto rw = pY[ m ] - pY[ j ];
585 const auto det = determinant( u, v, w );
586 const auto rdet = determinant( ru, rv, rw );
587 if ( rdet > eps ) nb_ok += ( det > 0 ) ? 1 : 0;
588 else if ( rdet < -eps ) nb_ok += ( det < 0 ) ? 1 : 0;
589 else nb_ok += ( det == 0 ) ? 1 : 0;
590 if ( nb_ok != i+1 )
591 {
592 std::cout << "eps=" << eps << " factor=" << factor << "\n";
593 std::cout << "u =" << u << " v =" << v << " w =" << w
594 << " det =" << det << "\n";
595 std::cout << "ru=" << ru << " rv=" << rv << " rw=" << rw
596 << " rdet=" << rdet << "\n";
597 break;
598 }
599 }
600 REQUIRE( nb_ok == nb );
601 }
602 }
603
604}
605
606
607
608
DGtal::PointVector
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition PointVector.h:593
DGtal::functions::computeAffineDimension
DGtal::int64_t computeAffineDimension(const std::vector< TPoint > &X, const double tolerance=1e-12)
Definition AffineGeometry.h:1068
DGtal::functions::computeIndependentVector
TPoint computeIndependentVector(const std::vector< TPoint > &basis, const double tolerance=1e-12)
Definition AffineGeometry.h:1202
DGtal::functions::computeOrthogonalVectorToBasis
static TPoint computeOrthogonalVectorToBasis(const std::vector< TPoint > &basis)
Definition AffineGeometry.h:1262
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
Definition ClosedIntegerHalfPlane.h:49
std
STL namespace.
DGtal::AffineBasis
Aim: Utility class to determine the affine geometry of an input set of points. It provides exact resu...
Definition AffineBasis.h:112
DGtal::AffineGeometry
Aim: Utility class to determine the affine geometry of an input set of points. It provides exact resu...
Definition AffineGeometry.h:454
g
std::mt19937 g(rd())
determinant
T determinant(const PointVector< 2, T, C > &u, const PointVector< 2, T, C > &v)
Definition testAffineBasis.cpp:115
rd
std::random_device rd
Definition testAffineBasis.cpp:45
makeRandomLatticePointsFromDirVectors
std::vector< Point > makeRandomLatticePointsFromDirVectors(int nb, const vector< Point > &V)
Definition testAffineBasis.cpp:66
sameDirection
bool sameDirection(const TPoint &a, const TPoint &b)
Definition testAffineBasis.cpp:91
boundingBox
std::pair< TPoint, TPoint > boundingBox(const std::vector< TPoint > &v)
Definition testAffineBasis.cpp:98
makeRandomVectors
std::vector< Point > makeRandomVectors(int nb, int amplitude)
Definition testAffineBasis.cpp:50
Point
MyPointD Point
Definition testClone2.cpp:381
CAPTURE
CAPTURE(thicknessHV)
GIVEN
GIVEN("A cubical complex with random 3-cells")
Definition testCubicalComplex.cpp:70
REQUIRE
REQUIRE(domain.isInside(aPoint))
RealPoint
PointVector< 3, double > RealPoint
Definition testTriangulatedSurface.cpp:51
SCENARIO
SCENARIO("UnorderedSetByBlock< PointVector< 2, int > unit tests with 32 bits blocks", "[unorderedsetbyblock][2d]")
Definition testUnorderedSetByBlock.cpp:58