git clone https://git.lucas.co/hou-control.git
vex/include/sidefxlabs_mathx.h (27.3K)
1 #ifndef __sidefxlabs_math_h__
2 #define __sidefxlabs_math_h__
3
4
5 #define LABS_2PI 6.28318530718
6 #define LABS_RAD360 6.28318530718
7 #define LABS_PI 3.14159265359
8 #define LABS_RAD180 3.14159265359
9 #define LABS_PI2 1.57079632679
10 #define LABS_RAD90 1.57079632679
11
12 #define LABS_E 2.71828182846
13 #define LABS_SR2 1.41421356237
14
15 // Defines 32-bit tolerances
16 #define LABS_TOL16 0.000001 // Accurate below 16
17 #define LABS_TOL128 0.00001 // Accurate below 128
18 #define LABS_TOL1K 0.0001 // Accurate below 1024
19 #define LABS_TOL 0.0001 // Accurate below 1024 (default tolerance)
20 #define LABS_TOL16K 0.001 // Accurate below 16384
21 #define LABS_TOL130K 0.01 // Accurate below 131072
22 #define LABS_TOL1M 0.1 // Accurate below 1048576
23
24
25 // Triangle Area
26
27 float labs_triarea(const vector pos_a, pos_b, pos_c)
28 {
29 return 0.5 * length(cross(pos_a - pos_c, pos_b - pos_c));
30 }
31
32 float labs_triarea(const int geometry, ptnum_a, ptnum_b, ptnum_c)
33 {
34 return labs_triarea(vector(point(geometry, "P", ptnum_a)),
35 vector(point(geometry, "P", ptnum_b)),
36 vector(point(geometry, "P", ptnum_c)));
37 }
38
39 float labs_triarea(const int geometry, primnum)
40 {
41 return labs_triarea(geometry,
42 primpoint(geometry, primnum, 0),
43 primpoint(geometry, primnum, 1),
44 primpoint(geometry, primnum, 2));
45 }
46
47
48 // Plane Normal
49
50 vector labs_planenormal(const vector pos_a, pos_b, pos_c; const int normalized)
51 {
52 // Uses left-hand rule (Houdini's default face winding order).
53 // If not normalized, the magnitude of the output vector is
54 // the area of the input triangle.
55 return normalized ?
56 normalize(cross(pos_b - pos_c, pos_a - pos_c)) :
57 0.5 * cross(pos_b - pos_c, pos_a - pos_c);
58 }
59
60 vector labs_planenormal(const int geometry, ptnum_a, ptnum_b, ptnum_c, normalized)
61 {
62 return labs_planenormal(vector(point(geometry, "P", ptnum_a)),
63 vector(point(geometry, "P", ptnum_b)),
64 vector(point(geometry, "P", ptnum_c)), normalized);
65 }
66
67 vector labs_planenormal(const int geometry, primnum, normalized)
68 {
69 return labs_planenormal(geometry,
70 primpoint(geometry, primnum, 0),
71 primpoint(geometry, primnum, 1),
72 primpoint(geometry, primnum, 2), normalized);
73 }
74
75
76 // Is to the Right Of
77
78 int labs_istotheright(const vector pivot, forward, up, pos)
79 {
80 return (int)sign(dot(cross(pos - pivot, forward), up));
81 }
82
83
84 // Is Above
85
86 int labs_isabove(const vector plane_pos, plane_normal, pos)
87 {
88 return (int)sign(dot(pos - plane_pos, plane_normal));
89 }
90
91 int labs_isabove(const vector pos_a, pos_b, pos_c, pos)
92 {
93 // Uses left-hand rule (Houdini's default face winding order)
94 return (int)sign(dot(cross(pos - pos_a, pos_c - pos_a), pos_b - pos_a));
95 }
96
97 int labs_isabove(const int geometry, ptnum_a, ptnum_b, ptnum_c; const vector pos)
98 {
99 return labs_isabove(vector(point(geometry, "P", ptnum_a)),
100 vector(point(geometry, "P", ptnum_b)),
101 vector(point(geometry, "P", ptnum_c)), pos);
102 }
103
104 int labs_isabove(const int geometry, primnum; const vector pos)
105 {
106 return labs_isabove(geometry,
107 primpoint(geometry, primnum, 0),
108 primpoint(geometry, primnum, 1),
109 primpoint(geometry, primnum, 2), pos);
110 }
111
112
113 // Point-Plane Distance
114
115 float labs_pointplanedist(const vector plane_pos, plane_unit_normal, pos)
116 {
117 return dot(pos - plane_pos, plane_unit_normal);
118 }
119
120 float labs_pointplanedist(const vector pos_a, pos_b, pos_c, pos)
121 {
122 return labs_pointplanedist(pos_a, labs_planenormal(pos_a, pos_b, pos_c, 1), pos);
123 }
124
125 float labs_pointplanedist(const int geometry, ptnum_a, ptnum_b, ptnum_c; const vector pos)
126 {
127 return labs_pointplanedist(vector(point(geometry, "P", ptnum_a)),
128 vector(point(geometry, "P", ptnum_b)),
129 vector(point(geometry, "P", ptnum_c)), pos);
130 }
131
132 float labs_pointplanedist(const int geometry, primnum; const vector pos)
133 {
134 return labs_pointplanedist(geometry,
135 primpoint(geometry, primnum, 0),
136 primpoint(geometry, primnum, 1),
137 primpoint(geometry, primnum, 2), pos);
138 }
139
140
141 // Point-Plane Projection
142
143 float labs_pointplaneproj(const vector plane_pos, plane_unit_normal, pos; export vector nearest_pos)
144 {
145 float signed_dist = dot(pos - plane_pos, plane_unit_normal);
146 nearest_pos = pos - signed_dist * plane_unit_normal;
147 return signed_dist;
148 }
149
150 float labs_pointplaneproj(const vector pos_a, pos_b, pos_c, pos; export vector nearest_pos)
151 {
152 return labs_pointplaneproj(pos_a, labs_planenormal(pos_a, pos_b, pos_c, 1), pos, nearest_pos);
153 }
154
155 float labs_pointplaneproj(const int geometry, ptnum_a, ptnum_b, ptnum_c; const vector pos; export vector nearest_pos)
156 {
157 return labs_pointplaneproj(vector(point(geometry, "P", ptnum_a)),
158 vector(point(geometry, "P", ptnum_b)),
159 vector(point(geometry, "P", ptnum_c)), pos, nearest_pos);
160 }
161
162 float labs_pointplaneproj(const int geometry, primnum; const vector pos; export vector nearest_pos)
163 {
164 return labs_pointplaneproj(geometry,
165 primpoint(geometry, primnum, 0),
166 primpoint(geometry, primnum, 1),
167 primpoint(geometry, primnum, 2), pos, nearest_pos);
168 }
169
170
171 // Point-Line Distance
172
173 float labs_pointlinedist(const vector line_pos, line_unit_direction, pos)
174 {
175 // This outperforms the cross product (triangle area) approach in accuracy
176 return distance(pos, line_pos + dot(pos - line_pos, line_unit_direction) * line_unit_direction);
177 }
178
179 float labs_pointlinedist(const int geometry, ptnum_a, ptnum_b; const vector pos)
180 {
181 vector line_pos = point(geometry, "P", ptnum_a);
182 return labs_pointlinedist(line_pos, normalize(vector(point(geometry, "P", ptnum_b)) - line_pos), pos);
183 }
184
185 float labs_pointlinedist(const int geometry, primnum; const vector pos)
186 {
187 return labs_pointlinedist(geometry,
188 primpoint(geometry, primnum, 0),
189 primpoint(geometry, primnum, 1), pos);
190 }
191
192
193 // Point-Line Projection
194
195 float labs_pointlineproj(const vector line_pos, line_unit_direction, pos; export vector nearest_pos)
196 {
197 nearest_pos = line_pos + dot(pos - line_pos, line_unit_direction) * line_unit_direction;
198 return distance(pos, nearest_pos);
199 }
200
201 float labs_pointlineproj(const int geometry, ptnum_a, ptnum_b; const vector pos; export vector nearest_pos)
202 {
203 vector line_pos = point(geometry, "P", ptnum_a);
204 return labs_pointlineproj(line_pos, normalize(vector(point(geometry, "P", ptnum_b)) - line_pos), pos, nearest_pos);
205 }
206
207 float labs_pointlineproj(const int geometry, primnum; const vector pos; export vector nearest_pos)
208 {
209 return labs_pointlineproj(geometry,
210 primpoint(geometry, primnum, 0),
211 primpoint(geometry, primnum, 1), pos, nearest_pos);
212 }
213
214
215 // Point-Segment Distance
216
217 float labs_pointsegmentdist(const vector pos_a, pos_b, pos)
218 {
219 vector vec_ab = pos_b - pos_a;
220
221 if (dot(pos - pos_a, vec_ab) <= 0)
222 return distance(pos, pos_a);
223 else if (dot(pos - pos_b, vec_ab) >= 0)
224 return distance(pos, pos_b);
225 else
226 return labs_pointlinedist(pos_a, normalize(vec_ab), pos);
227 }
228
229 float labs_pointsegmentdist(const int geometry, ptnum_a, ptnum_b; const vector pos)
230 {
231 return labs_pointsegmentdist(vector(point(geometry, "P", ptnum_a)),
232 vector(point(geometry, "P", ptnum_b)), pos);
233 }
234
235 float labs_pointsegmentdist(const int geometry, primnum; const vector pos)
236 {
237 return labs_pointsegmentdist(geometry,
238 primpoint(geometry, primnum, 0),
239 primpoint(geometry, primnum, 1), pos);
240 }
241
242
243 // Point-Segment Projection
244
245 float labs_pointsegmentproj(const vector pos_a, pos_b, pos; export vector nearest_pos)
246 {
247 vector vec_ab = pos_b - pos_a;
248
249 if (dot(pos - pos_a, vec_ab) <= 0)
250 {
251 nearest_pos = pos_a;
252 return distance(pos, nearest_pos);
253 }
254 else if (dot(pos - pos_b, vec_ab) >= 0)
255 {
256 nearest_pos = pos_b;
257 return distance(pos, nearest_pos);
258 }
259 else
260 {
261 return labs_pointlineproj(pos_a, normalize(vec_ab), pos, nearest_pos);
262 }
263 }
264
265 float labs_pointsegmentproj(const int geometry, ptnum_a, ptnum_b; const vector pos; export vector nearest_pos)
266 {
267 return labs_pointsegmentproj(vector(point(geometry, "P", ptnum_a)),
268 vector(point(geometry, "P", ptnum_b)), pos, nearest_pos);
269 }
270
271 float labs_pointsegmentproj(const int geometry, primnum; const vector pos; export vector nearest_pos)
272 {
273 return labs_pointsegmentproj(geometry,
274 primpoint(geometry, primnum, 0),
275 primpoint(geometry, primnum, 1), pos, nearest_pos);
276 }
277
278
279 // Is Outside Triangle
280
281 int labs_isoutsidetri(const vector pos_a, pos_b, pos_c, coplanar_pos)
282 {
283 vector vec_a = pos_a - coplanar_pos;
284 vector vec_b = pos_b - coplanar_pos;
285 vector vec_c = pos_c - coplanar_pos;
286 vector cross_ab = cross(vec_a, vec_b);
287 vector cross_bc = cross(vec_b, vec_c);
288 float dot_abbc = dot(cross_ab, cross_bc);
289
290 if (dot_abbc < 0)
291 {
292 return 1;
293 }
294 else
295 {
296 vector cross_ca = cross(vec_c, vec_a);
297 float dot_bcca = dot(cross_bc, cross_ca);
298
299 if (dot_bcca < 0)
300 {
301 return 1;
302 }
303 else
304 {
305 float dot_caab = dot(cross_ca, cross_ab);
306
307 if (dot_caab < 0)
308 return 1;
309 else if (dot_abbc * dot_bcca * dot_caab == 0)
310 return 0;
311 else
312 return -1;
313 }
314 }
315 }
316
317 int labs_isoutsidetri(const int geometry, ptnum_a, ptnum_b, ptnum_c; const vector coplanar_pos)
318 {
319 return labs_isoutsidetri(vector(point(geometry, "P", ptnum_a)),
320 vector(point(geometry, "P", ptnum_b)),
321 vector(point(geometry, "P", ptnum_c)), coplanar_pos);
322 }
323
324 int labs_isoutsidetri(const int geometry, primnum; const vector coplanar_pos)
325 {
326 return labs_isoutsidetri(geometry,
327 primpoint(geometry, primnum, 0),
328 primpoint(geometry, primnum, 1),
329 primpoint(geometry, primnum, 2), coplanar_pos);
330 }
331
332
333 // Point-Triangle-Edge Distance
334
335 float labs_pointtriedgedist(const vector pos_a, pos_b, pos_c, coplanar_pos)
336 {
337 vector vec_ab = pos_b - pos_a;
338 vector vec_bc = pos_c - pos_b;
339 vector normal = cross(vec_ab, vec_bc);
340 vector vec_ap = coplanar_pos - pos_a;
341 vector vec_bp = coplanar_pos - pos_b;
342 float dot_an = dot(cross(vec_ab, vec_ap), normal);
343
344 if (dot_an <= 0) // Outside of or on line AB:
345 {
346 if (dot(vec_ap, vec_ab) <= 0) // Behind AB:
347 {
348 vector vec_ca = pos_a - pos_c;
349
350 if (dot(vec_ap, vec_ca) >= 0) // In front of CA:
351 {
352 return distance(coplanar_pos, pos_a);
353 }
354 else // Between C and A (overlapping zone):
355 {
356 vector vec_ca_norm = normalize(vec_ca);
357 return distance(coplanar_pos, pos_c + max(0, dot(coplanar_pos - pos_c, vec_ca_norm)) * vec_ca_norm);
358 }
359 }
360 else if (dot(vec_bp, vec_ab) >= 0) // In front of AB:
361 {
362 if (dot(vec_bp, vec_bc) <= 0) // Behind BC:
363 {
364 return distance(coplanar_pos, pos_b);
365 }
366 else // Between B and C (overlapping zone):
367 {
368 vector vec_bc_norm = normalize(vec_bc);
369 return distance(coplanar_pos, pos_b + min(length(vec_bc), dot(vec_bp, vec_bc_norm)) * vec_bc_norm);
370 }
371 }
372 else // Between A and B:
373 {
374 vector vec_ab_norm = normalize(vec_ab);
375 return distance(coplanar_pos, pos_a + dot(vec_ap, vec_ab_norm) * vec_ab_norm);
376 }
377 }
378 else
379 {
380 float dot_bn = dot(cross(vec_bc, vec_bp), normal);
381 vector vec_cp = coplanar_pos - pos_c;
382
383 if (dot_bn <= 0) // Outside of or on line BC:
384 {
385 if (dot(vec_bp, vec_bc) <= 0) // Behind BC:
386 {
387 return distance(coplanar_pos, pos_b);
388 }
389 else if (dot(vec_cp, vec_bc) >= 0) // In front of BC:
390 {
391 vector vec_ca = pos_a - pos_c;
392
393 if (dot(vec_cp, vec_ca) <= 0) // Behind CA:
394 {
395 return distance(coplanar_pos, pos_c);
396 }
397 else // Between C and A (overlapping zone):
398 {
399 vector vec_ca_norm = normalize(vec_ca);
400 return distance(coplanar_pos, pos_c + min(length(vec_ca), dot(vec_cp, vec_ca_norm)) * vec_ca_norm);
401 }
402 }
403 else // Between B and C:
404 {
405 vector vec_bc_norm = normalize(vec_bc);
406 return distance(coplanar_pos, pos_b + dot(vec_bp, vec_bc_norm) * vec_bc_norm);
407 }
408 }
409 else
410 {
411 vector vec_ca = pos_a - pos_c;
412 float dot_cn = dot(cross(vec_ca, vec_cp), normal);
413
414 if (dot_cn <= 0) // Outside of or on line CA:
415 {
416 if (dot(vec_cp, vec_ca) <= 0) // Behind CA:
417 {
418 return distance(coplanar_pos, pos_c);
419 }
420 else if (dot(vec_ap, vec_ca) >= 0) // In front of CA:
421 {
422 return distance(coplanar_pos, pos_a);
423 }
424 else // Between C and A:
425 {
426 vector vec_ca_norm = normalize(vec_ca);
427 return distance(coplanar_pos, pos_c + dot(vec_cp, vec_ca_norm) * vec_ca_norm);
428 }
429 }
430 else // Inside the triangle:
431 {
432 vector vec_ab_norm = normalize(vec_ab);
433 vector vec_bc_norm = normalize(vec_bc);
434 vector vec_ca_norm = normalize(vec_ca);
435
436 return -min(distance(coplanar_pos, pos_a + dot(vec_ap, vec_ab_norm) * vec_ab_norm),
437 distance(coplanar_pos, pos_b + dot(vec_bp, vec_bc_norm) * vec_bc_norm),
438 distance(coplanar_pos, pos_c + dot(vec_cp, vec_ca_norm) * vec_ca_norm));
439 }
440 }
441 }
442 }
443
444 float labs_pointtriedgedist(const int geometry, ptnum_a, ptnum_b, ptnum_c; const vector coplanar_pos)
445 {
446 return labs_pointtriedgedist(vector(point(geometry, "P", ptnum_a)),
447 vector(point(geometry, "P", ptnum_b)),
448 vector(point(geometry, "P", ptnum_c)), coplanar_pos);
449 }
450
451 float labs_pointtriedgedist(const int geometry, primnum; const vector coplanar_pos)
452 {
453 return labs_pointtriedgedist(geometry,
454 primpoint(geometry, primnum, 0),
455 primpoint(geometry, primnum, 1),
456 primpoint(geometry, primnum, 2), coplanar_pos);
457 }
458
459
460 // Point-Triangle-Edge Projection
461
462 float labs_pointtriedgeproj(const vector pos_a, pos_b, pos_c, coplanar_pos; export vector nearest_pos)
463 {
464 vector vec_ab = pos_b - pos_a;
465 vector vec_bc = pos_c - pos_b;
466 vector normal = cross(vec_ab, vec_bc);
467 vector vec_ap = coplanar_pos - pos_a;
468 vector vec_bp = coplanar_pos - pos_b;
469 float dot_an = dot(cross(vec_ab, vec_ap), normal);
470
471 if (dot_an <= 0) // Outside of or on line AB:
472 {
473 if (dot(vec_ap, vec_ab) <= 0) // Behind AB:
474 {
475 vector vec_ca = pos_a - pos_c;
476
477 if (dot(vec_ap, vec_ca) >= 0) // In front of CA:
478 {
479 nearest_pos = pos_a;
480 return distance(coplanar_pos, nearest_pos);
481 }
482 else // Between C and A (overlapping zone):
483 {
484 vector vec_ca_norm = normalize(vec_ca);
485 nearest_pos = pos_c + max(0, dot(coplanar_pos - pos_c, vec_ca_norm)) * vec_ca_norm;
486 return distance(coplanar_pos, nearest_pos);
487 }
488 }
489 else if (dot(vec_bp, vec_ab) >= 0) // In front of AB:
490 {
491 if (dot(vec_bp, vec_bc) <= 0) // Behind BC:
492 {
493 nearest_pos = pos_b;
494 return distance(coplanar_pos, nearest_pos);
495 }
496 else // Between B and C (overlapping zone):
497 {
498 vector vec_bc_norm = normalize(vec_bc);
499 nearest_pos = pos_b + min(length(vec_bc), dot(vec_bp, vec_bc_norm)) * vec_bc_norm;
500 return distance(coplanar_pos, nearest_pos);
501 }
502 }
503 else // Between A and B:
504 {
505 vector vec_ab_norm = normalize(vec_ab);
506 nearest_pos = pos_a + dot(vec_ap, vec_ab_norm) * vec_ab_norm;
507 return distance(coplanar_pos, nearest_pos);
508 }
509 }
510 else
511 {
512 float dot_bn = dot(cross(vec_bc, vec_bp), normal);
513 vector vec_cp = coplanar_pos - pos_c;
514
515 if (dot_bn <= 0) // Outside of or on line BC:
516 {
517 if (dot(vec_bp, vec_bc) <= 0) // Behind BC:
518 {
519 nearest_pos = pos_b;
520 return distance(coplanar_pos, nearest_pos);
521 }
522 else if (dot(vec_cp, vec_bc) >= 0) // In front of BC:
523 {
524 vector vec_ca = pos_a - pos_c;
525
526 if (dot(vec_cp, vec_ca) <= 0) // Behind CA:
527 {
528 nearest_pos = pos_c;
529 return distance(coplanar_pos, nearest_pos);
530 }
531 else // Between C and A (overlapping zone):
532 {
533 vector vec_ca_norm = normalize(vec_ca);
534 nearest_pos = pos_c + min(length(vec_ca), dot(vec_cp, vec_ca_norm)) * vec_ca_norm;
535 return distance(coplanar_pos, nearest_pos);
536 }
537 }
538 else // Between B and C:
539 {
540 vector vec_bc_norm = normalize(vec_bc);
541 nearest_pos = pos_b + dot(vec_bp, vec_bc_norm) * vec_bc_norm;
542 return distance(coplanar_pos, nearest_pos);
543 }
544 }
545 else
546 {
547 vector vec_ca = pos_a - pos_c;
548 float dot_cn = dot(cross(vec_ca, vec_cp), normal);
549
550 if (dot_cn <= 0) // Outside of or on line CA:
551 {
552 if (dot(vec_cp, vec_ca) <= 0) // Behind CA:
553 {
554 nearest_pos = pos_c;
555 return distance(coplanar_pos, nearest_pos);
556 }
557 else if (dot(vec_ap, vec_ca) >= 0) // In front of CA:
558 {
559 nearest_pos = pos_a;
560 return distance(coplanar_pos, nearest_pos);
561 }
562 else // Between C and A:
563 {
564 vector vec_ca_norm = normalize(vec_ca);
565 nearest_pos = pos_c + dot(vec_cp, vec_ca_norm) * vec_ca_norm;
566 return distance(coplanar_pos, nearest_pos);
567 }
568 }
569 else // Inside the triangle:
570 {
571 vector vec_ab_norm = normalize(vec_ab);
572 vector vec_bc_norm = normalize(vec_bc);
573 vector vec_ca_norm = normalize(vec_ca);
574
575 vector nearest_pos_ab = pos_a + dot(vec_ap, vec_ab_norm) * vec_ab_norm;
576 vector nearest_pos_bc = pos_b + dot(vec_bp, vec_bc_norm) * vec_bc_norm;
577 vector nearest_pos_ca = pos_c + dot(vec_cp, vec_ca_norm) * vec_ca_norm;
578
579 float dist_ab = distance(coplanar_pos, nearest_pos_ab);
580 float dist_bc = distance(coplanar_pos, nearest_pos_bc);
581 float dist_ca = distance(coplanar_pos, nearest_pos_ca);
582
583 if (dist_ab <= dist_bc && dist_ab <= dist_ca)
584 {
585 nearest_pos = nearest_pos_ab;
586 return -dist_ab;
587 }
588 else if (dist_bc <= dist_ab && dist_bc <= dist_ca)
589 {
590 nearest_pos = nearest_pos_bc;
591 return -dist_bc;
592 }
593 else
594 {
595 nearest_pos = nearest_pos_ca;
596 return -dist_ca;
597 }
598 }
599 }
600 }
601 }
602
603 float labs_pointtriedgeproj(const int geometry, ptnum_a, ptnum_b, ptnum_c; const vector coplanar_pos; export vector nearest_pos)
604 {
605 return labs_pointtriedgeproj(vector(point(geometry, "P", ptnum_a)),
606 vector(point(geometry, "P", ptnum_b)),
607 vector(point(geometry, "P", ptnum_c)), coplanar_pos, nearest_pos);
608 }
609
610 float labs_pointtriedgeproj(const int geometry, primnum; const vector coplanar_pos; export vector nearest_pos)
611 {
612 return labs_pointtriedgeproj(geometry,
613 primpoint(geometry, primnum, 0),
614 primpoint(geometry, primnum, 1),
615 primpoint(geometry, primnum, 2), coplanar_pos, nearest_pos);
616 }
617
618
619 // Point-Triangle Distance
620
621 float labs_pointtridist(const vector pos_a, pos_b, pos_c, pos)
622 {
623 vector coplanar_pos;
624 float signed_pt_plane_dist = labs_pointplaneproj(pos_a, pos_b, pos_c, pos, coplanar_pos);
625
626 if (labs_isoutsidetri(pos_a, pos_b, pos_c, coplanar_pos) <= 0) // Coplanar position is not outside triangle
627 {
628 return abs(signed_pt_plane_dist);
629 }
630 else // Coplanar position is outside triangle
631 {
632 float pt_tri_edge_dist = labs_pointtriedgedist(pos_a, pos_b, pos_c, coplanar_pos);
633 return sqrt(signed_pt_plane_dist * signed_pt_plane_dist + pt_tri_edge_dist * pt_tri_edge_dist);
634 }
635 }
636
637 float labs_pointtridist(const int geometry; const int ptnum_a, ptnum_b, ptnum_c; const vector pos)
638 {
639 return labs_pointtridist(vector(point(geometry, "P", ptnum_a)),
640 vector(point(geometry, "P", ptnum_b)),
641 vector(point(geometry, "P", ptnum_c)), pos);
642 }
643
644 float labs_pointtridist(const int geometry, primnum; const vector pos)
645 {
646 return labs_pointtridist(geometry,
647 primpoint(geometry, primnum, 0),
648 primpoint(geometry, primnum, 1),
649 primpoint(geometry, primnum, 2), pos);
650 }
651
652
653 // Point-Triangle Projection
654
655 float labs_pointtriproj(const vector pos_a, pos_b, pos_c, pos; export vector nearest_pos)
656 {
657 vector coplanar_pos;
658 float signed_pt_plane_dist = labs_pointplaneproj(pos_a, pos_b, pos_c, pos, coplanar_pos);
659
660 if (labs_isoutsidetri(pos_a, pos_b, pos_c, coplanar_pos) <= 0) // Coplanar position is not outside triangle
661 {
662 nearest_pos = coplanar_pos;
663 return abs(signed_pt_plane_dist);
664 }
665 else // Coplanar position is outside triangle
666 {
667 float pt_tri_edge_dist = labs_pointtriedgeproj(pos_a, pos_b, pos_c, coplanar_pos, nearest_pos);
668 return sqrt(signed_pt_plane_dist * signed_pt_plane_dist + pt_tri_edge_dist * pt_tri_edge_dist);
669 }
670 }
671
672 float labs_pointtriproj(const int geometry; const int ptnum_a, ptnum_b, ptnum_c; const vector pos; export vector nearest_pos)
673 {
674 return labs_pointtriproj(vector(point(geometry, "P", ptnum_a)),
675 vector(point(geometry, "P", ptnum_b)),
676 vector(point(geometry, "P", ptnum_c)), pos, nearest_pos);
677 }
678
679 float labs_pointtriproj(const int geometry, primnum; const vector pos; export vector nearest_pos)
680 {
681 return labs_pointtriproj(geometry,
682 primpoint(geometry, primnum, 0),
683 primpoint(geometry, primnum, 1),
684 primpoint(geometry, primnum, 2), pos, nearest_pos);
685 }
686
687
688 // Triangle Circumcenter
689
690 vector labs_circumcenter_tri(const vector pos_a, pos_b, pos_c; export float radius)
691 {
692 vector vec_ab = pos_b - pos_a;
693 vector vec_ac = pos_c - pos_a;
694 vector cross_abac = cross(vec_ab, vec_ac);
695
696 vector a_to_center = (dot(vec_ac, vec_ac) * cross(cross_abac, vec_ab) + dot(vec_ab, vec_ab) * cross(vec_ac, cross_abac)) /
697 (2.0 * dot(cross_abac, cross_abac));
698
699 radius = length(a_to_center);
700 return pos_a + a_to_center;
701 }
702
703
704 vector labs_circumcenter_tri(const int geometry, ptnum_a, ptnum_b, ptnum_c; export float radius)
705 {
706 return labs_circumcenter_tri(vector(point(geometry, "P", ptnum_a)),
707 vector(point(geometry, "P", ptnum_b)),
708 vector(point(geometry, "P", ptnum_c)), radius);
709 }
710
711 vector labs_circumcenter_tri(const int geometry, primnum; export float radius)
712 {
713 return labs_circumcenter_tri(geometry,
714 primpoint(geometry, primnum, 0),
715 primpoint(geometry, primnum, 1),
716 primpoint(geometry, primnum, 2), radius);
717 }
718
719
720 // Triangle Bounding Sphere
721
722 vector labs_boundsphere_tri(const vector pos_a, pos_b, pos_c; export float radius)
723 {
724 float dist_ab = distance(pos_a, pos_b);
725 float dist_bc = distance(pos_b, pos_c);
726 float dist_ca = distance(pos_c, pos_a);
727
728 if (dist_ab >= dist_bc && dist_ab >= dist_ca)
729 {
730 radius = 0.5 * dist_ab;
731 return 0.5 * (pos_a + pos_b);
732 }
733 else if (dist_bc >= dist_ab && dist_bc >= dist_ca)
734 {
735 radius = 0.5 * dist_bc;
736 return 0.5 * (pos_b + pos_c);
737 }
738 else
739 {
740 radius = 0.5 * dist_ca;
741 return 0.5 * (pos_c + pos_a);
742 }
743 }
744
745
746 vector labs_boundsphere_tri(const int geometry, ptnum_a, ptnum_b, ptnum_c; export float radius)
747 {
748 return labs_boundsphere_tri(vector(point(geometry, "P", ptnum_a)),
749 vector(point(geometry, "P", ptnum_b)),
750 vector(point(geometry, "P", ptnum_c)), radius);
751 }
752
753 vector labs_boundsphere_tri(const int geometry, primnum; export float radius)
754 {
755 return labs_boundsphere_tri(geometry,
756 primpoint(geometry, primnum, 0),
757 primpoint(geometry, primnum, 1),
758 primpoint(geometry, primnum, 2), radius);
759 }
760
761
762 // Angle between Directions
763
764 float labs_anglebetween(const vector unit_vec_a, unit_vec_b; const int in_degrees)
765 {
766 return in_degrees ?
767 degrees(acos(dot(unit_vec_a, unit_vec_b))) :
768 acos(dot(unit_vec_a, unit_vec_b));
769 }
770
771
772 // Rotate Vector
773
774 vector labs_rotatevector(const vector vec; const float angle; const vector unit_axis; const int in_degrees)
775 {
776 return qrotate(quaternion(in_degrees ? radians(angle) : angle, unit_axis), vec);
777 }
778
779 vector labs_rotatevector(const vector vec, start_vec, end_vec)
780 {
781 return qrotate(dihedral(start_vec, end_vec), vec);
782 }
783
784
785 // Rotate Vector 2D
786
787 vector2 labs_rotatevector2d(const vector2 vec2d; const float angle; const int in_degrees)
788 {
789 float angle_rad = in_degrees ? radians(angle) : angle;
790 float cos_angle = cos(angle_rad);
791 float sin_angle = sin(angle_rad);
792
793 return set(cos_angle * vec2d.x - sin_angle * vec2d.y,
794 sin_angle * vec2d.x + cos_angle * vec2d.y);
795 }
796
797 vector2 labs_rotatevector2d(const vector2 vec2d, start_vec2d, end_vec2d)
798 {
799 return (vector2)labs_rotatevector((vector)vec2d, (vector)start_vec2d, (vector)end_vec2d);
800 }
801
802
803 // Get Non-collinear
804
805 vector labs_getnoncollinear(const vector vec)
806 {
807 return set(vec.z, vec.x, -vec.y);
808 }
809
810
811 #endif