JOGL v2.6.0-rc-20250712
JOGL, High-Performance Graphics Binding for Java™ (public API).
VectorUtil.java
Go to the documentation of this file.
1/**
2 * Copyright 2010-2024 JogAmp Community. All rights reserved.
3 *
4 * Redistribution and use in source and binary forms, with or without modification, are
5 * permitted provided that the following conditions are met:
6 *
7 * 1. Redistributions of source code must retain the above copyright notice, this list of
8 * conditions and the following disclaimer.
9 *
10 * 2. Redistributions in binary form must reproduce the above copyright notice, this list
11 * of conditions and the following disclaimer in the documentation and/or other materials
12 * provided with the distribution.
13 *
14 * THIS SOFTWARE IS PROVIDED BY JogAmp Community ``AS IS'' AND ANY EXPRESS OR IMPLIED
15 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
16 * FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL JogAmp Community OR
17 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
18 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
19 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
20 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
21 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
22 * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
23 *
24 * The views and conclusions contained in the software and documentation are those of the
25 * authors and should not be interpreted as representing official policies, either expressed
26 * or implied, of JogAmp Community.
27 */
28package com.jogamp.math;
29
30import java.util.List;
31
32import com.jogamp.graph.geom.Vertex;
33import com.jogamp.math.geom.plane.Winding;
34
35public final class VectorUtil {
36 /**
37 * Return true if 2D vector components are zero, no {@link FloatUtil#EPSILON} is taken into consideration.
38 */
39 public static boolean isVec2Zero(final Vec3f vec) {
40 return 0f == vec.x() && 0f == vec.y();
41 }
42
43 /**
44 * Return true if all three vector components are zero, i.e. it's their absolute value < <code>epsilon</code>.
45 * <p>
46 * Implementation uses {@link FloatUtil#isZero(float, float)}, see API doc for details.
47 * </p>
48 */
49 public static boolean isZero(final float x, final float y, final float z, final float epsilon) {
50 return FloatUtil.isZero(x, epsilon) &&
51 FloatUtil.isZero(y, epsilon) &&
52 FloatUtil.isZero(z, epsilon) ;
53 }
54
55 /**
56 * Return true if all three vector components are zero, i.e. it's their absolute value < {@link FloatUtil#EPSILON}.
57 * <p>
58 * Implementation uses {@link FloatUtil#isZero(float)}, see API doc for details.
59 * </p>
60 */
61 public static boolean isZero(final float x, final float y, final float z) {
62 return FloatUtil.isZero(x) &&
63 FloatUtil.isZero(y) &&
64 FloatUtil.isZero(z) ;
65 }
66
67 /**
68 * Return the squared distance between the given two points described vector v1 and v2.
69 * <p>
70 * When comparing the relative distance between two points it is usually sufficient to compare the squared
71 * distances, thus avoiding an expensive square root operation.
72 * </p>
73 */
74 public static float distSquareVec3(final float[] v1, final float[] v2) {
75 final float dx = v1[0] - v2[0];
76 final float dy = v1[1] - v2[1];
77 final float dz = v1[2] - v2[2];
78 return dx * dx + dy * dy + dz * dz;
79 }
80
81 /**
82 * Return the distance between the given two points described vector v1 and v2.
83 */
84 public static float distVec3(final float[] v1, final float[] v2) {
85 return FloatUtil.sqrt(distSquareVec3(v1, v2));
86 }
87
88 /**
89 * Return the squared length of a vector, a.k.a the squared <i>norm</i> or squared <i>magnitude</i>
90 */
91 public static float normSquareVec2(final float[] vec) {
92 return vec[0]*vec[0] + vec[1]*vec[1];
93 }
94
95 /**
96 * Return the squared length of a vector, a.k.a the squared <i>norm</i> or squared <i>magnitude</i>
97 */
98 public static float normSquareVec3(final float[] vec) {
99 return vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2];
100 }
101
102 /**
103 * Return the squared length of a vector, a.k.a the squared <i>norm</i> or squared <i>magnitude</i>
104 */
105 public static float normSquareVec3(final float[] vec, final int offset) {
106 float v = vec[0+offset];
107 float r = v*v;
108 v = vec[1+offset];
109 r += v*v;
110 v = vec[2+offset];
111 return r + v*v;
112 }
113
114 /**
115 * Return the length of a vector, a.k.a the <i>norm</i> or <i>magnitude</i>
116 */
117 public static float normVec2(final float[] vec) {
118 return FloatUtil.sqrt(normSquareVec2(vec));
119 }
120
121 /**
122 * Normalize a vector in place
123 * @param vector input vector
124 * @return normalized output vector
125 */
126 public static float[] normalizeVec3(final float[] vector) {
127 final float lengthSq = normSquareVec3(vector);
128 if ( FloatUtil.isZero(lengthSq) ) {
129 vector[0] = 0f;
130 vector[1] = 0f;
131 vector[2] = 0f;
132 } else {
133 final float invSqr = 1f / FloatUtil.sqrt(lengthSq);
134 vector[0] *= invSqr;
135 vector[1] *= invSqr;
136 vector[2] *= invSqr;
137 }
138 return vector;
139 }
140
141 /**
142 * Normalize a vector in place
143 * @param vector input vector
144 * @return normalized output vector
145 */
146 public static float[] normalizeVec3(final float[] vector, final int offset) {
147 final float lengthSq = normSquareVec3(vector, offset);
148 if ( FloatUtil.isZero(lengthSq) ) {
149 vector[0+offset] = 0f;
150 vector[1+offset] = 0f;
151 vector[2+offset] = 0f;
152 } else {
153 final float invSqr = 1f / FloatUtil.sqrt(lengthSq);
154 vector[0+offset] *= invSqr;
155 vector[1+offset] *= invSqr;
156 vector[2+offset] *= invSqr;
157 }
158 return vector;
159 }
160
161 /**
162 * Scales a vector by param using given result float[], result = vector * scale
163 * @param result vector for the result, may be vector (in-place)
164 * @param vector input vector
165 * @param scale single scale constant for all vector components
166 * @return result vector for chaining
167 */
168 public static float[] scaleVec2(final float[] result, final float[] vector, final float scale) {
169 result[0] = vector[0] * scale;
170 result[1] = vector[1] * scale;
171 return result;
172 }
173
174 /**
175 * Scales a vector by param using given result float[], result = vector * scale
176 * @param result vector for the result, may be vector (in-place)
177 * @param vector input vector
178 * @param scale 2 component scale constant for each vector component
179 * @return result vector for chaining
180 */
181 public static float[] scaleVec2(final float[] result, final float[] vector, final float[] scale)
182 {
183 result[0] = vector[0] * scale[0];
184 result[1] = vector[1] * scale[1];
185 return result;
186 }
187
188 /**
189 * Divides a vector by param using given result float[], result = vector / scale
190 * @param result vector for the result, may be vector (in-place)
191 * @param vector input vector
192 * @param scale single scale constant for all vector components
193 * @return result vector for chaining
194 */
195 public static float[] divVec2(final float[] result, final float[] vector, final float scale) {
196 result[0] = vector[0] / scale;
197 result[1] = vector[1] / scale;
198 return result;
199 }
200
201 /**
202 * Divides a vector by param using given result float[], result = vector / scale
203 * @param result vector for the result, may be vector (in-place)
204 * @param vector input vector
205 * @param scale 2 component scale constant for each vector component
206 * @return result vector for chaining
207 */
208 public static float[] divVec2(final float[] result, final float[] vector, final float[] scale)
209 {
210 result[0] = vector[0] / scale[0];
211 result[1] = vector[1] / scale[1];
212 return result;
213 }
214
215 /**
216 * Adds two vectors, result = v1 + v2
217 * @param result float[2] result vector, may be either v1 or v2 (in-place)
218 * @param v1 vector 1
219 * @param v2 vector 2
220 * @return result vector for chaining
221 */
222 public static float[] addVec2(final float[] result, final float[] v1, final float[] v2) {
223 result[0] = v1[0] + v2[0];
224 result[1] = v1[1] + v2[1];
225 return result;
226 }
227
228 /**
229 * Subtracts two vectors, result = v1 - v2
230 * @param result float[2] result vector, may be either v1 or v2 (in-place)
231 * @param v1 vector 1
232 * @param v2 vector 2
233 * @return result vector for chaining
234 */
235 public static float[] subVec2(final float[] result, final float[] v1, final float[] v2) {
236 result[0] = v1[0] - v2[0];
237 result[1] = v1[1] - v2[1];
238 return result;
239 }
240
241 /**
242 * cross product vec1 x vec2
243 * @param v1 vector 1
244 * @param v2 vector 2
245 * @return the resulting vector
246 */
247 public static float[] crossVec3(final float[] r, final int r_offset, final float[] v1, final int v1_offset, final float[] v2, final int v2_offset)
248 {
249 r[0+r_offset] = v1[1+v1_offset] * v2[2+v2_offset] - v1[2+v1_offset] * v2[1+v2_offset];
250 r[1+r_offset] = v1[2+v1_offset] * v2[0+v2_offset] - v1[0+v1_offset] * v2[2+v2_offset];
251 r[2+r_offset] = v1[0+v1_offset] * v2[1+v2_offset] - v1[1+v1_offset] * v2[0+v2_offset];
252 return r;
253 }
254
255 /**
256 * Calculate the midpoint of two points
257 * @param p1 first point vector
258 * @param p2 second point vector
259 * @return midpoint
260 */
261 public static Vec3f midpoint(final Vec3f result, final Vec3f p1, final Vec3f p2) {
262 result.set( (p1.x() + p2.x())*0.5f,
263 (p1.y() + p2.y())*0.5f,
264 (p1.z() + p2.z())*0.5f );
265 return result;
266 }
267
268 /**
269 * Return the determinant of 3 vectors
270 * @param a vector 1
271 * @param b vector 2
272 * @param c vector 3
273 * @return the determinant value
274 */
275 public static float determinant(final Vec3f a, final Vec3f b, final Vec3f c) {
276 return a.x()*b.y()*c.z() + a.y()*b.z()*c.x() + a.z()*b.x()*c.y() - a.x()*b.z()*c.y() - a.y()*b.x()*c.z() - a.z()*b.y()*c.x();
277 }
278
279 /**
280 * Check if three vertices are collinear
281 * @param v1 vertex 1
282 * @param v2 vertex 2
283 * @param v3 vertex 3
284 * @return true if collinear, false otherwise
285 */
286 public static boolean isCollinear(final Vec3f v1, final Vec3f v2, final Vec3f v3) {
287 return FloatUtil.isZero( determinant(v1, v2, v3) );
288 }
289
290 /**
291 * Check if vertices in triangle circumcircle given {@code d} vertex, from paper by Guibas and Stolfi (1985).
292 * <p>
293 * Implementation uses double precision.
294 * </p>
295 * @param a triangle vertex 1
296 * @param b triangle vertex 2
297 * @param c triangle vertex 3
298 * @param d vertex in question
299 * @return true if the vertex d is inside the circle defined by the vertices a, b, c.
300 */
301 public static boolean isInCircle(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c, final Vert2fImmutable d) {
302 return inCircleVal(a, b, c, d) > DoubleUtil.EPSILON;
303 }
304 public static double inCircleVal(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c, final Vert2fImmutable d) {
305 // Operation costs:
306 // - 4x (triAreaVec2: 5+, 2*) -> 20+, 8*
307 // - plus 7+, 12* -> 27+, 20*
308 return sqlend(a.x(), a.y()) * triArea(b, c, d) -
309 sqlend(b.x(), b.y()) * triArea(a, c, d) +
310 sqlend(c.x(), c.y()) * triArea(a, b, d) -
311 sqlend(d.x(), d.y()) * triArea(a, b, c);
312 }
313 private static double sqlend(final double x, final double y) {
314 return x*x + y*y;
315 }
316
317 /**
318 * Computes oriented double area of a triangle,
319 * i.e. the 2x2 determinant with b-a and c-a per column.
320 * <pre>
321 * | bx-ax, cx-ax |
322 * det = | by-ay, cy-ay |
323 * </pre>
324 * @param a first vertex
325 * @param b second vertex
326 * @param c third vertex
327 * @return area > 0 CCW, ..
328 */
329 public static double triArea(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c){
330 return triArea(a.x(), a.y(), b.x(), b.y(), c.x(), c.y());
331 }
332 private static double triArea(final double ax, final double ay, final double bx, final double by, final double cx, final double cy){
333 return (bx - ax) * (cy - ay) - (by - ay) * (cx - ax);
334 }
335
336 /**
337 * Check if a vertex is in triangle using barycentric coordinates computation.
338 * @param a first triangle vertex
339 * @param b second triangle vertex
340 * @param c third triangle vertex
341 * @param p the vertex in question
342 * @param ac temporary storage
343 * @param ab temporary storage
344 * @param ap temporary storage
345 * @return true if p is in triangle (a, b, c), false otherwise.
346 */
347 public static boolean isInTriangle(final Vec3f a, final Vec3f b, final Vec3f c,
348 final Vec3f p,
349 final Vec3f ac, final Vec3f ab, final Vec3f ap){
350 // Compute vectors
351 ac.minus( c, a); // v0
352 ab.minus( b, a); // v1
353 ap.minus( p, a); // v2
354
355 // Compute dot products
356 final float dotAC_AC = ac.dot(ac);
357 final float dotAC_AB = ac.dot(ab);
358 final float dotAB_AB = ab.dot(ab);
359 final float dotAC_AP = ac.dot(ap);
360 final float dotAB_AP = ab.dot(ap);
361
362 // Compute barycentric coordinates
363 final float invDenom = 1 / (dotAC_AC * dotAB_AB - dotAC_AB * dotAC_AB);
364 final float u = (dotAB_AB * dotAC_AP - dotAC_AB * dotAB_AP) * invDenom;
365 final float v = (dotAC_AC * dotAB_AP - dotAC_AB * dotAC_AP) * invDenom;
366
367 // Check if point is in triangle
368 return (u >= 0) && (v >= 0) && (u + v < 1);
369 }
370
371 /**
372 * Check if one of three vertices are in triangle using barycentric coordinates computation.
373 * @param a first triangle vertex
374 * @param b second triangle vertex
375 * @param c third triangle vertex
376 * @param p1 the vertex in question
377 * @param p2 the vertex in question
378 * @param p3 the vertex in question
379 * @param ac temporary storage
380 * @param ab temporary storage
381 * @param ap temporary storage
382 * @return true if p1 or p2 or p3 is in triangle (a, b, c), false otherwise.
383 */
384 public static boolean isInTriangle3(final Vec3f a, final Vec3f b, final Vec3f c,
385 final Vec3f p1, final Vec3f p2, final Vec3f p3,
386 final Vec3f ac, final Vec3f ab, final Vec3f ap){
387 // Compute vectors
388 ac.minus(c, a); // v0
389 ab.minus(b, a); // v1
390
391 // Compute dot products
392 final float dotAC_AC = ac.dot(ac);
393 final float dotAC_AB = ac.dot(ab);
394 final float dotAB_AB = ab.dot(ab);
395
396 // Compute barycentric coordinates
397 final float invDenom = 1 / (dotAC_AC * dotAB_AB - dotAC_AB * dotAC_AB);
398 {
399 ap.minus(p1, a); // v2
400 final float dotAC_AP1 = ac.dot(ap);
401 final float dotAB_AP1 = ab.dot(ap);
402 final float u = (dotAB_AB * dotAC_AP1 - dotAC_AB * dotAB_AP1) * invDenom;
403 final float v = (dotAC_AC * dotAB_AP1 - dotAC_AB * dotAC_AP1) * invDenom;
404
405 // Check if point is in triangle
406 if ( (u >= 0) && (v >= 0) && (u + v < 1) ) {
407 return true;
408 }
409 }
410
411 {
412 ap.minus(p2, a); // v2
413 final float dotAC_AP2 = ac.dot(ap);
414 final float dotAB_AP2 = ab.dot(ap);
415 final float u = (dotAB_AB * dotAC_AP2 - dotAC_AB * dotAB_AP2) * invDenom;
416 final float v = (dotAC_AC * dotAB_AP2 - dotAC_AB * dotAC_AP2) * invDenom;
417
418 // Check if point is in triangle
419 if ( (u >= 0) && (v >= 0) && (u + v < 1) ) {
420 return true;
421 }
422 }
423
424 {
425 ap.minus(p3, a); // v3
426 final float dotAC_AP3 = ac.dot(ap);
427 final float dotAB_AP3 = ab.dot(ap);
428 final float u = (dotAB_AB * dotAC_AP3 - dotAC_AB * dotAB_AP3) * invDenom;
429 final float v = (dotAC_AC * dotAB_AP3 - dotAC_AB * dotAC_AP3) * invDenom;
430
431 // Check if point is in triangle
432 if ( (u >= 0) && (v >= 0) && (u + v < 1) ) {
433 return true;
434 }
435 }
436 return false;
437 }
438 /**
439 * Check if one of three vertices are in triangle using
440 * barycentric coordinates computation, using given epsilon for comparison.
441 * @param a first triangle vertex
442 * @param b second triangle vertex
443 * @param c third triangle vertex
444 * @param p1 the vertex in question
445 * @param p2 the vertex in question
446 * @param p3 the vertex in question
447 * @param tmpAC
448 * @param tmpAB
449 * @param tmpAP
450 * @return true if p1 or p2 or p3 is in triangle (a, b, c), false otherwise.
451 */
452 public static boolean isInTriangle3(final Vec3f a, final Vec3f b, final Vec3f c,
453 final Vec3f p1, final Vec3f p2, final Vec3f p3,
454 final Vec3f ac, final Vec3f ab, final Vec3f ap,
455 final float epsilon) {
456 // Compute vectors
457 ac.minus(c, a); // v0
458 ab.minus(b, a); // v1
459
460 // Compute dot products
461 final float dotAC_AC = ac.dot(ac);
462 final float dotAC_AB = ac.dot(ab);
463 final float dotAB_AB = ab.dot(ab);
464
465 // Compute barycentric coordinates
466 final float invDenom = 1 / (dotAC_AC * dotAB_AB - dotAC_AB * dotAC_AB);
467 {
468 ap.minus(p1, a); // v2
469 final float dotAC_AP1 = ac.dot(ap);
470 final float dotAB_AP1 = ab.dot(ap);
471 final float u = (dotAB_AB * dotAC_AP1 - dotAC_AB * dotAB_AP1) * invDenom;
472 final float v = (dotAC_AC * dotAB_AP1 - dotAC_AB * dotAC_AP1) * invDenom;
473
474 // Check if point is in triangle
475 if( FloatUtil.compare(u, 0.0f, epsilon) >= 0 &&
476 FloatUtil.compare(v, 0.0f, epsilon) >= 0 &&
477 FloatUtil.compare(u+v, 1.0f, epsilon) < 0 ) {
478 return true;
479 }
480 }
481
482 {
483 ap.minus(p2, a); // v3
484 final float dotAC_AP2 = ac.dot(ap);
485 final float dotAB_AP2 = ab.dot(ap);
486 final float u = (dotAB_AB * dotAC_AP2 - dotAC_AB * dotAB_AP2) * invDenom;
487 final float v = (dotAC_AC * dotAB_AP2 - dotAC_AB * dotAC_AP2) * invDenom;
488
489 // Check if point is in triangle
490 if( FloatUtil.compare(u, 0.0f, epsilon) >= 0 &&
491 FloatUtil.compare(v, 0.0f, epsilon) >= 0 &&
492 FloatUtil.compare(u+v, 1.0f, epsilon) < 0 ) {
493 return true;
494 }
495 }
496
497 {
498 ap.minus(p3, a); // v4
499 final float dotAC_AP3 = ac.dot(ap);
500 final float dotAB_AP3 = ab.dot(ap);
501 final float u = (dotAB_AB * dotAC_AP3 - dotAC_AB * dotAB_AP3) * invDenom;
502 final float v = (dotAC_AC * dotAB_AP3 - dotAC_AB * dotAC_AP3) * invDenom;
503
504 // Check if point is in triangle
505 if( FloatUtil.compare(u, 0.0f, epsilon) >= 0 &&
506 FloatUtil.compare(v, 0.0f, epsilon) >= 0 &&
507 FloatUtil.compare(u+v, 1.0f, epsilon) < 0 ) {
508 return true;
509 }
510 }
511 return false;
512 }
513
514 /**
515 * Check if points are in ccw order
516 * <p>
517 * Consider using {@link #getWinding(List)} using the {@link #area(List)} function over all points
518 * on complex shapes for a reliable result!
519 * </p>
520 * @param a first vertex
521 * @param b second vertex
522 * @param c third vertex
523 * @return true if the points a,b,c are in a ccw order
524 * @see #getWinding(List)
525 */
526 public static boolean isCCW(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c){
527 return triArea(a,b,c) > DoubleUtil.EPSILON;
528 }
529
530 /**
531 * Compute the winding of the 3 given points
532 * <p>
533 * Consider using {@link #getWinding(List)} using the {@link #area(List)} function over all points
534 * on complex shapes for a reliable result!
535 * </p>
536 * @param a first vertex
537 * @param b second vertex
538 * @param c third vertex
539 * @return {@link Winding#CCW} or {@link Winding#CW}
540 * @see #getWinding(List)
541 */
542 public static Winding getWinding(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c) {
543 return triArea(a,b,c) > DoubleUtil.EPSILON ? Winding.CCW : Winding.CW ;
544 }
545
546 /**
547 * Computes the area of a list of vertices via shoelace formula.
548 * <p>
549 * This method is utilized e.g. to reliably compute the {@link Winding} of complex shapes.
550 * </p>
551 * <p>
552 * Implementation uses double precision.
553 * </p>
554 * @param vertices
555 * @return positive area if ccw else negative area value
556 * @see #getWinding(List)
557 */
558 public static double area(final List<? extends Vert2fImmutable> vertices) {
559 final int n = vertices.size();
560 double area = 0.0;
561 for (int p = n - 1, q = 0; q < n; p = q++) {
562 final Vert2fImmutable pCoord = vertices.get(p);
563 final Vert2fImmutable qCoord = vertices.get(q);
564 area += (double)pCoord.x() * (double)qCoord.y() - (double)qCoord.x() * (double)pCoord.y();
565 }
566 return area;
567 }
568
569 /**
570 * Compute the winding using the {@link #area(List)} function over all vertices for complex shapes.
571 * <p>
572 * Uses the {@link #area(List)} function over all points
573 * on complex shapes for a reliable result!
574 * </p>
575 * <p>
576 * Implementation uses double precision.
577 * </p>
578 * @param vertices array of Vertices
579 * @return {@link Winding#CCW} or {@link Winding#CW}
580 * @see #area(List)
581 */
582 public static Winding getWinding(final List<? extends Vert2fImmutable> vertices) {
583 return area(vertices) >= 0 ? Winding.CCW : Winding.CW ;
584 }
585
586 /**
587 * Finds the plane equation of a plane given its normal and a point on the plane.
588 *
589 * @param resultV4 vec4 plane equation
590 * @param normalVec3
591 * @param pVec3
592 * @return result for chaining
593 */
594 public static Vec4f getPlaneVec3(final Vec4f resultV4, final Vec3f normalVec3, final Vec3f pVec3) {
595 /**
596 Ax + By + Cz + D == 0 ;
597 D = - ( Ax + By + Cz )
598 = - ( A*a[0] + B*a[1] + C*a[2] )
599 = - vec3Dot ( normal, a ) ;
600 */
601 resultV4.set(normalVec3, -normalVec3.dot(pVec3));
602 return resultV4;
603 }
604
605 /**
606 * This finds the plane equation of a triangle given three vertices.
607 *
608 * @param resultVec4 vec4 plane equation
609 * @param v1 vec3
610 * @param v2 vec3
611 * @param v3 vec3
612 * @param temp1V3
613 * @param temp2V3
614 * @return result for chaining
615 */
616 public static Vec4f getPlaneVec3(final Vec4f resultVec4, final Vec3f v1, final Vec3f v2, final Vec3f v3,
617 final Vec3f temp1V3, final Vec3f temp2V3, final Vec3f temp3V3) {
618 /**
619 Ax + By + Cz + D == 0 ;
620 D = - ( Ax + By + Cz )
621 = - ( A*a[0] + B*a[1] + C*a[2] )
622 = - vec3Dot ( normal, a ) ;
623 */
624 temp3V3.cross(temp1V3.minus(v2, v1), temp2V3.minus(v3, v1)).normalize();
625 resultVec4.set(temp3V3, -temp3V3.dot(v1));
626 return resultVec4;
627 }
628
629 /**
630 * Return intersection of an infinite line with a plane if exists, otherwise null.
631 * <p>
632 * Thanks to <i>Norman Vine -- nhv@yahoo.com (with hacks by Steve)</i>
633 * </p>
634 *
635 * @param result vec3 result buffer for intersecting coords
636 * @param ray here representing an infinite line, origin and direction.
637 * @param plane vec4 plane equation
638 * @param epsilon
639 * @return resulting intersecting if exists, otherwise null
640 */
641 public static Vec3f line2PlaneIntersection(final Vec3f result, final Ray ray, final Vec4f plane, final float epsilon) {
642 final Vec3f plane3 = new Vec3f(plane);
643 final float tmp = ray.dir.dot(plane3);
644
645 if ( FloatUtil.isZero(tmp, epsilon) ) {
646 return null; // ray is parallel to plane
647 }
648 result.set( ray.dir );
649 return result.scale( -( ray.orig.dot(plane3) + plane.w() ) / tmp ).add(ray.orig);
650 }
651
652 /**
653 * Compute intersection between two lines
654 * @param a vertex 1 of first line
655 * @param b vertex 2 of first line
656 * @param c vertex 1 of second line
657 * @param d vertex 2 of second line
658 * @return the intersection coordinates if the lines intersect, otherwise
659 * returns null
660 * @deprecated analyze correctness
661 */
662 @Deprecated
663 public static Vec3f line2lineIntersection0(final Vec3f result,
664 final Vert2fImmutable a, final Vert2fImmutable b,
665 final Vert2fImmutable c, final Vert2fImmutable d) {
666 final float determinant = (a.x()-b.x())*(c.y()-d.y()) - (a.y()-b.y())*(c.x()-d.x());
667
668 if (determinant == 0)
669 return null;
670
671 final float alpha = (a.x()*b.y()-a.y()*b.x());
672 final float beta = (c.x()*d.y()-c.y()*d.y());
673 final float xi = ((c.x()-d.x())*alpha-(a.x()-b.x())*beta)/determinant;
674 final float yi = ((c.y()-d.y())*alpha-(a.y()-b.y())*beta)/determinant;
675
676 return result.set(xi, yi, 0);
677 }
678 /** Compute intersection between two segments
679 * @param a vertex 1 of first segment
680 * @param b vertex 2 of first segment
681 * @param c vertex 1 of second segment
682 * @param d vertex 2 of second segment
683 * @return the intersection coordinates if the segments intersect, otherwise returns null
684 * @deprecated use {@link #seg2SegIntersection(Vec3f, Vec2f, Vec2f, Vec2f, Vec2f, boolean)
685 */
686 @Deprecated
687 public static Vec3f seg2SegIntersection0(final Vec3f result, final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c, final Vert2fImmutable d) {
688 final float determinant = (a.x()-b.x())*(c.y()-d.y()) - (a.y()-b.y())*(c.x()-d.x());
689
690 if (determinant == 0)
691 return null;
692
693 final float alpha = (a.x()*b.y()-a.y()*b.x());
694 final float beta = (c.x()*d.y()-c.y()*d.y());
695 final float xi = ((c.x()-d.x())*alpha-(a.x()-b.x())*beta)/determinant;
696 final float yi = ((c.y()-d.y())*alpha-(a.y()-b.y())*beta)/determinant;
697
698 final float gamma = (xi - a.x())/(b.x() - a.x());
699 final float gamma1 = (xi - c.x())/(d.x() - c.x());
700 if(gamma <= 0 || gamma >= 1) return null;
701 if(gamma1 <= 0 || gamma1 >= 1) return null;
702
703 return result.set(xi, yi, 0);
704 }
705 /**
706 * Compute intersection between two segments
707 * @param a vertex 1 of first segment
708 * @param b vertex 2 of first segment
709 * @param c vertex 1 of second segment
710 * @param d vertex 2 of second segment
711 * @return true if the segments intersect, otherwise returns false
712 * @deprecated use {@link #testSeg2SegIntersection(Vert2fImmutable, Vert2fImmutable, Vert2fImmutable, Vert2fImmutable, float, boolean)}
713 */
714 @Deprecated
715 public static boolean testSeg2SegIntersection0(final Vert2fImmutable a, final Vert2fImmutable b,
716 final Vert2fImmutable c, final Vert2fImmutable d) {
717 // Operations: 14+, 11*, 2 branches
718 final float determinant = (a.x()-b.x())*(c.y()-d.y()) - (a.y()-b.y())*(c.x()-d.x());
719
720 if (determinant == 0) {
721 return false;
722 }
723
724 final float alpha = (a.x()*b.y()-a.y()*b.x());
725 final float beta = (c.x()*d.y()-c.y()*d.y());
726 final float xi = ((c.x()-d.x())*alpha-(a.x()-b.x())*beta)/determinant;
727
728 final float gamma0 = (xi - a.x())/(b.x() - a.x());
729 final float gamma1 = (xi - c.x())/(d.x() - c.x());
730 if(gamma0 <= 0 || gamma0 >= 1 || gamma1 <= 0 || gamma1 >= 1) {
731 return false;
732 }
733
734 return true;
735 }
736 /**
737 * Check if a segment intersects with a triangle
738 * @param a vertex 1 of the triangle
739 * @param b vertex 2 of the triangle
740 * @param c vertex 3 of the triangle
741 * @param d vertex 1 of first segment
742 * @param e vertex 2 of first segment
743 * @return true if the segment intersects at least one segment of the triangle, false otherwise
744 * @deprecated use {@link #testTri2SegIntersection(Vert2fImmutable, Vert2fImmutable, Vert2fImmutable, Vert2fImmutable, Vert2fImmutable, float, boolean)}
745 */
746 @Deprecated
747 public static boolean testTri2SegIntersection0(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c,
748 final Vert2fImmutable d, final Vert2fImmutable e){
749 return testSeg2SegIntersection0(a, b, d, e) ||
750 testSeg2SegIntersection0(b, c, d, e) ||
751 testSeg2SegIntersection0(a, c, d, e) ;
752 }
753
754 /**
755 * Line segment intersection test and returning the intersecting point.
756 * <p>
757 * See [p + t r = q + u s](https://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect/565282#565282)
758 * and [its terse C# implementation](https://www.codeproject.com/tips/862988/find-the-intersection-point-of-two-line-segments)
759 * </p>
760 * <p>
761 * Implementation uses float precision.
762 * </p>
763 * @param p vertex 1 of first segment
764 * @param p2 vertex 2 of first segment
765 * @param q vertex 1 of second segment
766 * @param q2 vertex 2 of second segment
767 * @param epsilon optional epsilon. If {@code 0} just compare against {@code 0}, otherwise compare against this {@code epsilon}. See {@link FloatUtil#EPSILON}
768 * @param doCollinear consider collinear case, i.e. overlapping segments as an intersection returning {@code true}
769 * @return true if line segments are intersecting, otherwise false
770 */
771 public static boolean seg2SegIntersection(final Vec3f result, final Vec2f p, final Vec2f p2, final Vec2f q, final Vec2f q2, final float epsilon, final boolean doCollinear)
772 {
773 final Vec2f r = p2.minus(p);
774 final Vec2f s = q2.minus(q);
775 final float rxs = r.cross(s);
776
777 if ( FloatUtil.isZero(rxs, epsilon) ) {
778 if( doCollinear ) {
779 final Vec2f q_p = q.minus(p);
780 final float qpxr = q_p.cross(r);
781 if ( FloatUtil.isZero(qpxr, epsilon) ) // disabled collinear case
782 {
783 // 1) r x s = 0 and (q - p) x r = 0, the two lines are collinear.
784 final Vec2f p_q = p.minus(q);
785 final float qp_dot_r = q_p.dot(r);
786 final float pq_dot_s = p_q.dot(s);
787 // if ( ( 0 <= qp_dot_r && qp_dot_r <= r.dot(r) ) ||
788 // ( 0 <= pq_dot_s && pq_dot_s <= s.dot(s) ) )
789 if ( ( epsilon <= qp_dot_r && qp_dot_r - r.dot(r) <= epsilon ) ||
790 ( epsilon <= pq_dot_s && pq_dot_s - s.dot(s) <= epsilon ) )
791 {
792 // 1.1) 0 <= (q - p) · r <= r · r or 0 <= (p - q) · s <= s · s, the two lines are overlapping
793 // FIXME: result set to q2 endpoint, OK?
794 result.set(q2, 0);
795 return true;
796 }
797 // 1.2 the two lines are collinear but disjoint.
798 return false;
799 } else {
800 // 2) r × s = 0 and (q − p) × r ≠ 0, the two lines are parallel and non-intersecting.
801 return false;
802 }
803 } else {
804 // Not considering collinear case as an intersection
805 return false;
806 }
807 } else {
808 // r x s != 0
809 final Vec2f q_p = q.minus(p);
810 final float qpxr = q_p.cross(r);
811
812 // p + t r = q + u s
813 // (p + t r) × s = (q + u s) × s
814 // t (r × s) = (q − p) × s, with s x s = 0
815 // t = (q - p) x s / (r x s)
816 final float t = q_p.cross(s) / rxs;
817
818 // u = (p − q) × r / (s × r) = (q - p) x r / (r x s), with s × r = − r × s
819 final float u = qpxr / rxs;
820
821 // if ( (0 <= t && t <= 1) && (0 <= u && u <= 1) )
822 if ( (epsilon <= t && t - 1 <= epsilon) && (epsilon <= u && u - 1 <= epsilon) )
823 {
824 // 3) r × s ≠ 0 and 0 ≤ t ≤ 1 and 0 ≤ u ≤ 1, the two line segments meet at the point p + t * r = q + u * s.
825 result.set( p.plus( r.mul( t ) ), 0 ); // == q + (u * s)
826 return true;
827 }
828 }
829
830 return false;
831 }
832 /**
833 * Line segment intersection test.
834 * <p>
835 * See [p + t r = q + u s](https://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect/565282#565282)
836 * and [its terse C# implementation](https://www.codeproject.com/tips/862988/find-the-intersection-point-of-two-line-segments)
837 * </p>
838 * <p>
839 * Implementation uses float precision.
840 * </p>
841 * @param p vertex 1 of first segment
842 * @param p2 vertex 2 of first segment
843 * @param q vertex 1 of second segment
844 * @param q2 vertex 2 of second segment
845 * @param epsilon optional epsilon. If {@code 0} just compare against {@code 0}, otherwise compare against this {@code epsilon}. See {@link FloatUtil#EPSILON}
846 * @param doCollinear consider collinear case, i.e. overlapping segments as an intersection returning {@code true}
847 * @return true if line segments are intersecting, otherwise false
848 */
849 public static boolean testSeg2SegIntersection(final Vec2f p, final Vec2f p2, final Vec2f q, final Vec2f q2, final float epsilon, final boolean doCollinear)
850 {
851 final Vec2f r = p2.minus(p);
852 final Vec2f s = q2.minus(q);
853 final float rxs = r.cross(s);
854
855 if ( FloatUtil.isZero(rxs, epsilon) ) {
856 if( doCollinear ) {
857 final Vec2f q_p = q.minus(p);
858 final float qpxr = q_p.cross(r);
859 if ( FloatUtil.isZero(qpxr, epsilon) ) // disabled collinear case
860 {
861 // 1) r x s = 0 and (q - p) x r = 0, the two lines are collinear.
862 final Vec2f p_q = p.minus(q);
863 final float qp_dot_r = q_p.dot(r);
864 final float pq_dot_s = p_q.dot(s);
865 // if ( ( 0 <= qp_dot_r && qp_dot_r <= r.dot(r) ) ||
866 // ( 0 <= pq_dot_s && pq_dot_s <= s.dot(s) ) )
867 if ( ( epsilon <= qp_dot_r && qp_dot_r - r.dot(r) <= epsilon ) ||
868 ( epsilon <= pq_dot_s && pq_dot_s - s.dot(s) <= epsilon ) )
869 {
870 // 1.1) 0 <= (q - p) · r <= r · r or 0 <= (p - q) · s <= s · s, the two lines are overlapping
871 return true;
872 }
873 // 1.2 the two lines are collinear but disjoint.
874 return false;
875 } else {
876 // 2) r × s = 0 and (q − p) × r ≠ 0, the two lines are parallel and non-intersecting.
877 return false;
878 }
879 } else {
880 // Not considering collinear case as an intersection
881 return false;
882 }
883 } else {
884 // r x s != 0
885 final Vec2f q_p = q.minus(p);
886 final float qpxr = q_p.cross(r);
887
888 // p + t r = q + u s
889 // (p + t r) × s = (q + u s) × s
890 // t (r × s) = (q − p) × s, with s x s = 0
891 // t = (q - p) x s / (r x s)
892 final float t = q_p.cross(s) / rxs;
893
894 // u = (p − q) × r / (s × r) = (q - p) x r / (r x s), with s × r = − r × s
895 final float u = qpxr / rxs;
896
897 // if ( (0 <= t && t <= 1) && (0 <= u && u <= 1) )
898 if ( (epsilon <= t && t - 1 <= epsilon) && (epsilon <= u && u - 1 <= epsilon) )
899 {
900 // 3) r × s ≠ 0 and 0 ≤ t ≤ 1 and 0 ≤ u ≤ 1, the two line segments meet at the point p + t * r = q + u * s.
901 return true;
902 }
903 }
904 return false;
905 }
906 /**
907 * Line segment intersection test
908 * <p>
909 * See [p + t r = q + u s](https://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect/565282#565282)
910 * and [its terse C# implementation](https://www.codeproject.com/tips/862988/find-the-intersection-point-of-two-line-segments)
911 * </p>
912 * <p>
913 * Implementation uses float precision.
914 * </p>
915 * @param p vertex 1 of first segment
916 * @param p2 vertex 2 of first segment
917 * @param q vertex 1 of second segment
918 * @param q2 vertex 2 of second segment
919 * @param epsilon optional epsilon. If {@code 0} just compare against {@code 0}, otherwise compare against this {@code epsilon}. See {@link FloatUtil#EPSILON}
920 * @param doCollinear consider collinear case, i.e. overlapping segments as an intersection returning {@code true}
921 * @return true if line segments are intersecting, otherwise false
922 */
923 public static boolean testSeg2SegIntersection(final Vert2fImmutable p, final Vert2fImmutable p2, final Vert2fImmutable q, final Vert2fImmutable q2,
924 final float epsilon, final boolean doCollinear)
925 {
926 // Operations: 11+, 8*, 2 branches
927 final float rx = p2.x() - p.x(); // p2.minus(p)
928 final float ry = p2.y() - p.y();
929 final float sx = q2.x() - q.x(); // q2.minus(q)
930 final float sy = q2.y() - q.y();
931 final float rxs = rx * sy - ry * sx; // r.cross(s)
932
933 if ( FloatUtil.isZero(rxs, epsilon) ) {
934 if( doCollinear ) {
935 final float q_px = q.x() - p.x(); // q.minus(p)
936 final float q_py = q.y() - p.y();
937 final float qpxr = q_px * ry - q_py * rx; // q_p.cross(r)
938 if ( FloatUtil.isZero(qpxr, epsilon) ) // disabled collinear case
939 {
940 // 1) r x s = 0 and (q - p) x r = 0, the two lines are collinear.
941 final float p_qx = p.x() - q.x(); // p.minus(q)
942 final float p_qy = p.y() - q.y();
943 final float qp_dot_r = q_px * rx + q_py * ry; // q_p.dot(r);
944 final float pq_dot_s = p_qx * sx + p_qy * sy; // p_q.dot(s);
945 final float r_dot_r = rx * rx + ry * ry; // r.dot(r);
946 final float s_dot_s = sx * sx + sy * sy; // r.dot(r);
947 // if ( ( 0 <= qp_dot_r && qp_dot_r <= r.dot(r) ) ||
948 // ( 0 <= pq_dot_s && pq_dot_s <= s.dot(s) ) )
949 if ( ( epsilon <= qp_dot_r && qp_dot_r - r_dot_r <= epsilon ) ||
950 ( epsilon <= pq_dot_s && pq_dot_s - s_dot_s <= epsilon ) )
951 {
952 // 1.1) 0 <= (q - p) · r <= r · r or 0 <= (p - q) · s <= s · s, the two lines are overlapping
953 return true;
954 }
955 // 1.2 the two lines are collinear but disjoint.
956 return false;
957 } else {
958 // 2) r × s = 0 and (q − p) × r ≠ 0, the two lines are parallel and non-intersecting.
959 return false;
960 }
961 } else {
962 // Not considering collinear case as an intersection
963 return false;
964 }
965 } else {
966 // r x s != 0
967 final float q_px = q.x() - p.x(); // q.minus(p)
968 final float q_py = q.y() - p.y();
969 final float qpxr = q_px * ry - q_py * rx; // q_p.cross(r)
970
971 // p + t r = q + u s
972 // (p + t r) × s = (q + u s) × s
973 // t (r × s) = (q − p) × s, with s x s = 0
974 // t = (q - p) x s / (r x s)
975 final float t = ( q_px * sy - q_py * sx ) / rxs; // q_p.cross(s) / rxs
976
977 // u = (p − q) × r / (s × r) = (q - p) x r / (r x s), with s × r = − r × s
978 final float u = qpxr / rxs;
979
980 // if ( (0 <= t && t <= 1) && (0 <= u && u <= 1) )
981 if ( (epsilon <= t && t - 1 <= epsilon) && (epsilon <= u && u - 1 <= epsilon) )
982 {
983 // 3) r × s ≠ 0 and 0 ≤ t ≤ 1 and 0 ≤ u ≤ 1, the two line segments meet at the point p + t * r = q + u * s.
984 return true;
985 }
986 }
987 return false;
988 }
989 /**
990 * Check if a segment intersects with a triangle.
991 * <p>
992 * Implementation uses float precision.
993 * </p>
994 * @param a vertex 1 of the triangle
995 * @param b vertex 2 of the triangle
996 * @param c vertex 3 of the triangle
997 * @param d vertex 1 of first segment
998 * @param e vertex 2 of first segment
999 * @param epsilon optional epsilon. If {@code 0} just compare against {@code 0}, otherwise compare against this {@code epsilon}. See {@link FloatUtil#EPSILON}
1000 * @param doCollinear consider collinear case, i.e. overlapping segments as an intersection returning {@code true}
1001 * @return true if the segment intersects at least one segment of the triangle, false otherwise
1002 * @see #testSeg2SegIntersection(Vert2fImmutable, Vert2fImmutable, Vert2fImmutable, Vert2fImmutable, float, boolean)
1003 */
1004 public static boolean testTri2SegIntersection(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c,
1005 final Vert2fImmutable d, final Vert2fImmutable e, final float epsilon, final boolean doCollinear) {
1006 return testSeg2SegIntersection(a, b, d, e, epsilon, doCollinear) ||
1007 testSeg2SegIntersection(b, c, d, e, epsilon, doCollinear) ||
1008 testSeg2SegIntersection(a, c, d, e, epsilon, doCollinear) ;
1009 }
1010
1011 /**
1012 * Line segment intersection test using {@link FloatUtil#EPSILON} w/o considering collinear-case
1013 * <p>
1014 * See [p + t r = q + u s](https://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect/565282#565282)
1015 * and [its terse C# implementation](https://www.codeproject.com/tips/862988/find-the-intersection-point-of-two-line-segments)
1016 * </p>
1017 * <p>
1018 * Implementation uses float precision.
1019 * </p>
1020 * @param p vertex 1 of first segment
1021 * @param p2 vertex 2 of first segment
1022 * @param q vertex 1 of second segment
1023 * @param q2 vertex 2 of second segment
1024 * @return true if line segments are intersecting, otherwise false
1025 */
1026 public static boolean testSeg2SegIntersection(final Vert2fImmutable p, final Vert2fImmutable p2, final Vert2fImmutable q, final Vert2fImmutable q2)
1027 {
1028 // Operations: 11+, 8*, 2 branches
1029 final float rx = p2.x() - p.x(); // p2.minus(p)
1030 final float ry = p2.y() - p.y();
1031 final float sx = q2.x() - q.x(); // q2.minus(q)
1032 final float sy = q2.y() - q.y();
1033 final float rxs = rx * sy - ry * sx; // r.cross(s)
1034
1035 final float epsilon = FloatUtil.EPSILON;
1036
1037 if ( FloatUtil.isZero(rxs) ) {
1038 // Not considering collinear case as an intersection
1039 return false;
1040 } else {
1041 // r x s != 0
1042 final float q_px = q.x() - p.x(); // q.minus(p)
1043 final float q_py = q.y() - p.y();
1044 final float qpxr = q_px * ry - q_py * rx; // q_p.cross(r)
1045
1046 // p + t r = q + u s
1047 // (p + t r) × s = (q + u s) × s
1048 // t (r × s) = (q − p) × s, with s x s = 0
1049 // t = (q - p) x s / (r x s)
1050 final float t = ( q_px * sy - q_py * sx ) / rxs; // q_p.cross(s) / rxs
1051
1052 // u = (p − q) × r / (s × r) = (q - p) x r / (r x s), with s × r = − r × s
1053 final float u = qpxr / rxs;
1054
1055 // if ( (0 <= t && t <= 1) && (0 <= u && u <= 1) )
1056 if ( (epsilon <= t && t - 1 <= epsilon) && (epsilon <= u && u - 1 <= epsilon) )
1057 {
1058 // 3) r × s ≠ 0 and 0 ≤ t ≤ 1 and 0 ≤ u ≤ 1, the two line segments meet at the point p + t * r = q + u * s.
1059 return true;
1060 }
1061 }
1062 return false;
1063 }
1064 /**
1065 * Check if a segment intersects with a triangle using {@link FloatUtil#EPSILON} w/o considering collinear-case
1066 * <p>
1067 * Implementation uses float precision.
1068 * </p>
1069 * @param a vertex 1 of the triangle
1070 * @param b vertex 2 of the triangle
1071 * @param c vertex 3 of the triangle
1072 * @param d vertex 1 of first segment
1073 * @param e vertex 2 of first segment
1074 * @return true if the segment intersects at least one segment of the triangle, false otherwise
1075 * @see #testSeg2SegIntersection(Vert2fImmutable, Vert2fImmutable, Vert2fImmutable, Vert2fImmutable, float, boolean)
1076 */
1077 public static boolean testTri2SegIntersection(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c,
1078 final Vert2fImmutable d, final Vert2fImmutable e) {
1079 return testSeg2SegIntersection(a, b, d, e) ||
1080 testSeg2SegIntersection(b, c, d, e) ||
1081 testSeg2SegIntersection(a, c, d, e) ;
1082 }
1083
1084 /**
1085 * Returns whether the given {@code polyline} denotes a convex shape with O(n) complexity.
1086 * <p>
1087 * See [Determine whether a polygon is convex based on its vertices](https://math.stackexchange.com/questions/1743995/determine-whether-a-polygon-is-convex-based-on-its-vertices/1745427#1745427)
1088 * </p>
1089 * @param polyline connected {@link Vert2fImmutable}, i.e. a poly-line
1090 * @param shortIsConvex return value if {@code vertices} have less than three elements, allows simplification
1091 */
1092 public static boolean isConvex0(final List<? extends Vert2fImmutable> polyline, final boolean shortIsConvex) {
1093 final int polysz = polyline.size();
1094 if( polysz < 3 ) {
1095 return shortIsConvex;
1096 }
1097 final float eps = FloatUtil.EPSILON;
1098
1099 float wSign = 0; // First nonzero orientation (positive or negative)
1100
1101 int xSign = 0;
1102 int xFirstSign = 0; // Sign of first nonzero edge vector x
1103 int xFlips = 0; // Number of sign changes in x
1104
1105 int ySign = 0;
1106 int yFirstSign = 0; // Sign of first nonzero edge vector y
1107 int yFlips = 0; // Number of sign changes in y
1108
1109 Vert2fImmutable curr = polyline.get(polysz-2); // Second-to-last vertex
1110 Vert2fImmutable next = polyline.get(polysz-1); // Last vertex
1111
1112 for(int i=0; i<polysz; ++i) { // Each vertex, in order
1113 final Vert2fImmutable prev = curr; // Previous vertex
1114 curr = next; // Current vertex
1115 next = polyline.get(i); // Next vertex
1116
1117 // Previous edge vector ("before"):
1118 final float bx = curr.x() - prev.x();
1119 final float by = curr.y() - prev.y();
1120
1121 // Next edge vector ("after"):
1122 final float ax = next.x() - curr.x();
1123 final float ay = next.y() - curr.y();
1124
1125 // Calculate sign flips using the next edge vector ("after"),
1126 // recording the first sign.
1127 if( ax > eps ) {
1128 if( xSign == 0 ) {
1129 xFirstSign = +1;
1130 } else if( xSign < 0 ) {
1131 xFlips = xFlips + 1;
1132 }
1133 xSign = +1;
1134 } else if( ax < -eps ) {
1135 if( xSign == 0 ) {
1136 xFirstSign = -1;
1137 } else if ( xSign > 0 ) {
1138 xFlips = xFlips + 1;
1139 }
1140 xSign = -1;
1141 }
1142 if( xFlips > 2 ) {
1143 return false;
1144 }
1145
1146 if( ay > eps ) {
1147 if( ySign == 0 ) {
1148 yFirstSign = +1;
1149 } else if( ySign < 0 ) {
1150 yFlips = yFlips + 1;
1151 }
1152 ySign = +1;
1153 } else if( ay < -eps ) {
1154 if( ySign == 0 ) {
1155 yFirstSign = -1;
1156 } else if( ySign > 0 ) {
1157 yFlips = yFlips + 1;
1158 }
1159 ySign = -1;
1160 }
1161 if( yFlips > 2 ) {
1162 return false;
1163 }
1164
1165 // Find out the orientation of this pair of edges,
1166 // and ensure it does not differ from previous ones.
1167 final float w = bx*ay - ax*by;
1168 if( FloatUtil.isZero(wSign) && !FloatUtil.isZero(w) ) {
1169 wSign = w;
1170 } else if( wSign > eps && w < -eps ) {
1171 return false;
1172 } else if( wSign < -eps && w > eps ) {
1173 return false;
1174 }
1175 }
1176
1177 // Final/wraparound sign flips:
1178 if( xSign != 0 && xFirstSign != 0 && xSign != xFirstSign ) {
1179 xFlips = xFlips + 1;
1180 }
1181 if( ySign != 0 && yFirstSign != 0 && ySign != yFirstSign ) {
1182 yFlips = yFlips + 1;
1183 }
1184
1185 // Concave polygons have two sign flips along each axis.
1186 if( xFlips != 2 || yFlips != 2 ) {
1187 return false;
1188 }
1189
1190 // This is a convex polygon.
1191 return true;
1192 }
1193
1194
1195 private static int cmod(final int i, final int count) {
1196 if( i >= 0 ) {
1197 return i % count;
1198 } else {
1199 return i + count;
1200 }
1201 }
1202 /**
1203 * Returns whether the given on-curve {@code polyline} points denotes a convex shape with O(n) complexity.
1204 * <p>
1205 * See [Determine whether a polygon is convex based on its vertices](https://math.stackexchange.com/questions/1743995/determine-whether-a-polygon-is-convex-based-on-its-vertices/1745427#1745427)
1206 * </p>
1207 * <p>
1208 * All off-curve points are ignored.
1209 * </p>
1210 * @param polyline connected {@link Vert2fImmutable}, i.e. a poly-line
1211 * @param shortIsConvex return value if {@code vertices} have less than three elements, allows simplification
1212 */
1213 public static boolean isConvex1(final List<Vertex> polyline, final boolean shortIsConvex) {
1214 final int polysz = polyline.size();
1215 if( polysz < 3 ) {
1216 return shortIsConvex;
1217 }
1218 final float eps = FloatUtil.EPSILON;
1219
1220 float wSign = 0; // First nonzero orientation (positive or negative)
1221
1222 int xSign = 0;
1223 int xFirstSign = 0; // Sign of first nonzero edge vector x
1224 int xFlips = 0; // Number of sign changes in x
1225
1226 int ySign = 0;
1227 int yFirstSign = 0; // Sign of first nonzero edge vector y
1228 int yFlips = 0; // Number of sign changes in y
1229
1230 int offset=-3;
1231 Vertex v0, v1;
1232 {
1233 do {
1234 ++offset; // -2
1235 v0 = polyline.get(cmod(offset, polysz)); // current, polyline[-2] if on-curve
1236 } while( !v0.isOnCurve() && offset < polysz );
1237 if( offset >= polysz ) {
1238 return shortIsConvex;
1239 }
1240 do {
1241 ++offset; // -1
1242 v1 = polyline.get(cmod(offset, polysz)); // next, polyline[-1] if both on-curve
1243 } while( !v1.isOnCurve() && offset < polysz );
1244 if( offset >= polysz ) {
1245 return shortIsConvex;
1246 }
1247 }
1248
1249 while( offset < polysz ) {
1250 final Vertex vp = v0; // previous on-curve vertex
1251 v0 = v1; // current on-curve vertex
1252 do {
1253 ++offset; // 0, ...
1254 v1 = polyline.get(cmod(offset, polysz)); // next on-curve vertex
1255 } while( !v1.isOnCurve() && offset < polysz );
1256 if( offset >= polysz ) {
1257 break;
1258 }
1259
1260 // Previous edge vector ("before"):
1261 final float bx = v0.x() - vp.x();
1262 final float by = v0.y() - vp.y();
1263
1264 // Next edge vector ("after"):
1265 final float ax = v1.x() - v0.x();
1266 final float ay = v1.y() - v0.y();
1267
1268 // Calculate sign flips using the next edge vector ("after"),
1269 // recording the first sign.
1270 if( ax > eps ) {
1271 if( xSign == 0 ) {
1272 xFirstSign = +1;
1273 } else if( xSign < 0 ) {
1274 xFlips = xFlips + 1;
1275 }
1276 xSign = +1;
1277 } else if( ax < -eps ) {
1278 if( xSign == 0 ) {
1279 xFirstSign = -1;
1280 } else if ( xSign > 0 ) {
1281 xFlips = xFlips + 1;
1282 }
1283 xSign = -1;
1284 }
1285 if( xFlips > 2 ) {
1286 return false;
1287 }
1288
1289 if( ay > eps ) {
1290 if( ySign == 0 ) {
1291 yFirstSign = +1;
1292 } else if( ySign < 0 ) {
1293 yFlips = yFlips + 1;
1294 }
1295 ySign = +1;
1296 } else if( ay < -eps ) {
1297 if( ySign == 0 ) {
1298 yFirstSign = -1;
1299 } else if( ySign > 0 ) {
1300 yFlips = yFlips + 1;
1301 }
1302 ySign = -1;
1303 }
1304 if( yFlips > 2 ) {
1305 return false;
1306 }
1307
1308 // Find out the orientation of this pair of edges,
1309 // and ensure it does not differ from previous ones.
1310 final float w = bx*ay - ax*by;
1311 if( FloatUtil.isZero(wSign) && !FloatUtil.isZero(w) ) {
1312 wSign = w;
1313 } else if( wSign > eps && w < -eps ) {
1314 return false;
1315 } else if( wSign < -eps && w > eps ) {
1316 return false;
1317 }
1318 }
1319
1320 // Final/wraparound sign flips:
1321 if( xSign != 0 && xFirstSign != 0 && xSign != xFirstSign ) {
1322 xFlips = xFlips + 1;
1323 }
1324 if( ySign != 0 && yFirstSign != 0 && ySign != yFirstSign ) {
1325 yFlips = yFlips + 1;
1326 }
1327
1328 // Concave polygons have two sign flips along each axis.
1329 if( xFlips != 2 || yFlips != 2 ) {
1330 return false;
1331 }
1332
1333 // This is a convex polygon.
1334 return true;
1335 }
1336
1337 /**
1338 * Returns whether the given {@code polyline} is self-intersecting shape (crossing over) with O(n*n) complexity.
1339 */
1340 public static boolean isSelfIntersecting1(final List<Vertex> polyline) {
1341 final int polysz = polyline.size();
1342 if (polysz < 4) {
1343 return false;
1344 }
1345 for (int i = 0; i < polysz-1; i++) {
1346 final Vertex p0 = polyline.get(i);
1347 final Vertex p1 = polyline.get(i+1);
1348 for (int j = i+2; j < polysz; j++) {
1349 if( i != 0 || j != (polysz-1) ) {
1350 final Vertex q0 = polyline.get(j);
1351 final Vertex q1 = polyline.get((j+1)%polysz);
1352
1353 if( VectorUtil.testSeg2SegIntersection(p0, p1, q0, q1, 0, false) ) {
1354 return true;
1355 }
1356 }
1357 }
1358 }
1359 return false;
1360 }
1361}
A Vertex exposing Vec3f vertex- and texture-coordinates.
Definition: Vertex.java:37
final boolean isOnCurve()
Definition: Vertex.java:146
Basic Double math utility functions.
Definition: DoubleUtil.java:33
static final double EPSILON
Epsilon for floating point {@value}, as once computed via getMachineEpsilon() on an AMD-64 CPU.
Basic Float math utility functions.
Definition: FloatUtil.java:83
static int compare(final float a, final float b)
Returns -1, 0 or 1 if a is less, equal or greater than b, disregarding epsilon but considering NaN,...
static final float EPSILON
Epsilon for floating point {@value}, as once computed via getMachineEpsilon() on an AMD-64 CPU.
static boolean isZero(final float a, final float epsilon)
Returns true if value is zero, i.e.
static float sqrt(final float a)
Simple compound denoting a ray.
Definition: Ray.java:49
final Vec3f dir
Normalized direction vector of ray.
Definition: Ray.java:54
final Vec3f orig
Origin of Ray.
Definition: Ray.java:51
2D Vector based upon two float components.
Definition: Vec2f.java:37
float cross(final Vec2f o)
Returns cross product of this vectors and the given one, i.e.
Definition: Vec2f.java:335
Vec2f plus(final Vec2f arg)
Returns this + arg; creates new vector.
Definition: Vec2f.java:198
Vec2f minus(final Vec2f arg)
Returns this - arg; creates new vector.
Definition: Vec2f.java:224
float dot(final Vec2f arg)
Return the dot product of this vector and the given one.
Definition: Vec2f.java:324
Vec2f mul(final float val)
Returns this * val; creates new vector.
Definition: Vec2f.java:155
3D Vector based upon three float components.
Definition: Vec3f.java:37
Vec3f scale(final float s)
this = this * s, returns this.
Definition: Vec3f.java:218
Vec3f normalize()
Normalize this vector in place.
Definition: Vec3f.java:297
float dot(final Vec3f o)
Return the dot product of this vector and the given one.
Definition: Vec3f.java:338
Vec3f minus(final Vec3f arg)
Returns this - arg; creates new vector.
Definition: Vec3f.java:255
Vec3f cross(final Vec3f arg)
Returns this cross arg; creates new vector.
Definition: Vec3f.java:343
Vec3f set(final Vec3f o)
this = o, returns this.
Definition: Vec3f.java:79
Vec3f add(final float dx, final float dy, final float dz)
this = this + { dx, dy, dz }, returns this.
Definition: Vec3f.java:239
4D Vector based upon four float components.
Definition: Vec4f.java:37
Vec4f set(final Vec4f o)
this = o, returns this.
Definition: Vec4f.java:67
static float[] crossVec3(final float[] r, final int r_offset, final float[] v1, final int v1_offset, final float[] v2, final int v2_offset)
cross product vec1 x vec2
static float determinant(final Vec3f a, final Vec3f b, final Vec3f c)
Return the determinant of 3 vectors.
static Vec3f midpoint(final Vec3f result, final Vec3f p1, final Vec3f p2)
Calculate the midpoint of two points.
static boolean isInTriangle3(final Vec3f a, final Vec3f b, final Vec3f c, final Vec3f p1, final Vec3f p2, final Vec3f p3, final Vec3f ac, final Vec3f ab, final Vec3f ap)
Check if one of three vertices are in triangle using barycentric coordinates computation.
static boolean testSeg2SegIntersection0(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c, final Vert2fImmutable d)
Compute intersection between two segments.
static float[] subVec2(final float[] result, final float[] v1, final float[] v2)
Subtracts two vectors, result = v1 - v2.
static boolean testSeg2SegIntersection(final Vert2fImmutable p, final Vert2fImmutable p2, final Vert2fImmutable q, final Vert2fImmutable q2, final float epsilon, final boolean doCollinear)
Line segment intersection test.
static boolean testSeg2SegIntersection(final Vec2f p, final Vec2f p2, final Vec2f q, final Vec2f q2, final float epsilon, final boolean doCollinear)
Line segment intersection test.
static Winding getWinding(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c)
Compute the winding of the 3 given points.
static float[] scaleVec2(final float[] result, final float[] vector, final float scale)
Scales a vector by param using given result float[], result = vector * scale.
static boolean seg2SegIntersection(final Vec3f result, final Vec2f p, final Vec2f p2, final Vec2f q, final Vec2f q2, final float epsilon, final boolean doCollinear)
Line segment intersection test and returning the intersecting point.
static boolean isCCW(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c)
Check if points are in ccw order.
static boolean isConvex0(final List<? extends Vert2fImmutable > polyline, final boolean shortIsConvex)
Returns whether the given polyline denotes a convex shape with O(n) complexity.
static double triArea(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c)
Computes oriented double area of a triangle, i.e.
static float normSquareVec2(final float[] vec)
Return the squared length of a vector, a.k.a the squared norm or squared magnitude
Definition: VectorUtil.java:91
static Winding getWinding(final List<? extends Vert2fImmutable > vertices)
Compute the winding using the area(List) function over all vertices for complex shapes.
static float[] divVec2(final float[] result, final float[] vector, final float scale)
Divides a vector by param using given result float[], result = vector / scale.
static Vec3f seg2SegIntersection0(final Vec3f result, final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c, final Vert2fImmutable d)
Compute intersection between two segments.
static boolean testTri2SegIntersection(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c, final Vert2fImmutable d, final Vert2fImmutable e)
Check if a segment intersects with a triangle using FloatUtil#EPSILON w/o considering collinear-case.
static float[] addVec2(final float[] result, final float[] v1, final float[] v2)
Adds two vectors, result = v1 + v2.
static float[] divVec2(final float[] result, final float[] vector, final float[] scale)
Divides a vector by param using given result float[], result = vector / scale.
static boolean testTri2SegIntersection0(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c, final Vert2fImmutable d, final Vert2fImmutable e)
Check if a segment intersects with a triangle.
static float[] normalizeVec3(final float[] vector)
Normalize a vector in place.
static boolean isInTriangle3(final Vec3f a, final Vec3f b, final Vec3f c, final Vec3f p1, final Vec3f p2, final Vec3f p3, final Vec3f ac, final Vec3f ab, final Vec3f ap, final float epsilon)
Check if one of three vertices are in triangle using barycentric coordinates computation,...
static boolean isSelfIntersecting1(final List< Vertex > polyline)
Returns whether the given polyline is self-intersecting shape (crossing over) with O(n*n) complexity.
static double inCircleVal(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c, final Vert2fImmutable d)
static boolean isInCircle(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c, final Vert2fImmutable d)
Check if vertices in triangle circumcircle given d vertex, from paper by Guibas and Stolfi (1985).
static float normSquareVec3(final float[] vec)
Return the squared length of a vector, a.k.a the squared norm or squared magnitude
Definition: VectorUtil.java:98
static boolean isInTriangle(final Vec3f a, final Vec3f b, final Vec3f c, final Vec3f p, final Vec3f ac, final Vec3f ab, final Vec3f ap)
Check if a vertex is in triangle using barycentric coordinates computation.
static boolean isZero(final float x, final float y, final float z)
Return true if all three vector components are zero, i.e.
Definition: VectorUtil.java:61
static boolean isZero(final float x, final float y, final float z, final float epsilon)
Return true if all three vector components are zero, i.e.
Definition: VectorUtil.java:49
static float normVec2(final float[] vec)
Return the length of a vector, a.k.a the norm or magnitude
static boolean testSeg2SegIntersection(final Vert2fImmutable p, final Vert2fImmutable p2, final Vert2fImmutable q, final Vert2fImmutable q2)
Line segment intersection test using FloatUtil#EPSILON w/o considering collinear-case.
static Vec3f line2PlaneIntersection(final Vec3f result, final Ray ray, final Vec4f plane, final float epsilon)
Return intersection of an infinite line with a plane if exists, otherwise null.
static boolean isVec2Zero(final Vec3f vec)
Return true if 2D vector components are zero, no FloatUtil#EPSILON is taken into consideration.
Definition: VectorUtil.java:39
static float distVec3(final float[] v1, final float[] v2)
Return the distance between the given two points described vector v1 and v2.
Definition: VectorUtil.java:84
static double area(final List<? extends Vert2fImmutable > vertices)
Computes the area of a list of vertices via shoelace formula.
static Vec3f line2lineIntersection0(final Vec3f result, final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c, final Vert2fImmutable d)
Compute intersection between two lines.
static boolean testTri2SegIntersection(final Vert2fImmutable a, final Vert2fImmutable b, final Vert2fImmutable c, final Vert2fImmutable d, final Vert2fImmutable e, final float epsilon, final boolean doCollinear)
Check if a segment intersects with a triangle.
static boolean isConvex1(final List< Vertex > polyline, final boolean shortIsConvex)
Returns whether the given on-curve polyline points denotes a convex shape with O(n) complexity.
static Vec4f getPlaneVec3(final Vec4f resultV4, final Vec3f normalVec3, final Vec3f pVec3)
Finds the plane equation of a plane given its normal and a point on the plane.
static float distSquareVec3(final float[] v1, final float[] v2)
Return the squared distance between the given two points described vector v1 and v2.
Definition: VectorUtil.java:74
static float[] scaleVec2(final float[] result, final float[] vector, final float[] scale)
Scales a vector by param using given result float[], result = vector * scale.
static float normSquareVec3(final float[] vec, final int offset)
Return the squared length of a vector, a.k.a the squared norm or squared magnitude
static float[] normalizeVec3(final float[] vector, final int offset)
Normalize a vector in place.
static Vec4f getPlaneVec3(final Vec4f resultVec4, final Vec3f v1, final Vec3f v2, final Vec3f v3, final Vec3f temp1V3, final Vec3f temp2V3, final Vec3f temp3V3)
This finds the plane equation of a triangle given three vertices.
static boolean isCollinear(final Vec3f v1, final Vec3f v2, final Vec3f v3)
Check if three vertices are collinear.
Winding direction, either clockwise (CW) or counter-clockwise (CCW).
Definition: Winding.java:6