%% options copyright owner = Dirk Krause copyright year = 2017-xxxx SPDX-License-Identifier: BSD-3-Clause %% header #include "dk4types.h" #ifdef __cplusplus extern "C" { #endif /** Arc calculation. For 3 given points (x1;y1), (x2;y2) and (x3;y3) calculate center point, radius, start and end angles and direction. If there is no solution or an infinite number of solutions, set the direction to 0. In this case drawing functions will draw two lines instead of an arc. @param pxc Address of result variable for center point x. @param pyc Address of result variable for center point y. @param pr Address of result variable for radius. @param pa Address of result variable for start angle in radians (-pi to pi). @param pb Address of result variable for end angle in radians (-5*pi to 5*pi). @param pd Address of result variable for direction. @param x1 Point 1 x. @param y1 Point 1 y. @param x2 Point 2 x. @param y2 Point 2 y. @param x3 Point 3 x. @param y3 Point 3 y. */ void wxdarc_calculation( double *pxc, double *pyc, double *pr, double *pa, double *pb, int8_t *pd, int32_t x1, int32_t y1, int32_t x2, int32_t y2, int32_t x3, int32_t y3 ); #ifdef __cplusplus } #endif %% module #include "dk4conf.h" #if DK4_HAVE_MATH_H #ifndef MATH_H_INCLUDED #define _USE_MATH_DEFINES #include #define MATH_H_INCLUDED 1 #endif #endif #ifndef DK4ERROR_H_INCLUDED #include "dk4error.h" #endif #ifndef DK4MATH_H_INCLUDED #include "dk4math.h" #endif #ifndef DK4MAADI_H_INCLUDED #include "dk4maadi.h" #endif #ifndef WXDARC_H_INCLUDED #include "wxdarc.h" #endif $!trace-include /* The solution was found using Maxima: eq1: r^2=(x1-xc)^2+(y1-yc)^2; eq2: r^2=(x2-xc)^2+(y2-yc)^2; eq3: r^2=(x3-xc)^2+(y3-yc)^2; solve([eq1,eq2,eq3],[xc,yc,r]); From the two solution triples, we use the one with non-negative radius. The calculation can have a different number of solutions: - No solution if all points are placed on a line. - One solution for normal arcs specified by 3 different points. - An infinite number of solutions if two or all 3 points are equal. The special cases (no solution, infinite number of solutions) will result in n=0 in the function below, thus in inifinite xc and yc. In the direction calculation we use two vectors (first vector from point 1 to point 2, second vector from point 2 to point 3). The z component of the vector product of these vectors is 0 if the 3 points are on a line or if 2 or more points are equal. The z component is positive if the rotation of the arc is positive (counterclockwise). For performance reasons the vector product is calculated in integer arithmetics by default. If integer arithmetics results in overflow, the calculation uses double values. */ void wxdarc_calculation( double *pxc, double *pyc, double *pr, double *pa, double *pb, int8_t *pd, int32_t x1, int32_t y1, int32_t x2, int32_t y2, int32_t x3, int32_t y3 ) { dk4_er_t er; /* Error report to recognize overflow */ double zx; /* Divident to calculate xc */ double zy; /* Divident to calculate yc */ double n; /* Divisor to calculate xc and yc */ double dx1; /* Point 1 X as double */ double dy1; /* Point 1 Y as double */ double dx2; /* Point 2 X as double */ double dy2; /* Point 2 Y as double */ double dx3; /* Point 3 X as double */ double dy3; /* Point 3 Y as double */ double sx1; /* Square of point 1 X */ double sy1; /* Square of point 1 Y */ double sx2; /* Square of point 2 X */ double sy2; /* Square of point 2 Y */ double sx3; /* Square of point 3 X */ double sy3; /* Square of point 3 Y */ double xc = 0.0; /* Center point x */ double yc = 0.0; /* Center point y */ double r = 0.0; /* Radius */ double a = 0.0; /* Angle in radians, point 1 */ double b = 0.0; /* Angle in radians, point 3 */ double g = 0.0; /* Angle in radians, point 2 */ double dvp; /* Z component of vector product */ dk4_im_t deltax1; /* X component of first vector */ dk4_im_t deltay1; /* Y component of first vector */ dk4_im_t deltax2; /* X component of second vector */ dk4_im_t deltay2; /* Y component of second vector */ dk4_im_t vp; /* Z component of vector product */ int8_t d = (int8_t)0; /* Positive or negative arc direction */ $? "+ wxdarc_calculation" /* Show arguments in debugging */ $? ". x1 = %d", (int)x1 $? ". y1 = %d", (int)y1 $? ". x2 = %d", (int)x2 $? ". y2 = %d", (int)y2 $? ". x3 = %d", (int)x3 $? ". y3 = %d", (int)y3 /* Initialize result variables */ #if 0 /* 2020-03-24 Modification Turned into initialization. */ xc = 0.0; yc = 0.0; r = 0.0; a = 0.0; b = 0.0; g = 0.0; d = (int8_t)0; #endif /* Convert integers to doubles */ dx1 = (double)x1; dx2 = (double)x2; dx3 = (double)x3; dy1 = (double)y1; dy2 = (double)y2; dy3 = (double)y3; /* Find direction Use integer arithmetics by default, attempt with double arithmetics if integer arithmetics results in overflow. */ dk4error_init(&er); deltax1 = dk4ma_im_sub((dk4_im_t)x2, (dk4_im_t)x1, &er); deltay1 = dk4ma_im_sub((dk4_im_t)y2, (dk4_im_t)y1, &er); deltax2 = dk4ma_im_sub((dk4_im_t)x3, (dk4_im_t)x2, &er); deltay2 = dk4ma_im_sub((dk4_im_t)y3, (dk4_im_t)y2, &er); vp = dk4ma_im_sub( dk4ma_im_mul(deltax1, deltay2, &er), dk4ma_im_mul(deltax2, deltay1, &er), &er ); if (DK4_E_NONE == er.ec) { /* No overflow occured, can use values */ if ((dk4_im_t)0L < vp) { d = (int8_t)1; } else { if ((dk4_im_t)0L > vp) { d = (int8_t)(-1); } } } else { /* On integer overflow attempt calculation in double */ if (DK4_E_MATH_OVERFLOW == er.ec) { dvp = (dx2 - dx1) * (dy3 - dy2) - (dx3 - dx2) * (dy2 - dy1); if (0 != dk4ma_is_finite(dvp)) { if (1.0e-8 < fabs(dvp)) { if (0.0 < dvp) { d = (int8_t)1; } else { d = (int8_t)(-1); } } } } } /* Attempt further calculations only if d is not 0 */ if ((int8_t)0 != d) { /* Calculate squares values */ sx1 = dx1 * dx1; sx2 = dx2 * dx2; sx3 = dx3 * dx3; sy1 = dy1 * dy1; sy2 = dy2 * dy2; sy3 = dy3 * dy3; /* Calculate determinants */ n = 2.0 * ((dx2 - dx1) * dy3 + (dx1 - dx3) * dy2 + (dx3 - dx2) * dy1); zx = (dy2-dy1)*sy3 + (sx1-sx2+sy1-sy2)*dy3 + dy1*sy2 + (sx3-sx1-sy1)*dy2 + (sx2-sx3)*dy1; zx *= -1.0; zy = (dx2-dx1)*sy3 + (dx1-dx3)*sy2 + (dx3-dx2)*sy1 + (dx2-dx1)*sx3 + (sx1-sx2)*dx3 + dx1*sx2 - sx1*dx2; /* Calculate center point */ xc = zx / n; yc = zy / n; /* Radius and angles */ r = sqrt((dx1-xc)*(dx1-xc)+(dy1-yc)*(dy1-yc)); a = atan2((dy1-yc), (dx1-xc)); b = atan2((dy3-yc), (dx3-xc)); g = atan2((dy2-yc), (dx2-xc)); /* Check for infinite and undefined values (no solution or multiple solutions will result in division by 0) */ if (0 == dk4ma_is_finite(xc)) { d = (int8_t)0; } if (0 == dk4ma_is_finite(yc)) { d = (int8_t)0; } if (0 == dk4ma_is_finite(r)) { d = (int8_t)0; } if (0 == dk4ma_is_finite(a)) { d = (int8_t)0; } if (0 == dk4ma_is_finite(b)) { d = (int8_t)0; } if (0 == dk4ma_is_finite(g)) { d = (int8_t)0; } /* Angle corrections depending on arc direction */ if ((int8_t)0 < d) { while (g < a) { g += (2.0 * M_PI); } while (b < g) { b += (2.0 * M_PI); } } else { if ((int8_t)0 > d) { while (g > a) { g -= (2.0 * M_PI); } while (b > g) { b -= (2.0 * M_PI); } } else { xc = 0.0; yc = 0.0; r = 0.0; a = 0.0; b = 0.0; } } } /* Transfer results to destination variables */ if (NULL != pxc) { *pxc = xc; } if (NULL != pyc) { *pyc = yc; } if (NULL != pr) { *pr = r; } if (NULL != pa) { *pa = a; } if (NULL != pb) { *pb = b; } if (NULL != pd) { *pd = d; } /* Show results in debugging */ $? ". xc = %g", xc $? ". yc = %g", yc $? ". r = %g", r $? ". a = %g", a $? ". b = %g", b $? ". g = %g", g $? "- wxdarc_calculation" } #if TEST_WXDARC /* To build test program: gcc -DTEST_WXDARC=1 -o wxdarc -I . wxdarc.c -ldk4c -ldk4ma -ldk4base -lm */ int main(int argc, char *argv[]) { int iv[6]; double xc = 1024.0; double yc = 1024.0; double r = -0.5; double a = 27.0; double b = 128.0; int8_t d = -100; int32_t x1 = 1024L; int32_t y1 = 0L; int32_t x2 = 724L; int32_t y2 = 724L; int32_t x3 = 0L; int32_t y3 = 1024L; int ok = 1; int i; if (7 == argc) { for (i = 0; i < 6; i++) { if (0 == sscanf(argv[i + 1], "%d", &(iv[i]))) { ok = 0; } } if (0 != ok) { x1 = (int32_t)(iv[0]); y1 = (int32_t)(iv[1]); x2 = (int32_t)(iv[2]); y2 = (int32_t)(iv[3]); x3 = (int32_t)(iv[4]); y3 = (int32_t)(iv[5]); } } wxdarc_calculation( &xc, &yc, &r, &a, &b, &d, x1, y1, x2, y2, x3, y3 ); printf("x1 = %d\n", (int)x1); printf("y1 = %d\n", (int)y1); printf("x2 = %d\n", (int)x2); printf("y2 = %d\n", (int)y2); printf("x3 = %d\n", (int)x3); printf("y3 = %d\n", (int)y3); printf("x = %g\n", xc); printf("y = %g\n", yc); printf("r = %g\n", r); printf("a = %g\n", a); printf("b = %g\n", b); printf("d = %d\n", (int)d); return 0; } #endif /* vim: set ai sw=4 ts=4 : */