agg_rust/
math.rs

1//! Geometric math utilities.
2//!
3//! Port of `agg_math.h` and `agg_sqrt_tables.cpp` — distances, intersections,
4//! cross products, triangle operations, fast integer sqrt, and Bessel functions.
5
6// ============================================================================
7// Constants
8// ============================================================================
9
10/// Coinciding points maximal distance (epsilon).
11pub const VERTEX_DIST_EPSILON: f64 = 1e-14;
12
13/// Epsilon for intersection calculations.
14pub const INTERSECTION_EPSILON: f64 = 1.0e-30;
15
16// ============================================================================
17// Cross product and point-in-triangle
18// ============================================================================
19
20/// Cross product of vectors (x2-x1, y2-y1) and (x-x2, y-y2).
21/// The sign indicates which side of the line (x1,y1)→(x2,y2) the point (x,y) is on.
22#[inline]
23pub fn cross_product(x1: f64, y1: f64, x2: f64, y2: f64, x: f64, y: f64) -> f64 {
24    (x - x2) * (y2 - y1) - (y - y2) * (x2 - x1)
25}
26
27/// Test if point (x, y) is inside triangle (x1,y1), (x2,y2), (x3,y3).
28#[inline]
29#[allow(clippy::too_many_arguments)]
30pub fn point_in_triangle(
31    x1: f64,
32    y1: f64,
33    x2: f64,
34    y2: f64,
35    x3: f64,
36    y3: f64,
37    x: f64,
38    y: f64,
39) -> bool {
40    let cp1 = cross_product(x1, y1, x2, y2, x, y) < 0.0;
41    let cp2 = cross_product(x2, y2, x3, y3, x, y) < 0.0;
42    let cp3 = cross_product(x3, y3, x1, y1, x, y) < 0.0;
43    cp1 == cp2 && cp2 == cp3 && cp3 == cp1
44}
45
46// ============================================================================
47// Distance calculations
48// ============================================================================
49
50/// Euclidean distance between two points.
51#[inline]
52pub fn calc_distance(x1: f64, y1: f64, x2: f64, y2: f64) -> f64 {
53    let dx = x2 - x1;
54    let dy = y2 - y1;
55    (dx * dx + dy * dy).sqrt()
56}
57
58/// Squared Euclidean distance between two points.
59#[inline]
60pub fn calc_sq_distance(x1: f64, y1: f64, x2: f64, y2: f64) -> f64 {
61    let dx = x2 - x1;
62    let dy = y2 - y1;
63    dx * dx + dy * dy
64}
65
66/// Signed distance from point (x, y) to the infinite line through (x1,y1)→(x2,y2).
67/// Positive means left side, negative means right side.
68/// If the line segment is degenerate (length < VERTEX_DIST_EPSILON), returns
69/// the distance from (x,y) to (x1,y1).
70#[inline]
71pub fn calc_line_point_distance(x1: f64, y1: f64, x2: f64, y2: f64, x: f64, y: f64) -> f64 {
72    let dx = x2 - x1;
73    let dy = y2 - y1;
74    let d = (dx * dx + dy * dy).sqrt();
75    if d < VERTEX_DIST_EPSILON {
76        return calc_distance(x1, y1, x, y);
77    }
78    ((x - x2) * dy - (y - y2) * dx) / d
79}
80
81/// Compute the parameter `u` for the projection of point (x, y) onto
82/// the line segment (x1,y1)→(x2,y2). Returns 0 if the segment is degenerate.
83#[inline]
84pub fn calc_segment_point_u(x1: f64, y1: f64, x2: f64, y2: f64, x: f64, y: f64) -> f64 {
85    let dx = x2 - x1;
86    let dy = y2 - y1;
87
88    if dx == 0.0 && dy == 0.0 {
89        return 0.0;
90    }
91
92    let pdx = x - x1;
93    let pdy = y - y1;
94
95    (pdx * dx + pdy * dy) / (dx * dx + dy * dy)
96}
97
98/// Squared distance from point (x, y) to the closest point on segment
99/// (x1,y1)→(x2,y2), given pre-computed parameter `u`.
100#[inline]
101pub fn calc_segment_point_sq_distance_with_u(
102    x1: f64,
103    y1: f64,
104    x2: f64,
105    y2: f64,
106    x: f64,
107    y: f64,
108    u: f64,
109) -> f64 {
110    if u <= 0.0 {
111        calc_sq_distance(x, y, x1, y1)
112    } else if u >= 1.0 {
113        calc_sq_distance(x, y, x2, y2)
114    } else {
115        calc_sq_distance(x, y, x1 + u * (x2 - x1), y1 + u * (y2 - y1))
116    }
117}
118
119/// Squared distance from point (x, y) to the closest point on segment
120/// (x1,y1)→(x2,y2).
121#[inline]
122pub fn calc_segment_point_sq_distance(x1: f64, y1: f64, x2: f64, y2: f64, x: f64, y: f64) -> f64 {
123    calc_segment_point_sq_distance_with_u(
124        x1,
125        y1,
126        x2,
127        y2,
128        x,
129        y,
130        calc_segment_point_u(x1, y1, x2, y2, x, y),
131    )
132}
133
134// ============================================================================
135// Intersection
136// ============================================================================
137
138/// Calculate the intersection point of two line segments:
139/// (ax,ay)→(bx,by) and (cx,cy)→(dx,dy).
140/// Returns `Some((x, y))` if they intersect, `None` if parallel.
141#[inline]
142#[allow(clippy::too_many_arguments)]
143pub fn calc_intersection(
144    ax: f64,
145    ay: f64,
146    bx: f64,
147    by: f64,
148    cx: f64,
149    cy: f64,
150    dx: f64,
151    dy: f64,
152) -> Option<(f64, f64)> {
153    let num = (ay - cy) * (dx - cx) - (ax - cx) * (dy - cy);
154    let den = (bx - ax) * (dy - cy) - (by - ay) * (dx - cx);
155    if den.abs() < INTERSECTION_EPSILON {
156        return None;
157    }
158    let r = num / den;
159    Some((ax + r * (bx - ax), ay + r * (by - ay)))
160}
161
162/// Quick check whether two line segments (x1,y1)→(x2,y2) and
163/// (x3,y3)→(x4,y4) intersect (boundary excluded).
164#[inline]
165#[allow(clippy::too_many_arguments)]
166pub fn intersection_exists(
167    x1: f64,
168    y1: f64,
169    x2: f64,
170    y2: f64,
171    x3: f64,
172    y3: f64,
173    x4: f64,
174    y4: f64,
175) -> bool {
176    let dx1 = x2 - x1;
177    let dy1 = y2 - y1;
178    let dx2 = x4 - x3;
179    let dy2 = y4 - y3;
180    ((x3 - x2) * dy1 - (y3 - y2) * dx1 < 0.0) != ((x4 - x2) * dy1 - (y4 - y2) * dx1 < 0.0)
181        && ((x1 - x4) * dy2 - (y1 - y4) * dx2 < 0.0) != ((x2 - x4) * dy2 - (y2 - y4) * dx2 < 0.0)
182}
183
184// ============================================================================
185// Orthogonal and triangle operations
186// ============================================================================
187
188/// Calculate the orthogonal displacement vector of magnitude `thickness`
189/// perpendicular to the line (x1,y1)→(x2,y2).
190#[inline]
191pub fn calc_orthogonal(thickness: f64, x1: f64, y1: f64, x2: f64, y2: f64) -> (f64, f64) {
192    let dx = x2 - x1;
193    let dy = y2 - y1;
194    let d = (dx * dx + dy * dy).sqrt();
195    (thickness * dy / d, -thickness * dx / d)
196}
197
198/// Dilate a triangle by distance `d`, producing 6 output points
199/// (two per edge). Returns `([x0..x5], [y0..y5])`.
200pub fn dilate_triangle(
201    x1: f64,
202    y1: f64,
203    x2: f64,
204    y2: f64,
205    x3: f64,
206    y3: f64,
207    d: f64,
208) -> ([f64; 6], [f64; 6]) {
209    let mut dx1 = 0.0;
210    let mut dy1 = 0.0;
211    let mut dx2 = 0.0;
212    let mut dy2 = 0.0;
213    let mut dx3 = 0.0;
214    let mut dy3 = 0.0;
215    let mut d = d;
216    let loc = cross_product(x1, y1, x2, y2, x3, y3);
217    if loc.abs() > INTERSECTION_EPSILON {
218        if cross_product(x1, y1, x2, y2, x3, y3) > 0.0 {
219            d = -d;
220        }
221        let o1 = calc_orthogonal(d, x1, y1, x2, y2);
222        dx1 = o1.0;
223        dy1 = o1.1;
224        let o2 = calc_orthogonal(d, x2, y2, x3, y3);
225        dx2 = o2.0;
226        dy2 = o2.1;
227        let o3 = calc_orthogonal(d, x3, y3, x1, y1);
228        dx3 = o3.0;
229        dy3 = o3.1;
230    }
231    let x = [x1 + dx1, x2 + dx1, x2 + dx2, x3 + dx2, x3 + dx3, x1 + dx3];
232    let y = [y1 + dy1, y2 + dy1, y2 + dy2, y3 + dy2, y3 + dy3, y1 + dy3];
233    (x, y)
234}
235
236/// Signed area of triangle (x1,y1), (x2,y2), (x3,y3).
237#[inline]
238pub fn calc_triangle_area(x1: f64, y1: f64, x2: f64, y2: f64, x3: f64, y3: f64) -> f64 {
239    (x1 * y2 - x2 * y1 + x2 * y3 - x3 * y2 + x3 * y1 - x1 * y3) * 0.5
240}
241
242/// Signed area of a polygon defined by a slice of points with `x` and `y` fields.
243pub fn calc_polygon_area(vertices: &[crate::basics::PointD]) -> f64 {
244    if vertices.is_empty() {
245        return 0.0;
246    }
247    let mut sum = 0.0;
248    let mut x = vertices[0].x;
249    let mut y = vertices[0].y;
250    let xs = x;
251    let ys = y;
252
253    for v in &vertices[1..] {
254        sum += x * v.y - y * v.x;
255        x = v.x;
256        y = v.y;
257    }
258    (sum + x * ys - y * xs) * 0.5
259}
260
261/// Signed area of a polygon defined by a slice of `VertexDist`.
262/// Same algorithm as `calc_polygon_area` but works with `VertexDist` slices
263/// (needed by `vcgen_contour`).
264pub fn calc_polygon_area_vd(vertices: &[crate::array::VertexDist]) -> f64 {
265    if vertices.is_empty() {
266        return 0.0;
267    }
268    let mut sum = 0.0;
269    let mut x = vertices[0].x;
270    let mut y = vertices[0].y;
271    let xs = x;
272    let ys = y;
273
274    for v in &vertices[1..] {
275        sum += x * v.y - y * v.x;
276        x = v.x;
277        y = v.y;
278    }
279    (sum + x * ys - y * xs) * 0.5
280}
281
282// ============================================================================
283// Fast integer square root (lookup table based)
284// ============================================================================
285
286/// Lookup table for fast integer sqrt. 1024 entries.
287/// Port of `g_sqrt_table` from `agg_sqrt_tables.cpp`.
288#[rustfmt::skip]
289static SQRT_TABLE: [u16; 1024] = [
290    0,
291    2048,2896,3547,4096,4579,5017,5418,5793,6144,6476,6792,7094,7384,7663,7932,8192,8444,
292    8689,8927,9159,9385,9606,9822,10033,10240,10443,10642,10837,11029,11217,11403,11585,
293    11765,11942,12116,12288,12457,12625,12790,12953,13114,13273,13430,13585,13738,13890,
294    14040,14189,14336,14482,14626,14768,14910,15050,15188,15326,15462,15597,15731,15864,
295    15995,16126,16255,16384,16512,16638,16764,16888,17012,17135,17257,17378,17498,17618,
296    17736,17854,17971,18087,18203,18318,18432,18545,18658,18770,18882,18992,19102,19212,
297    19321,19429,19537,19644,19750,19856,19961,20066,20170,20274,20377,20480,20582,20684,
298    20785,20886,20986,21085,21185,21283,21382,21480,21577,21674,21771,21867,21962,22058,
299    22153,22247,22341,22435,22528,22621,22713,22806,22897,22989,23080,23170,23261,23351,
300    23440,23530,23619,23707,23796,23884,23971,24059,24146,24232,24319,24405,24491,24576,
301    24661,24746,24831,24915,24999,25083,25166,25249,25332,25415,25497,25580,25661,25743,
302    25824,25905,25986,26067,26147,26227,26307,26387,26466,26545,26624,26703,26781,26859,
303    26937,27015,27092,27170,27247,27324,27400,27477,27553,27629,27705,27780,27856,27931,
304    28006,28081,28155,28230,28304,28378,28452,28525,28599,28672,28745,28818,28891,28963,
305    29035,29108,29180,29251,29323,29394,29466,29537,29608,29678,29749,29819,29890,29960,
306    30030,30099,30169,30238,30308,30377,30446,30515,30583,30652,30720,30788,30856,30924,
307    30992,31059,31127,31194,31261,31328,31395,31462,31529,31595,31661,31727,31794,31859,
308    31925,31991,32056,32122,32187,32252,32317,32382,32446,32511,32575,32640,32704,32768,
309    32832,32896,32959,33023,33086,33150,33213,33276,33339,33402,33465,33527,33590,33652,
310    33714,33776,33839,33900,33962,34024,34086,34147,34208,34270,34331,34392,34453,34514,
311    34574,34635,34695,34756,34816,34876,34936,34996,35056,35116,35176,35235,35295,35354,
312    35413,35472,35531,35590,35649,35708,35767,35825,35884,35942,36001,36059,36117,36175,
313    36233,36291,36348,36406,36464,36521,36578,36636,36693,36750,36807,36864,36921,36978,
314    37034,37091,37147,37204,37260,37316,37372,37429,37485,37540,37596,37652,37708,37763,
315    37819,37874,37929,37985,38040,38095,38150,38205,38260,38315,38369,38424,38478,38533,
316    38587,38642,38696,38750,38804,38858,38912,38966,39020,39073,39127,39181,39234,39287,
317    39341,39394,39447,39500,39553,39606,39659,39712,39765,39818,39870,39923,39975,40028,
318    40080,40132,40185,40237,40289,40341,40393,40445,40497,40548,40600,40652,40703,40755,
319    40806,40857,40909,40960,41011,41062,41113,41164,41215,41266,41317,41368,41418,41469,
320    41519,41570,41620,41671,41721,41771,41821,41871,41922,41972,42021,42071,42121,42171,
321    42221,42270,42320,42369,42419,42468,42518,42567,42616,42665,42714,42763,42813,42861,
322    42910,42959,43008,43057,43105,43154,43203,43251,43300,43348,43396,43445,43493,43541,
323    43589,43637,43685,43733,43781,43829,43877,43925,43972,44020,44068,44115,44163,44210,
324    44258,44305,44352,44400,44447,44494,44541,44588,44635,44682,44729,44776,44823,44869,
325    44916,44963,45009,45056,45103,45149,45195,45242,45288,45334,45381,45427,45473,45519,
326    45565,45611,45657,45703,45749,45795,45840,45886,45932,45977,46023,46069,46114,46160,
327    46205,46250,46296,46341,46386,46431,46477,46522,46567,46612,46657,46702,46746,46791,
328    46836,46881,46926,46970,47015,47059,47104,47149,47193,47237,47282,47326,47370,47415,
329    47459,47503,47547,47591,47635,47679,47723,47767,47811,47855,47899,47942,47986,48030,
330    48074,48117,48161,48204,48248,48291,48335,48378,48421,48465,48508,48551,48594,48637,
331    48680,48723,48766,48809,48852,48895,48938,48981,49024,49067,49109,49152,49195,49237,
332    49280,49322,49365,49407,49450,49492,49535,49577,49619,49661,49704,49746,49788,49830,
333    49872,49914,49956,49998,50040,50082,50124,50166,50207,50249,50291,50332,50374,50416,
334    50457,50499,50540,50582,50623,50665,50706,50747,50789,50830,50871,50912,50954,50995,
335    51036,51077,51118,51159,51200,51241,51282,51323,51364,51404,51445,51486,51527,51567,
336    51608,51649,51689,51730,51770,51811,51851,51892,51932,51972,52013,52053,52093,52134,
337    52174,52214,52254,52294,52334,52374,52414,52454,52494,52534,52574,52614,52654,52694,
338    52734,52773,52813,52853,52892,52932,52972,53011,53051,53090,53130,53169,53209,53248,
339    53287,53327,53366,53405,53445,53484,53523,53562,53601,53640,53679,53719,53758,53797,
340    53836,53874,53913,53952,53991,54030,54069,54108,54146,54185,54224,54262,54301,54340,
341    54378,54417,54455,54494,54532,54571,54609,54647,54686,54724,54762,54801,54839,54877,
342    54915,54954,54992,55030,55068,55106,55144,55182,55220,55258,55296,55334,55372,55410,
343    55447,55485,55523,55561,55599,55636,55674,55712,55749,55787,55824,55862,55900,55937,
344    55975,56012,56049,56087,56124,56162,56199,56236,56273,56311,56348,56385,56422,56459,
345    56497,56534,56571,56608,56645,56682,56719,56756,56793,56830,56867,56903,56940,56977,
346    57014,57051,57087,57124,57161,57198,57234,57271,57307,57344,57381,57417,57454,57490,
347    57527,57563,57599,57636,57672,57709,57745,57781,57817,57854,57890,57926,57962,57999,
348    58035,58071,58107,58143,58179,58215,58251,58287,58323,58359,58395,58431,58467,58503,
349    58538,58574,58610,58646,58682,58717,58753,58789,58824,58860,58896,58931,58967,59002,
350    59038,59073,59109,59144,59180,59215,59251,59286,59321,59357,59392,59427,59463,59498,
351    59533,59568,59603,59639,59674,59709,59744,59779,59814,59849,59884,59919,59954,59989,
352    60024,60059,60094,60129,60164,60199,60233,60268,60303,60338,60373,60407,60442,60477,
353    60511,60546,60581,60615,60650,60684,60719,60753,60788,60822,60857,60891,60926,60960,
354    60995,61029,61063,61098,61132,61166,61201,61235,61269,61303,61338,61372,61406,61440,
355    61474,61508,61542,61576,61610,61644,61678,61712,61746,61780,61814,61848,61882,61916,
356    61950,61984,62018,62051,62085,62119,62153,62186,62220,62254,62287,62321,62355,62388,
357    62422,62456,62489,62523,62556,62590,62623,62657,62690,62724,62757,62790,62824,62857,
358    62891,62924,62957,62991,63024,63057,63090,63124,63157,63190,63223,63256,63289,63323,
359    63356,63389,63422,63455,63488,63521,63554,63587,63620,63653,63686,63719,63752,63785,
360    63817,63850,63883,63916,63949,63982,64014,64047,64080,64113,64145,64178,64211,64243,
361    64276,64309,64341,64374,64406,64439,64471,64504,64536,64569,64601,64634,64666,64699,
362    64731,64763,64796,64828,64861,64893,64925,64957,64990,65022,65054,65086,65119,65151,
363    65183,65215,65247,65279,65312,65344,65376,65408,65440,65472,65504,
364];
365
366/// Lookup table mapping byte values to their most significant bit position.
367/// Port of `g_elder_bit_table` from `agg_sqrt_tables.cpp`.
368#[rustfmt::skip]
369static ELDER_BIT_TABLE: [i8; 256] = [
370    0,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,
371    5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
372    6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
373    6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
374    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
375    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
376    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
377    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
378];
379
380/// Fast integer square root using lookup tables.
381/// No divisions, multiplications, or loops — just bit shifts and table lookups.
382/// Port of C++ `fast_sqrt` (portable C path).
383pub fn fast_sqrt(val: u32) -> u32 {
384    let t = val;
385    let mut shift: i32 = 11;
386
387    let bit: i32;
388    let b = (t >> 24) as u8;
389    if b != 0 {
390        bit = ELDER_BIT_TABLE[b as usize] as i32 + 24;
391    } else {
392        let b = ((t >> 16) & 0xFF) as u8;
393        if b != 0 {
394            bit = ELDER_BIT_TABLE[b as usize] as i32 + 16;
395        } else {
396            let b = ((t >> 8) & 0xFF) as u8;
397            if b != 0 {
398                bit = ELDER_BIT_TABLE[b as usize] as i32 + 8;
399            } else {
400                bit = ELDER_BIT_TABLE[t as u8 as usize] as i32;
401            }
402        }
403    }
404
405    let mut val = val;
406    let bit = bit - 9;
407    if bit > 0 {
408        let half_bit = (bit >> 1) + (bit & 1);
409        shift -= half_bit;
410        val >>= (half_bit << 1) as u32;
411    }
412    (SQRT_TABLE[val as usize] as u32) >> shift as u32
413}
414
415// ============================================================================
416// Bessel function
417// ============================================================================
418
419/// Bessel function of the first kind of order `n`.
420///
421/// Adapted for AGG by Andy Wilk (castor.vulgaris@gmail.com).
422/// Originally from C++ Mathematical Library (Gareth Walker).
423pub fn besj(x: f64, n: i32) -> f64 {
424    if n < 0 {
425        return 0.0;
426    }
427    let d = 1e-6;
428    let mut b = 0.0;
429    if x.abs() <= d {
430        if n != 0 {
431            return 0.0;
432        }
433        return 1.0;
434    }
435    let mut b1 = 0.0;
436    let mut m1 = (x.abs() + 6.0) as i32;
437    if x.abs() > 5.0 {
438        m1 = (1.4 * x.abs() + 60.0 / x.abs()) as i32;
439    }
440    let mut m2 = (n as f64 + 2.0 + x.abs() / 4.0) as i32;
441    if m1 > m2 {
442        m2 = m1;
443    }
444
445    loop {
446        let mut c3 = 0.0;
447        let mut c2 = 1e-30;
448        let mut c4 = 0.0;
449        let mut m8 = 1;
450        if m2 / 2 * 2 == m2 {
451            m8 = -1;
452        }
453        let imax = m2 - 2;
454        for i in 1..=imax {
455            let c6 = 2.0 * (m2 - i) as f64 * c2 / x - c3;
456            c3 = c2;
457            c2 = c6;
458            if m2 - i - 1 == n {
459                b = c6;
460            }
461            m8 = -m8;
462            if m8 > 0 {
463                c4 += 2.0 * c6;
464            }
465        }
466        let c6 = 2.0 * c2 / x - c3;
467        if n == 0 {
468            b = c6;
469        }
470        c4 += c6;
471        b /= c4;
472        if (b - b1).abs() < d {
473            return b;
474        }
475        b1 = b;
476        m2 += 3;
477    }
478}
479
480// ============================================================================
481// Tests
482// ============================================================================
483
484#[cfg(test)]
485mod tests {
486    use super::*;
487
488    const EPSILON: f64 = 1e-10;
489
490    #[test]
491    fn test_cross_product() {
492        // Point on the line: cross product should be 0
493        let cp = cross_product(0.0, 0.0, 1.0, 0.0, 2.0, 0.0);
494        assert!(cp.abs() < EPSILON);
495
496        // Point above the line (left side): negative cross product
497        // Formula: (x-x2)*(y2-y1) - (y-y2)*(x2-x1) = (-0.5)*0 - (1.0)*(1.0) = -1.0
498        let cp = cross_product(0.0, 0.0, 1.0, 0.0, 0.5, 1.0);
499        assert!(cp < 0.0);
500
501        // Point below the line (right side): positive cross product
502        let cp = cross_product(0.0, 0.0, 1.0, 0.0, 0.5, -1.0);
503        assert!(cp > 0.0);
504    }
505
506    #[test]
507    fn test_point_in_triangle() {
508        // Center of unit triangle
509        assert!(point_in_triangle(0.0, 0.0, 1.0, 0.0, 0.5, 1.0, 0.5, 0.3));
510        // Outside
511        assert!(!point_in_triangle(0.0, 0.0, 1.0, 0.0, 0.5, 1.0, 2.0, 2.0));
512    }
513
514    #[test]
515    fn test_calc_distance() {
516        assert!((calc_distance(0.0, 0.0, 3.0, 4.0) - 5.0).abs() < EPSILON);
517        assert!((calc_distance(0.0, 0.0, 0.0, 0.0)).abs() < EPSILON);
518        assert!((calc_distance(1.0, 1.0, 1.0, 1.0)).abs() < EPSILON);
519    }
520
521    #[test]
522    fn test_calc_sq_distance() {
523        assert!((calc_sq_distance(0.0, 0.0, 3.0, 4.0) - 25.0).abs() < EPSILON);
524    }
525
526    #[test]
527    fn test_calc_line_point_distance() {
528        // Point (0, 1) relative to line (0,0)→(1,0):
529        // ((x-x2)*dy - (y-y2)*dx) / d = ((0-1)*0 - (1-0)*1) / 1 = -1.0
530        let d = calc_line_point_distance(0.0, 0.0, 1.0, 0.0, 0.0, 1.0);
531        assert!((d - (-1.0)).abs() < EPSILON);
532
533        // Point (0, -1) below the line: positive
534        let d = calc_line_point_distance(0.0, 0.0, 1.0, 0.0, 0.0, -1.0);
535        assert!((d - 1.0).abs() < EPSILON);
536    }
537
538    #[test]
539    fn test_calc_segment_point_u() {
540        // Midpoint of segment
541        let u = calc_segment_point_u(0.0, 0.0, 2.0, 0.0, 1.0, 0.0);
542        assert!((u - 0.5).abs() < EPSILON);
543
544        // Before start
545        let u = calc_segment_point_u(0.0, 0.0, 2.0, 0.0, -1.0, 0.0);
546        assert!(u < 0.0);
547
548        // After end
549        let u = calc_segment_point_u(0.0, 0.0, 2.0, 0.0, 3.0, 0.0);
550        assert!(u > 1.0);
551
552        // Degenerate segment
553        let u = calc_segment_point_u(0.0, 0.0, 0.0, 0.0, 1.0, 1.0);
554        assert_eq!(u, 0.0);
555    }
556
557    #[test]
558    fn test_calc_segment_point_sq_distance() {
559        // Distance to midpoint of horizontal segment from point above
560        let d = calc_segment_point_sq_distance(0.0, 0.0, 2.0, 0.0, 1.0, 1.0);
561        assert!((d - 1.0).abs() < EPSILON);
562
563        // Distance to start (before segment)
564        let d = calc_segment_point_sq_distance(0.0, 0.0, 2.0, 0.0, -1.0, 0.0);
565        assert!((d - 1.0).abs() < EPSILON);
566    }
567
568    #[test]
569    fn test_calc_intersection() {
570        // Two perpendicular lines crossing at (1, 1)
571        let result = calc_intersection(0.0, 1.0, 2.0, 1.0, 1.0, 0.0, 1.0, 2.0);
572        assert!(result.is_some());
573        let (x, y) = result.unwrap();
574        assert!((x - 1.0).abs() < EPSILON);
575        assert!((y - 1.0).abs() < EPSILON);
576
577        // Parallel lines
578        let result = calc_intersection(0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0);
579        assert!(result.is_none());
580    }
581
582    #[test]
583    fn test_intersection_exists() {
584        // Crossing segments
585        assert!(intersection_exists(0.0, 0.0, 2.0, 2.0, 0.0, 2.0, 2.0, 0.0));
586
587        // Non-crossing segments
588        assert!(!intersection_exists(0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0));
589    }
590
591    #[test]
592    fn test_calc_triangle_area() {
593        // Unit right triangle: area = 0.5
594        let area = calc_triangle_area(0.0, 0.0, 1.0, 0.0, 0.0, 1.0);
595        assert!((area - 0.5).abs() < EPSILON);
596
597        // Reversed winding: negative area
598        let area = calc_triangle_area(0.0, 0.0, 0.0, 1.0, 1.0, 0.0);
599        assert!((area - (-0.5)).abs() < EPSILON);
600    }
601
602    #[test]
603    fn test_calc_polygon_area() {
604        use crate::basics::PointD;
605        // Unit square: area = 1.0
606        let square = vec![
607            PointD::new(0.0, 0.0),
608            PointD::new(1.0, 0.0),
609            PointD::new(1.0, 1.0),
610            PointD::new(0.0, 1.0),
611        ];
612        let area = calc_polygon_area(&square);
613        assert!((area - 1.0).abs() < EPSILON);
614    }
615
616    #[test]
617    fn test_calc_polygon_area_vd() {
618        use crate::array::VertexDist;
619        // Unit square: area = 1.0
620        let square = vec![
621            VertexDist::new(0.0, 0.0),
622            VertexDist::new(1.0, 0.0),
623            VertexDist::new(1.0, 1.0),
624            VertexDist::new(0.0, 1.0),
625        ];
626        let area = calc_polygon_area_vd(&square);
627        assert!((area - 1.0).abs() < EPSILON);
628
629        // Empty: area = 0
630        let empty: Vec<VertexDist> = vec![];
631        assert_eq!(calc_polygon_area_vd(&empty), 0.0);
632    }
633
634    #[test]
635    fn test_calc_polygon_area_vd_ccw() {
636        use crate::array::VertexDist;
637        // CCW square: area = -1.0
638        let square = vec![
639            VertexDist::new(0.0, 0.0),
640            VertexDist::new(0.0, 1.0),
641            VertexDist::new(1.0, 1.0),
642            VertexDist::new(1.0, 0.0),
643        ];
644        let area = calc_polygon_area_vd(&square);
645        assert!((area - (-1.0)).abs() < EPSILON);
646    }
647
648    #[test]
649    fn test_fast_sqrt() {
650        // Test a few known values
651        assert_eq!(fast_sqrt(0), 0);
652        assert_eq!(fast_sqrt(1), 1);
653        assert_eq!(fast_sqrt(4), 2);
654        assert_eq!(fast_sqrt(9), 3);
655        assert_eq!(fast_sqrt(16), 4);
656        assert_eq!(fast_sqrt(100), 10);
657        assert_eq!(fast_sqrt(10000), 100);
658    }
659
660    #[test]
661    fn test_fast_sqrt_accuracy() {
662        // fast_sqrt should be reasonably accurate for values used in AGG
663        for val in [25, 49, 64, 81, 144, 225, 400, 625, 900, 1600, 2500, 10000] {
664            let expected = (val as f64).sqrt().round() as u32;
665            let result = fast_sqrt(val);
666            assert_eq!(
667                result, expected,
668                "fast_sqrt({}) = {}, expected {}",
669                val, result, expected
670            );
671        }
672    }
673
674    #[test]
675    fn test_besj_order_zero() {
676        // J_0(0) = 1
677        assert!((besj(0.0, 0) - 1.0).abs() < 1e-5);
678        // J_0(2.4048...) ≈ 0 (first zero)
679        assert!(besj(2.4048, 0).abs() < 0.001);
680    }
681
682    #[test]
683    fn test_besj_order_one() {
684        // J_1(0) = 0
685        assert!((besj(0.0, 1)).abs() < 1e-5);
686        // J_1(3.8317...) ≈ 0 (first zero)
687        assert!(besj(3.8317, 1).abs() < 0.001);
688    }
689
690    #[test]
691    fn test_besj_negative_order() {
692        assert_eq!(besj(1.0, -1), 0.0);
693    }
694
695    #[test]
696    fn test_calc_orthogonal() {
697        let (dx, dy) = calc_orthogonal(1.0, 0.0, 0.0, 1.0, 0.0);
698        assert!((dx).abs() < EPSILON);
699        assert!((dy - (-1.0)).abs() < EPSILON);
700    }
701
702    #[test]
703    fn test_dilate_triangle() {
704        // Just verify it doesn't panic and produces 6 points
705        let (x, y) = dilate_triangle(0.0, 0.0, 1.0, 0.0, 0.5, 1.0, 0.1);
706        assert_eq!(x.len(), 6);
707        assert_eq!(y.len(), 6);
708    }
709}
agg_rust/math.rs

agg_rust/
math.rs
