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