KiCad PCB EDA Suite
Loading...
Searching...
No Matches
test_kimath.cpp
Go to the documentation of this file.
1/*
2 * This program source code file is part of KiCad, a free EDA CAD application.
3 *
4 * Copyright (C) 2020 KiCad Developers, see AUTHORS.TXT for contributors.
5 *
6 * This program is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU General Public License
8 * as published by the Free Software Foundation; either version 2
9 * of the License, or (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, you may find one here:
18 * http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
19 * or you may search the http://www.gnu.org website for the version 2 license,
20 * or you may write to the Free Software Foundation, Inc.,
21 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
22 */
23
29
30// Code under test
31#include <trigo.h>
32#include <hash_128.h>
33#include <geometry/shape_arc.h>
34
35#include <iomanip>
36
37
38/*
39CircleCenterFrom3Points calculate the center of a circle defined by 3 points
40It is similar to CalcArcCenter( const VECTOR2D& aStart, const VECTOR2D& aMid, const VECTOR2D& aEnd )
41but it was needed to debug CalcArcCenter, so I keep it available for other issues in CalcArcCenter
42
43The perpendicular bisector of the segment between two points is the
44set of all points equidistant from both. So if you take the
45perpendicular bisector of (x1,y1) and (x2,y2) and the perpendicular
46bisector of the segment from (x2,y2) to (x3,y3) and find the
47intersection of those lines, that point will be the center.
48
49To find the equation of the perpendicular bisector of (x1,y1) to (x2,y2),
50you know that it passes through the midpoint of the segment:
51((x1+x2)/2,(y1+y2)/2), and if the slope of the line
52connecting (x1,y1) to (x2,y2) is m, the slope of the perpendicular
53bisector is -1/m. Work out the equations for the two lines, find
54their intersection, and bingo! You've got the coordinates of the center.
55
56An error should occur if the three points lie on a line, and you'll
57need special code to check for the case where one of the slopes is zero.
58
59see https://web.archive.org/web/20171223103555/http://mathforum.org/library/drmath/view/54323.html
60*/
61static bool Ref0CircleCenterFrom3Points( const VECTOR2D& p1, const VECTOR2D& p2, const VECTOR2D& p3,
62 VECTOR2D* aCenter )
63{
64 // Move coordinate origin to p2, to simplify calculations
65 VECTOR2D b = p1 - p2;
66 VECTOR2D d = p3 - p2;
67 double bc = ( b.x * b.x + b.y * b.y ) / 2.0;
68 double cd = ( -d.x * d.x - d.y * d.y ) / 2.0;
69 double det = -b.x * d.y + d.x * b.y;
70
71 if( fabs( det ) < 1.0e-6 ) // arbitrary limit to avoid divide by 0
72 return false;
73
74 det = 1 / det;
75 aCenter->x = ( -bc * d.y - cd * b.y ) * det;
76 aCenter->y = ( b.x * cd + d.x * bc ) * det;
77 *aCenter += p2;
78
79 return true;
80}
81
82
83// Based on https://paulbourke.net/geometry/circlesphere/
84// Uses the Y'b equation too, depending on slope values
85static const VECTOR2D Ref1CalcArcCenter( const VECTOR2D& aStart, const VECTOR2D& aMid, const VECTOR2D& aEnd )
86{
87 VECTOR2D center;
88 double yDelta_21 = aMid.y - aStart.y;
89 double xDelta_21 = aMid.x - aStart.x;
90 double yDelta_32 = aEnd.y - aMid.y;
91 double xDelta_32 = aEnd.x - aMid.x;
92
93 // This is a special case for aMid as the half-way point when aSlope = 0 and bSlope = inf
94 // or the other way around. In that case, the center lies in a straight line between
95 // aStart and aEnd
96 if( ( ( xDelta_21 == 0.0 ) && ( yDelta_32 == 0.0 ) )
97 || ( ( yDelta_21 == 0.0 ) && ( xDelta_32 == 0.0 ) ) )
98 {
99 center.x = ( aStart.x + aEnd.x ) / 2.0;
100 center.y = ( aStart.y + aEnd.y ) / 2.0;
101 return center;
102 }
103
104 // Prevent div=0 errors
105 /*if( xDelta_21 == 0.0 )
106 xDelta_21 = std::numeric_limits<double>::epsilon();
107
108 if( xDelta_32 == 0.0 )
109 xDelta_32 = -std::numeric_limits<double>::epsilon();*/
110
111 double aSlope = yDelta_21 / xDelta_21;
112 double bSlope = yDelta_32 / xDelta_32;
113
114 if( aSlope == bSlope )
115 {
116 if( aStart == aEnd )
117 {
118 // This is a special case for a 360 degrees arc. In this case, the center is halfway between
119 // the midpoint and either newEnd point
120 center.x = ( aStart.x + aMid.x ) / 2.0;
121 center.y = ( aStart.y + aMid.y ) / 2.0;
122 return center;
123 }
124 else
125 {
126 // If the points are colinear, the center is at infinity, so offset
127 // the slope by a minimal amount
128 // Warning: This will induce a small error in the center location
129 /*aSlope += std::numeric_limits<double>::epsilon();
130 bSlope -= std::numeric_limits<double>::epsilon();*/
131 }
132 }
133
134 center.x = ( aSlope * bSlope * ( aStart.y - aEnd.y ) + bSlope * ( aStart.x + aMid.x )
135 - aSlope * ( aMid.x + aEnd.x ) )
136 / ( 2 * ( bSlope - aSlope ) );
137
138 if( std::abs( aSlope ) > std::abs( bSlope ) )
139 {
140 center.y = ( ( ( aStart.x + aMid.x ) / 2.0 - center.x ) / aSlope
141 + ( aStart.y + aMid.y ) / 2.0 );
142 }
143 else
144 {
145 center.y =
146 ( ( ( aMid.x + aEnd.x ) / 2.0 - center.x ) / bSlope + ( aMid.y + aEnd.y ) / 2.0 );
147 }
148
149 return center;
150}
151
152
154{
158};
159
160
164BOOST_AUTO_TEST_SUITE( KiMath )
165
166
167BOOST_AUTO_TEST_CASE( TestCalcArcCenter3Pts )
168{
169 // Order: start, mid, end
170 std::vector<TEST_CALC_ARC_CENTER_CASE> calc_center_cases = {
171 //{ { 183000000, 89000000 }, { 185333332, 89000000 }, { 185333333, 91333333 } } // Fails currently
172 };
173
174 const double tolerance = 1.0;
175
176 // Check random points (fails currently)
177 /*for( int i = 0; i < 100; i++ )
178 {
179 TEST_CALC_ARC_CENTER_CASE cas;
180
181 cas.istart.x = rand();
182 cas.istart.y = rand();
183
184 cas.imid.x = rand();
185 cas.imid.y = rand();
186
187 cas.iend.x = rand();
188 cas.iend.y = rand();
189
190 calc_center_cases.push_back( cas );
191 }*/
192
193 for( const TEST_CALC_ARC_CENTER_CASE& entry : calc_center_cases )
194 {
195 wxString msg;
196
197 VECTOR2D start( entry.istart );
198 VECTOR2D mid( entry.imid );
199 VECTOR2D end( entry.iend );
200
201 VECTOR2D calcCenter = CalcArcCenter( start, mid, end );
202
203 double crs = ( calcCenter - start ).EuclideanNorm();
204 double crm = ( calcCenter - mid ).EuclideanNorm();
205 double cre = ( calcCenter - end ).EuclideanNorm();
206
207 double cavg = ( crs + crm + cre ) / 3.0;
208
209 if( std::abs( crs - cavg ) > tolerance || std::abs( crm - cavg ) > tolerance
210 || std::abs( cre - cavg ) > tolerance )
211 {
212 msg << "CalcArcCenter failed.";
213
214 msg << "\nstart: " << entry.istart.Format();
215 msg << "\nmid: " << entry.imid.Format();
216 msg << "\nend: " << entry.iend.Format();
217
218 {
219 msg << "\nCalculated center: " << wxString::Format( "%.15f", calcCenter.x ) << ", "
220 << wxString::Format( "%.15f", calcCenter.y );
221
222 msg << "\n Avg radius: " << wxString::Format( "%.15f", cavg );
223 msg << "\nStart radius: " << wxString::Format( "%.15f", crs );
224 msg << "\n Mid radius: " << wxString::Format( "%.15f", crm );
225 msg << "\n End radius: " << wxString::Format( "%.15f", cre );
226 msg << "\n";
227
228 // Check mid/end points using the calculated center (like SHAPE_ARC)
229 EDA_ANGLE angStart( start - calcCenter );
230 EDA_ANGLE angMid( mid - calcCenter );
231 EDA_ANGLE angEnd( end - calcCenter );
232
233 EDA_ANGLE angCenter = angEnd - angStart;
234
235 VECTOR2D newMid = start;
236 VECTOR2D newEnd = start;
237
238 RotatePoint( newMid, calcCenter, -angCenter / 2.0 );
239 RotatePoint( newEnd, calcCenter, -angCenter );
240
241 msg << "\nNew mid: " << wxString::Format( "%.15f", newMid.x ) << ", "
242 << wxString::Format( "%.15f", newMid.y );
243
244 msg << "\nNew end: " << wxString::Format( "%.15f", newEnd.x ) << ", "
245 << wxString::Format( "%.15f", newEnd.y );
246 msg << "\n";
247
248 double endsDist = ( newEnd - end ).EuclideanNorm();
249
250 msg << "\nNew end is off by " << wxString::Format( "%.15f", endsDist );
251 msg << "\n";
252 }
253
254 {
255 VECTOR2D ref0Center;
256 Ref0CircleCenterFrom3Points( start, mid, end, &ref0Center );
257
258 double r0_rs = ( ref0Center - start ).EuclideanNorm();
259 double r0_rm = ( ref0Center - mid ).EuclideanNorm();
260 double r0_rre = ( ref0Center - end ).EuclideanNorm();
261
262 double r0_ravg = ( r0_rs + r0_rm + r0_rre ) / 3.0;
263
264 msg << "\nReference0 center: " << wxString::Format( "%.15f", ref0Center.x ) << ", "
265 << wxString::Format( "%.15f", ref0Center.y );
266
267 msg << "\nRef0 Avg radius: " << wxString::Format( "%.15f", r0_ravg );
268 msg << "\nRef0 Start radius: " << wxString::Format( "%.15f", r0_rs );
269 msg << "\nRef0 Mid radius: " << wxString::Format( "%.15f", r0_rm );
270 msg << "\nRef0 End radius: " << wxString::Format( "%.15f", r0_rre );
271 msg << "\n";
272 }
273
274 {
275 VECTOR2D ref1Center = Ref1CalcArcCenter( start, mid, end );
276
277 double r1_rs = ( ref1Center - start ).EuclideanNorm();
278 double r1_rm = ( ref1Center - mid ).EuclideanNorm();
279 double r1_rre = ( ref1Center - end ).EuclideanNorm();
280
281 double r1_ravg = ( r1_rs + r1_rm + r1_rre ) / 3.0;
282
283 msg << "\nReference1 center: " << wxString::Format( "%.15f", ref1Center.x ) << ", "
284 << wxString::Format( "%.15f", ref1Center.y );
285
286 msg << "\nRef1 Avg radius: " << wxString::Format( "%.15f", r1_ravg );
287 msg << "\nRef1 Start radius: " << wxString::Format( "%.15f", r1_rs );
288 msg << "\nRef1 Mid radius: " << wxString::Format( "%.15f", r1_rm );
289 msg << "\nRef1 End radius: " << wxString::Format( "%.15f", r1_rre );
290 msg << "\n";
291 }
292
293
294 BOOST_CHECK_MESSAGE( false, msg );
295 }
296 }
297}
298
299
300BOOST_AUTO_TEST_CASE( TestInterceptsPositiveX )
301{
302 BOOST_CHECK( !InterceptsPositiveX( 10.0, 20.0 ) );
303 BOOST_CHECK( !InterceptsPositiveX( 10.0, 120.0 ) );
304 BOOST_CHECK( !InterceptsPositiveX( 10.0, 220.0 ) );
305 BOOST_CHECK( !InterceptsPositiveX( 10.0, 320.0 ) );
306 BOOST_CHECK( InterceptsPositiveX( 20.0, 10.0 ) );
307 BOOST_CHECK( InterceptsPositiveX( 345.0, 15.0 ) );
308}
309
310
311BOOST_AUTO_TEST_CASE( TestInterceptsNegativeX )
312{
313 BOOST_CHECK( !InterceptsNegativeX( 10.0, 20.0 ) );
314 BOOST_CHECK( !InterceptsNegativeX( 10.0, 120.0 ) );
315 BOOST_CHECK( InterceptsNegativeX( 10.0, 220.0 ) );
316 BOOST_CHECK( InterceptsNegativeX( 10.0, 320.0 ) );
317 BOOST_CHECK( InterceptsNegativeX( 20.0, 10.0 ) );
318 BOOST_CHECK( !InterceptsNegativeX( 345.0, 15.0 ) );
319 BOOST_CHECK( InterceptsNegativeX( 145.0, 225.0 ) );
320}
321
322
324{
325 const HASH_128 zero{};
326
327 for( size_t i = 0; i < 16; i++ )
328 BOOST_CHECK( zero.Value8[i] == 0 );
329
330 HASH_128* zero_h = new HASH_128;
331
332 for( size_t i = 0; i < 16; i++ )
333 BOOST_CHECK( zero_h->Value8[i] == 0 );
334
335 delete zero_h;
336
337 HASH_128 h;
338 h.Value64[0] = 0x00CDEF0123456789ULL;
339 h.Value64[1] = 0x56789ABCDEF01234ULL;
340
341 BOOST_CHECK( h != zero );
342
343 HASH_128 b = h;
344 BOOST_CHECK( b != zero );
345
346 BOOST_CHECK( b == h );
347 BOOST_CHECK( h == b );
348 BOOST_CHECK_EQUAL( h.ToString(), std::string( "00CDEF012345678956789ABCDEF01234" ) );
349}
350
351
352//-----------------------------------------------------------------------------
353// Test MMH3_HASH against the original MurmurHash3_x64_128 implementation
354
355#include <mmh3_hash.h>
356
357void MurmurHash3_x64_128( const void* key, const int len, const uint32_t seed, void* out );
358
359BOOST_AUTO_TEST_CASE( TestMMH3Hash )
360{
361 std::vector<std::vector<int32_t>> data;
362
363 for( size_t i = 0; i < 10; i++ )
364 {
365 std::vector<int32_t>& vec = data.emplace_back();
366
367 size_t vecSize = rand() % 128;
368
369 for( size_t j = 0; j < vecSize; j++ )
370 vec.emplace_back( rand() );
371 }
372
373 int64_t seed = 0; // Test 0 seed first
374
375 for( const std::vector<int32_t>& vec : data )
376 {
377 MMH3_HASH mmh3( seed );
378
379 for( const int32_t val : vec )
380 mmh3.add( val );
381
382 HASH_128 h128 = mmh3.digest();
383
384 HASH_128 orig128;
385 MurmurHash3_x64_128( (void*) vec.data(), vec.size() * sizeof( vec[0] ), seed,
386 (void*) orig128.Value64 );
387
388 BOOST_CHECK( h128 == orig128 );
389
390 // Generate new random seed
391 seed = rand();
392 }
393}
394
395
396//-----------------------------------------------------------------------------
397// The original MurmurHash3_x64_128 implementation
398
399
400// Microsoft Visual Studio
401#if defined( _MSC_VER )
402 #define FORCE_INLINE __forceinline
403 #include <stdlib.h>
404 #define ROTL64( x, y ) _rotl64( x, y )
405 #define BIG_CONSTANT( x ) ( x )
406 // Other compilers
407#else // defined(_MSC_VER)
408 #define FORCE_INLINE inline __attribute__( ( always_inline ) )
409 inline uint64_t rotl64( uint64_t x, int8_t r )
410 {
411 return ( x << r ) | ( x >> ( 64 - r ) );
412 }
413 #define ROTL64( x, y ) rotl64( x, y )
414 #define BIG_CONSTANT( x ) ( x##LLU )
415#endif // !defined(_MSC_VER)
416
417
418FORCE_INLINE uint64_t getblock64 ( const uint64_t * p, int i )
419{
420 return p[i];
421}
422
423FORCE_INLINE uint64_t fmix64 ( uint64_t k )
424{
425 k ^= k >> 33;
426 k *= BIG_CONSTANT(0xff51afd7ed558ccd);
427 k ^= k >> 33;
428 k *= BIG_CONSTANT(0xc4ceb9fe1a85ec53);
429 k ^= k >> 33;
430
431 return k;
432}
433
434void MurmurHash3_x64_128 ( const void * key, const int len,
435 const uint32_t seed, void * out )
436{
437 const uint8_t * data = (const uint8_t*)key;
438 const int nblocks = len / 16;
439
440 uint64_t h1 = seed;
441 uint64_t h2 = seed;
442
443 const uint64_t c1 = BIG_CONSTANT(0x87c37b91114253d5);
444 const uint64_t c2 = BIG_CONSTANT(0x4cf5ad432745937f);
445
446 //----------
447 // body
448
449 const uint64_t * blocks = (const uint64_t *)(data);
450
451 for(int i = 0; i < nblocks; i++)
452 {
453 uint64_t k1 = getblock64(blocks,i*2+0);
454 uint64_t k2 = getblock64(blocks,i*2+1);
455
456 k1 *= c1; k1 = ROTL64(k1,31); k1 *= c2; h1 ^= k1;
457
458 h1 = ROTL64(h1,27); h1 += h2; h1 = h1*5+0x52dce729;
459
460 k2 *= c2; k2 = ROTL64(k2,33); k2 *= c1; h2 ^= k2;
461
462 h2 = ROTL64(h2,31); h2 += h1; h2 = h2*5+0x38495ab5;
463 }
464
465 //----------
466 // tail
467
468 const uint8_t * tail = (const uint8_t*)(data + nblocks*16);
469
470 uint64_t k1 = 0;
471 uint64_t k2 = 0;
472
473 switch(len & 15)
474 {
475 case 15: k2 ^= ((uint64_t)tail[14]) << 48;
476 case 14: k2 ^= ((uint64_t)tail[13]) << 40;
477 case 13: k2 ^= ((uint64_t)tail[12]) << 32;
478 case 12: k2 ^= ((uint64_t)tail[11]) << 24;
479 case 11: k2 ^= ((uint64_t)tail[10]) << 16;
480 case 10: k2 ^= ((uint64_t)tail[ 9]) << 8;
481 case 9: k2 ^= ((uint64_t)tail[ 8]) << 0;
482 k2 *= c2; k2 = ROTL64(k2,33); k2 *= c1; h2 ^= k2;
483
484 case 8: k1 ^= ((uint64_t)tail[ 7]) << 56;
485 case 7: k1 ^= ((uint64_t)tail[ 6]) << 48;
486 case 6: k1 ^= ((uint64_t)tail[ 5]) << 40;
487 case 5: k1 ^= ((uint64_t)tail[ 4]) << 32;
488 case 4: k1 ^= ((uint64_t)tail[ 3]) << 24;
489 case 3: k1 ^= ((uint64_t)tail[ 2]) << 16;
490 case 2: k1 ^= ((uint64_t)tail[ 1]) << 8;
491 case 1: k1 ^= ((uint64_t)tail[ 0]) << 0;
492 k1 *= c1; k1 = ROTL64(k1,31); k1 *= c2; h1 ^= k1;
493 };
494
495 //----------
496 // finalization
497
498 h1 ^= len; h2 ^= len;
499
500 h1 += h2;
501 h2 += h1;
502
503 h1 = fmix64(h1);
504 h2 = fmix64(h2);
505
506 h1 += h2;
507 h2 += h1;
508
509 ((uint64_t*)out)[0] = h1;
510 ((uint64_t*)out)[1] = h2;
511}
512
513//-----------------------------------------------------------------------------
514
515#undef FORCE_INLINE
516#undef ROTL64
517#undef BIG_CONSTANT
518
519
A streaming C++ equivalent for MurmurHash3_x64_128.
Definition: mmh3_hash.h:60
FORCE_INLINE void add(const std::string &input)
Definition: mmh3_hash.h:95
FORCE_INLINE HASH_128 digest()
Definition: mmh3_hash.h:114
EDA_ANGLE abs(const EDA_ANGLE &aAngle)
Definition: eda_angle.h:390
A storage class for 128-bit hash value.
Definition: hash_128.h:36
uint8_t Value8[16]
Definition: hash_128.h:59
std::string ToString() const
Definition: hash_128.h:47
uint64_t Value64[2]
Definition: hash_128.h:61
BOOST_AUTO_TEST_SUITE(CadstarPartParser)
BOOST_AUTO_TEST_SUITE_END()
FORCE_INLINE uint64_t fmix64(uint64_t k)
uint64_t rotl64(uint64_t x, int8_t r)
static const VECTOR2D Ref1CalcArcCenter(const VECTOR2D &aStart, const VECTOR2D &aMid, const VECTOR2D &aEnd)
Definition: test_kimath.cpp:85
FORCE_INLINE uint64_t getblock64(const uint64_t *p, int i)
#define ROTL64(x, y)
void MurmurHash3_x64_128(const void *key, const int len, const uint32_t seed, void *out)
static bool Ref0CircleCenterFrom3Points(const VECTOR2D &p1, const VECTOR2D &p2, const VECTOR2D &p3, VECTOR2D *aCenter)
Test suite for KiCad math code.
Definition: test_kimath.cpp:61
#define BIG_CONSTANT(x)
BOOST_AUTO_TEST_CASE(TestCalcArcCenter3Pts)
Declare the test suite.
#define FORCE_INLINE
void RotatePoint(int *pX, int *pY, const EDA_ANGLE &aAngle)
Calculate the new point of coord coord pX, pY, for a rotation center 0, 0.
Definition: trigo.cpp:229
bool InterceptsPositiveX(double aStartAngle, double aEndAngle)
Test if an arc from aStartAngle to aEndAngle crosses the positive X axis (0 degrees).
Definition: trigo.h:215
bool InterceptsNegativeX(double aStartAngle, double aEndAngle)
Test if an arc from aStartAngle to aEndAngle crosses the negative X axis (180 degrees).
Definition: trigo.h:233
const VECTOR2I CalcArcCenter(const VECTOR2I &aStart, const VECTOR2I &aMid, const VECTOR2I &aEnd)
Determine the center of an arc or circle given three points on its circumference.
Definition: trigo.cpp:521