1// This file is part of SmallBASIC
2//
3// SmallBASIC LIBRARY - extra geometry algorithms
4//
5// This program is distributed under the terms of the GPL v2.0 or later
6// Download the GNU Public License (GPL) from www.gnu.org
7//
8// Copyright(C) 2000 Nicholas Christopoulos
9
10#include "common/geom.h"
11#include "common/kw.h"
12
13/**
14 * length of line segment A-B
15 */
16double geo_seglen(double Ax, double Ay, double Bx, double By) {
17 double dx = Bx - Ax;
18 double dy = By - Ay;
19
20 return sqrt(dx * dx + dy * dy);
21}
22
23/**
24 * distance of point C from line (A,B)
25 */
26double geo_distfromline(double Ax, double Ay, double Bx, double By, double Cx, double Cy) {
27 double L, s;
28// double r;
29// double Px, Py; // projection of C on AB
30
31 L = geo_seglen(Ax, Ay, Bx, By);
32
33 if (L == 0.0) // A,B are the same point, TODO: tolerance
34 return geo_seglen(Ax, Ay, Cx, Cy);
35
36// r = ((Ay-Cy) * (Ay-By) - (Ax-Cx) * (Bx-Ax)) / (L * L);
37// Px = Ax + r * (Bx - Ax);
38// Py = Ay + r * (By - Ay);
39
40 s = ((Ay - Cy) * (Bx - Ax) - (Ax - Cx) * (By - Ay)) / (L * L);
41 return s * L;
42}
43
44/**
45 * distance of point C from line segment (A->B)
46 */
47double geo_distfromseg(double Ax, double Ay, double Bx, double By, double Cx, double Cy) {
48// double Px, Py; // projection of C on AB
49 double L, r, s;
50 double ca, cb;
51
52 L = geo_seglen(Ax, Ay, Bx, By);
53
54 if (L == 0.0) // A,B are the same point, TODO: tolerance
55 return geo_seglen(Ax, Ay, Cx, Cy);
56
57 r = ((Ay - Cy) * (Ay - By) - (Ax - Cx) * (Bx - Ax)) / (L * L);
58
59 if (r >= 0.0 && r <= 1.0) { // TODO: tolerance
60// Px = Ax + r * (Bx - Ax);
61// Py = Ay + r * (By - Ay);
62
63 s = ((Ay - Cy) * (Bx - Ax) - (Ax - Cx) * (By - Ay)) / (L * L);
64 return s * L;
65 }
66
67 // else the P is out of A-B
68 ca = geo_seglen(Ax, Ay, Cx, Cy);
69 cb = geo_seglen(Bx, By, Cx, Cy);
70 return (ca < cb) ? ca : cb;
71}
72
73/*
74 * COS/SIN of two line segments
75 */
76double geo_segangle(int type, double Adx, double Ady, double Bdx, double Bdy) {
77 double la, lb;
78 double UAdx, UAdy, UBdx, UBdy;
79
80 la = sqrt(Adx * Adx + Ady * Ady);
81 lb = sqrt(Bdx * Bdx + Bdy * Bdy);
82
83 if (la != 0 && lb != 0) {
84 UAdx = Adx / la;
85 UAdy = Ady / la;
86 UBdx = Bdx / lb;
87 UBdy = Bdy / lb;
88 if (type == kwSEGCOS)
89 return UAdx * UBdx + UAdy * UBdy;
90 else
91 return UAdy * UBdx + UAdx * UBdy;
92 }
93
94 return 0; // ��� ��� �� ���������� ���� ����� 0.
95}
96
97/*********************************************************************
98 ANSI C code from the article "Centroid of a Polygon"
99 by Gerard Bashein and Paul R. Detmer
100
101 polyCentroid: Calculates the centroid (xCentroid, yCentroid) and area
102 of a polygon, given its vertices (x[0], y[0]) ... (x[n-1], y[n-1]). It
103 is assumed that the contour is closed, i.e., that the vertex following
104 (x[n-1], y[n-1]) is (x[0], y[0]). The algebraic sign of the area is
105 positive for counterclockwise ordering of vertices in x-y plane;
106 otherwise negative.
107
108 Returned values: 0 for normal execution; 1 if the polygon is
109 degenerate (number of vertices < 3); and 2 if area = 0 (and the
110 centroid is undefined).
111 **********************************************************************/
112int geo_polycentroid(pt_t * poly, int n, var_num_t *xCentroid, var_num_t *yCentroid, var_num_t *area) {
113 int i, j;
114 var_num_t ai, atmp = 0, xtmp = 0, ytmp = 0;
115
116 if (n < 3) {
117 return 1;
118 }
119
120 for (i = n - 1, j = 0; j < n; i = j, j++) {
121 ai = poly[i].x * poly[j].y - poly[j].x * poly[i].y;
122 atmp += ai;
123 xtmp += (poly[j].x + poly[i].x) * ai;
124 ytmp += (poly[j].y + poly[i].y) * ai;
125 }
126
127 *area = atmp / 2;
128 if (atmp != 0) {
129 *xCentroid = xtmp / (3 * atmp);
130 *yCentroid = ytmp / (3 * atmp);
131 return 0;
132 }
133
134 return 2;
135}
136