1%% options 2copyright owner = Dirk Krause 3copyright year = 2017-xxxx 4SPDX-License-Identifier: BSD-3-Clause 5 6%% header 7 8#include <libdk4base/dk4types.h> 9 10#ifdef __cplusplus 11extern "C" { 12#endif 13 14/** Arc calculation. 15 For 3 given points (x1;y1), (x2;y2) and (x3;y3) calculate 16 center point, radius, start and end angles and direction. 17 If there is no solution or an infinite number of solutions, 18 set the direction to 0. In this case drawing functions will draw 19 two lines instead of an arc. 20 @param pxc Address of result variable for center point x. 21 @param pyc Address of result variable for center point y. 22 @param pr Address of result variable for radius. 23 @param pa Address of result variable for start angle in radians 24 (-pi to pi). 25 @param pb Address of result variable for end angle in radians 26 (-5*pi to 5*pi). 27 @param pd Address of result variable for direction. 28 @param x1 Point 1 x. 29 @param y1 Point 1 y. 30 @param x2 Point 2 x. 31 @param y2 Point 2 y. 32 @param x3 Point 3 x. 33 @param y3 Point 3 y. 34*/ 35void 36wxdarc_calculation( 37 double *pxc, 38 double *pyc, 39 double *pr, 40 double *pa, 41 double *pb, 42 int8_t *pd, 43 int32_t x1, 44 int32_t y1, 45 int32_t x2, 46 int32_t y2, 47 int32_t x3, 48 int32_t y3 49); 50 51#ifdef __cplusplus 52} 53#endif 54 55 56%% module 57 58#include "dk4conf.h" 59 60#if DK4_HAVE_MATH_H 61#ifndef MATH_H_INCLUDED 62#define _USE_MATH_DEFINES 63#include <math.h> 64#define MATH_H_INCLUDED 1 65#endif 66#endif 67 68#ifndef DK4ERROR_H_INCLUDED 69#include <libdk4base/dk4error.h> 70#endif 71 72#ifndef DK4MATH_H_INCLUDED 73#include <libdk4c/dk4math.h> 74#endif 75 76#ifndef DK4MAADI_H_INCLUDED 77#include <libdk4ma/dk4maadi.h> 78#endif 79 80#ifndef WXDARC_H_INCLUDED 81#include <wxdkdraw/wxdarc.h> 82#endif 83 84 85 86$!trace-include 87 88 89/* 90 The solution was found using Maxima: 91 92 eq1: r^2=(x1-xc)^2+(y1-yc)^2; 93 eq2: r^2=(x2-xc)^2+(y2-yc)^2; 94 eq3: r^2=(x3-xc)^2+(y3-yc)^2; 95 solve([eq1,eq2,eq3],[xc,yc,r]); 96 97 From the two solution triples, we use the one with non-negative radius. 98 99 The calculation can have a different number of solutions: 100 101 - No solution 102 if all points are placed on a line. 103 104 - One solution 105 for normal arcs specified by 3 different points. 106 107 - An infinite number of solutions 108 if two or all 3 points are equal. 109 110 The special cases (no solution, infinite number of solutions) will 111 result in n=0 in the function below, thus in inifinite xc and yc. 112 113 In the direction calculation we use two vectors (first vector from 114 point 1 to point 2, second vector from point 2 to point 3). 115 The z component of the vector product of these vectors is 0 if 116 the 3 points are on a line or if 2 or more points are equal. 117 The z component is positive if the rotation of the arc is positive 118 (counterclockwise). 119 For performance reasons the vector product is calculated in integer 120 arithmetics by default. If integer arithmetics results in overflow, 121 the calculation uses double values. 122 123*/ 124 125 126void 127wxdarc_calculation( 128 double *pxc, 129 double *pyc, 130 double *pr, 131 double *pa, 132 double *pb, 133 int8_t *pd, 134 int32_t x1, 135 int32_t y1, 136 int32_t x2, 137 int32_t y2, 138 int32_t x3, 139 int32_t y3 140) 141{ 142 dk4_er_t er; /* Error report to recognize overflow */ 143 double zx; /* Divident to calculate xc */ 144 double zy; /* Divident to calculate yc */ 145 double n; /* Divisor to calculate xc and yc */ 146 double dx1; /* Point 1 X as double */ 147 double dy1; /* Point 1 Y as double */ 148 double dx2; /* Point 2 X as double */ 149 double dy2; /* Point 2 Y as double */ 150 double dx3; /* Point 3 X as double */ 151 double dy3; /* Point 3 Y as double */ 152 double sx1; /* Square of point 1 X */ 153 double sy1; /* Square of point 1 Y */ 154 double sx2; /* Square of point 2 X */ 155 double sy2; /* Square of point 2 Y */ 156 double sx3; /* Square of point 3 X */ 157 double sy3; /* Square of point 3 Y */ 158 double xc = 0.0; /* Center point x */ 159 double yc = 0.0; /* Center point y */ 160 double r = 0.0; /* Radius */ 161 double a = 0.0; /* Angle in radians, point 1 */ 162 double b = 0.0; /* Angle in radians, point 3 */ 163 double g = 0.0; /* Angle in radians, point 2 */ 164 double dvp; /* Z component of vector product */ 165 dk4_im_t deltax1; /* X component of first vector */ 166 dk4_im_t deltay1; /* Y component of first vector */ 167 dk4_im_t deltax2; /* X component of second vector */ 168 dk4_im_t deltay2; /* Y component of second vector */ 169 dk4_im_t vp; /* Z component of vector product */ 170 int8_t d = (int8_t)0; /* Positive or negative arc direction */ 171 $? "+ wxdarc_calculation" 172 173 /* Show arguments in debugging 174 */ 175 $? ". x1 = %d", (int)x1 176 $? ". y1 = %d", (int)y1 177 $? ". x2 = %d", (int)x2 178 $? ". y2 = %d", (int)y2 179 $? ". x3 = %d", (int)x3 180 $? ". y3 = %d", (int)y3 181 182 /* Initialize result variables 183 */ 184#if 0 185 /* 2020-03-24 Modification 186 Turned into initialization. 187 */ 188 xc = 0.0; 189 yc = 0.0; 190 r = 0.0; 191 a = 0.0; 192 b = 0.0; 193 g = 0.0; 194 d = (int8_t)0; 195#endif 196 197 /* Convert integers to doubles 198 */ 199 dx1 = (double)x1; dx2 = (double)x2; dx3 = (double)x3; 200 dy1 = (double)y1; dy2 = (double)y2; dy3 = (double)y3; 201 202 /* Find direction 203 Use integer arithmetics by default, attempt with 204 double arithmetics if integer arithmetics results in overflow. 205 */ 206 dk4error_init(&er); 207 deltax1 = dk4ma_im_sub((dk4_im_t)x2, (dk4_im_t)x1, &er); 208 deltay1 = dk4ma_im_sub((dk4_im_t)y2, (dk4_im_t)y1, &er); 209 deltax2 = dk4ma_im_sub((dk4_im_t)x3, (dk4_im_t)x2, &er); 210 deltay2 = dk4ma_im_sub((dk4_im_t)y3, (dk4_im_t)y2, &er); 211 vp = dk4ma_im_sub( 212 dk4ma_im_mul(deltax1, deltay2, &er), 213 dk4ma_im_mul(deltax2, deltay1, &er), 214 &er 215 ); 216 if (DK4_E_NONE == er.ec) { 217 /* 218 No overflow occured, can use values 219 */ 220 if ((dk4_im_t)0L < vp) { 221 d = (int8_t)1; 222 } 223 else { 224 if ((dk4_im_t)0L > vp) { 225 d = (int8_t)(-1); 226 } 227 } 228 } 229 else { 230 /* On integer overflow attempt calculation in double 231 */ 232 if (DK4_E_MATH_OVERFLOW == er.ec) { 233 dvp = (dx2 - dx1) * (dy3 - dy2) - (dx3 - dx2) * (dy2 - dy1); 234 if (0 != dk4ma_is_finite(dvp)) { 235 if (1.0e-8 < fabs(dvp)) { 236 if (0.0 < dvp) { 237 d = (int8_t)1; 238 } 239 else { 240 d = (int8_t)(-1); 241 } 242 } 243 } 244 } 245 } 246 247 /* Attempt further calculations only if d is not 0 248 */ 249 if ((int8_t)0 != d) { 250 /* 251 Calculate squares values 252 */ 253 sx1 = dx1 * dx1; sx2 = dx2 * dx2; sx3 = dx3 * dx3; 254 sy1 = dy1 * dy1; sy2 = dy2 * dy2; sy3 = dy3 * dy3; 255 /* 256 Calculate determinants 257 */ 258 n = 2.0 * ((dx2 - dx1) * dy3 + (dx1 - dx3) * dy2 + (dx3 - dx2) * dy1); 259 zx = (dy2-dy1)*sy3 + (sx1-sx2+sy1-sy2)*dy3 260 + dy1*sy2 + (sx3-sx1-sy1)*dy2 + (sx2-sx3)*dy1; 261 zx *= -1.0; 262 zy = (dx2-dx1)*sy3 + (dx1-dx3)*sy2 + (dx3-dx2)*sy1 + (dx2-dx1)*sx3 263 + (sx1-sx2)*dx3 + dx1*sx2 - sx1*dx2; 264 /* 265 Calculate center point 266 */ 267 xc = zx / n; 268 yc = zy / n; 269 /* 270 Radius and angles 271 */ 272 r = sqrt((dx1-xc)*(dx1-xc)+(dy1-yc)*(dy1-yc)); 273 a = atan2((dy1-yc), (dx1-xc)); 274 b = atan2((dy3-yc), (dx3-xc)); 275 g = atan2((dy2-yc), (dx2-xc)); 276 /* 277 Check for infinite and undefined values 278 (no solution or multiple solutions will result in division by 0) 279 */ 280 if (0 == dk4ma_is_finite(xc)) { d = (int8_t)0; } 281 if (0 == dk4ma_is_finite(yc)) { d = (int8_t)0; } 282 if (0 == dk4ma_is_finite(r)) { d = (int8_t)0; } 283 if (0 == dk4ma_is_finite(a)) { d = (int8_t)0; } 284 if (0 == dk4ma_is_finite(b)) { d = (int8_t)0; } 285 if (0 == dk4ma_is_finite(g)) { d = (int8_t)0; } 286 /* 287 Angle corrections depending on arc direction 288 */ 289 if ((int8_t)0 < d) { 290 while (g < a) { g += (2.0 * M_PI); } 291 while (b < g) { b += (2.0 * M_PI); } 292 } 293 else { 294 if ((int8_t)0 > d) { 295 while (g > a) { g -= (2.0 * M_PI); } 296 while (b > g) { b -= (2.0 * M_PI); } 297 } 298 else { 299 xc = 0.0; yc = 0.0; r = 0.0; a = 0.0; b = 0.0; 300 } 301 } 302 } 303 304 /* Transfer results to destination variables 305 */ 306 if (NULL != pxc) { *pxc = xc; } 307 if (NULL != pyc) { *pyc = yc; } 308 if (NULL != pr) { *pr = r; } 309 if (NULL != pa) { *pa = a; } 310 if (NULL != pb) { *pb = b; } 311 if (NULL != pd) { *pd = d; } 312 313 /* Show results in debugging 314 */ 315 $? ". xc = %g", xc 316 $? ". yc = %g", yc 317 $? ". r = %g", r 318 $? ". a = %g", a 319 $? ". b = %g", b 320 $? ". g = %g", g 321 $? "- wxdarc_calculation" 322} 323 324 325 326#if TEST_WXDARC 327 328/* To build test program: 329 gcc -DTEST_WXDARC=1 -o wxdarc -I . wxdarc.c -ldk4c -ldk4ma -ldk4base -lm 330*/ 331 332int 333main(int argc, char *argv[]) 334{ 335 int iv[6]; 336 double xc = 1024.0; 337 double yc = 1024.0; 338 double r = -0.5; 339 double a = 27.0; 340 double b = 128.0; 341 int8_t d = -100; 342 int32_t x1 = 1024L; 343 int32_t y1 = 0L; 344 int32_t x2 = 724L; 345 int32_t y2 = 724L; 346 int32_t x3 = 0L; 347 int32_t y3 = 1024L; 348 int ok = 1; 349 int i; 350 351 if (7 == argc) { 352 for (i = 0; i < 6; i++) { 353 if (0 == sscanf(argv[i + 1], "%d", &(iv[i]))) { 354 ok = 0; 355 } 356 } 357 if (0 != ok) { 358 x1 = (int32_t)(iv[0]); 359 y1 = (int32_t)(iv[1]); 360 x2 = (int32_t)(iv[2]); 361 y2 = (int32_t)(iv[3]); 362 x3 = (int32_t)(iv[4]); 363 y3 = (int32_t)(iv[5]); 364 } 365 } 366 367 wxdarc_calculation( 368 &xc, &yc, &r, &a, &b, &d, 369 x1, y1, x2, y2, x3, y3 370 ); 371 372 printf("x1 = %d\n", (int)x1); 373 printf("y1 = %d\n", (int)y1); 374 printf("x2 = %d\n", (int)x2); 375 printf("y2 = %d\n", (int)y2); 376 printf("x3 = %d\n", (int)x3); 377 printf("y3 = %d\n", (int)y3); 378 printf("x = %g\n", xc); 379 printf("y = %g\n", yc); 380 printf("r = %g\n", r); 381 printf("a = %g\n", a); 382 printf("b = %g\n", b); 383 printf("d = %d\n", (int)d); 384 385 return 0; 386} 387 388 389#endif 390 391 392 393/* vim: set ai sw=4 ts=4 : */ 394