1 /**
2  * Efficiently performs frustum intersection tests by caching the frustum planes of an arbitrary transformation {@link Matrix4d matrix}.
3  * <p>
4  * This class is preferred over the frustum intersection methods in {@link Matrix4d} when many objects need to be culled by the same static frustum.
5  * 
6  * @author Kai Burjack
7  */
8 module doml.frustum_intersection;
9 
10 import Math = doml.math;
11 
12 import doml.matrix_4d;
13 
14 import doml.vector_2d;
15 import doml.vector_3d;
16 import doml.vector_4d;
17 
18 /*
19  * The MIT License
20  *
21  * Copyright (c) 2015-2021 Kai Burjack
22  *
23  * Permission is hereby granted, free of charge, to any person obtaining a copy
24  * of this software and associated documentation files (the "Software"), to deal
25  * in the Software without restriction, including without limitation the rights
26  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
27  * copies of the Software, and to permit persons to whom the Software is
28  * furnished to do so, subject to the following conditions:
29  *
30  * The above copyright notice and this permission notice shall be included in
31  * all copies or substantial portions of the Software.
32  *
33  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
34  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
35  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
36  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
37  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
38  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
39  * THE SOFTWARE.
40  */
41 
42 /**
43  * Efficiently performs frustum intersection tests by caching the frustum planes of an arbitrary transformation {@link Matrix4d matrix}.
44  * <p>
45  * This class is preferred over the frustum intersection methods in {@link Matrix4d} when many objects need to be culled by the same static frustum.
46  * 
47  * @author Kai Burjack
48  */
49 struct FrustumIntersection {
50 
51     /**
52      * Return value of {@link #intersectAab(double, double, double, double, double, double) intersectAab()}
53      * and its different overloads identifying the plane with equation <code>x=-1</code> when using the identity frustum.
54      */
55     static immutable int PLANE_NX = 0;
56     /**
57      * Return value of {@link #intersectAab(double, double, double, double, double, double) intersectAab()}
58      * and its different overloads identifying the plane with equation <code>x=1</code> when using the identity frustum.
59      */
60     static immutable int PLANE_PX = 1;
61     /**
62      * Return value of {@link #intersectAab(double, double, double, double, double, double) intersectAab()}
63      * and its different overloads identifying the plane with equation <code>y=-1</code> when using the identity frustum.
64      */
65     static immutable int PLANE_NY= 2;
66     /**
67      * Return value of {@link #intersectAab(double, double, double, double, double, double) intersectAab()}
68      * and its different overloads identifying the plane with equation <code>y=1</code> when using the identity frustum.
69      */
70     static immutable int PLANE_PY = 3;
71     /**
72      * Return value of {@link #intersectAab(double, double, double, double, double, double) intersectAab()}
73      * and its different overloads identifying the plane with equation <code>z=-1</code> when using the identity frustum.
74      */
75     static immutable int PLANE_NZ = 4;
76     /**
77      * Return value of {@link #intersectAab(double, double, double, double, double, double) intersectAab()}
78      * and its different overloads identifying the plane with equation <code>z=1</code> when using the identity frustum.
79      */
80     static immutable int PLANE_PZ = 5;
81 
82     /**
83      * Return value of {@link #intersectAab(double, double, double, double, double, double) intersectAab()}
84      * and its different overloads indicating that the axis-aligned box intersects the frustum.
85      */
86     static immutable int INTERSECT = -1;
87     /**
88      * Return value of {@link #intersectAab(double, double, double, double, double, double) intersectAab()}
89      * and its different overloads indicating that the axis-aligned box is fully inside of the frustum.
90      */
91     static immutable int INSIDE = -2;
92     /**
93      * Return value of {@link #intersectSphere(Vector3d, double)} or {@link #intersectSphere(double, double, double, double)}
94      * indicating that the sphere is completely outside of the frustum.
95      */
96     static immutable int OUTSIDE = -3;
97 
98     /**
99      * The value in a bitmask for
100      * {@link #intersectAab(double, double, double, double, double, double, int) intersectAab()}
101      * that identifies the plane with equation <code>x=-1</code> when using the identity frustum.
102      */
103     static immutable int PLANE_MASK_NX = 1<<PLANE_NX;
104     /**
105      * The value in a bitmask for
106      * {@link #intersectAab(double, double, double, double, double, double, int) intersectAab()}
107      * that identifies the plane with equation <code>x=1</code> when using the identity frustum.
108      */
109     static immutable int PLANE_MASK_PX = 1<<PLANE_PX;
110     /**
111      * The value in a bitmask for
112      * {@link #intersectAab(double, double, double, double, double, double, int) intersectAab()}
113      * that identifies the plane with equation <code>y=-1</code> when using the identity frustum.
114      */
115     static immutable int PLANE_MASK_NY = 1<<PLANE_NY;
116     /**
117      * The value in a bitmask for
118      * {@link #intersectAab(double, double, double, double, double, double, int) intersectAab()}
119      * that identifies the plane with equation <code>y=1</code> when using the identity frustum.
120      */
121     static immutable int PLANE_MASK_PY = 1<<PLANE_PY;
122     /**
123      * The value in a bitmask for
124      * {@link #intersectAab(double, double, double, double, double, double, int) intersectAab()}
125      * that identifies the plane with equation <code>z=-1</code> when using the identity frustum.
126      */
127     static immutable int PLANE_MASK_NZ = 1<<PLANE_NZ;
128     /**
129      * The value in a bitmask for
130      * {@link #intersectAab(double, double, double, double, double, double, int) intersectAab()}
131      * that identifies the plane with equation <code>z=1</code> when using the identity frustum.
132      */
133     static immutable int PLANE_MASK_PZ = 1<<PLANE_PZ;
134 
135     private double nxX = 0.0;
136     private double nxY = 0.0; 
137     private double nxZ = 0.0;
138     private double nxW = 0.0;
139 
140     private double pxX = 0.0; 
141     private double pxY = 0.0; 
142     private double pxZ = 0.0; 
143     private double pxW = 0.0;
144 
145     private double nyX = 0.0; 
146     private double nyY = 0.0; 
147     private double nyZ = 0.0; 
148     private double nyW = 0.0;
149 
150     private double pyX = 0.0;
151     private double pyY = 0.0;
152     private double pyZ = 0.0;
153     private double pyW = 0.0;
154 
155     private double nzX = 0.0;
156     private double nzY = 0.0;
157     private double nzZ = 0.0;
158     private double nzW = 0.0;
159 
160     private double pzX = 0.0;
161     private double pzY = 0.0;
162     private double pzZ = 0.0;
163     private double pzW = 0.0;
164 
165     private Vector4d[] planes = {
166         Vector4d[6] temp;
167         for (int i = 0; i < 6; i++) {
168             temp[i] = Vector4d();
169         }
170         return temp;
171     } ();
172 
173     /**
174      * Create a new {@link FrustumIntersection} from the given {@link Matrix4d matrix} by extracing the matrix's frustum planes.
175      * <p>
176      * In order to update the compute frustum planes later on, call {@link #set(Matrix4d)}.
177      * 
178      * @see #set(Matrix4d)
179      * 
180      * @param m
181      *          the {@link Matrix4d} to create the frustum culler from
182      */
183     this(Matrix4d m) {
184         set(m, true);
185     }
186 
187     /**
188      * Create a new {@link FrustumIntersection} from the given {@link Matrix4d matrix} by extracing the matrix's frustum planes.
189      * <p>
190      * In order to update the compute frustum planes later on, call {@link #set(Matrix4d)}.
191      * 
192      * @see #set(Matrix4d)
193      * 
194      * @param m
195      *          the {@link Matrix4d} to create the frustum culler from
196      * @param allowTestSpheres
197      *          whether the methods {@link #testSphere(Vector3d, double)}, {@link #testSphere(double, double, double, double)},
198      *          {@link #intersectSphere(Vector3d, double)} or {@link #intersectSphere(double, double, double, double)} will used.
199      *          If no spheres need to be tested, then <code>false</code> should be used 
200      */
201     this(Matrix4d m, bool allowTestSpheres) {
202         set(m, allowTestSpheres);
203     }
204 
205     /**
206      * Update the stored frustum planes of <code>this</code> {@link FrustumIntersection} with the given {@link Matrix4d matrix}.
207      * <p>
208      * Reference: <a href="http://gamedevs.org/uploads/fast-extraction-viewing-frustum-planes-from-world-view-projection-matrix.pdf">
209      * Fast Extraction of Viewing Frustum Planes from the World-View-Projection Matrix</a>
210      * 
211      * @param m
212      *          the {@link Matrix4d matrix} to update <code>this</code> frustum culler's frustum planes from
213      * @return this
214      */
215     ref public FrustumIntersection set(Matrix4d m) return {
216         return set(m, true);
217     }
218 
219     /**
220      * Update the stored frustum planes of <code>this</code> {@link FrustumIntersection} with the given {@link Matrix4d matrix} and
221      * allow to optimize the frustum plane extraction in the case when no intersection test is needed for spheres.
222      * <p>
223      * Reference: <a href="http://gamedevs.org/uploads/fast-extraction-viewing-frustum-planes-from-world-view-projection-matrix.pdf">
224      * Fast Extraction of Viewing Frustum Planes from the World-View-Projection Matrix</a>
225      * 
226      * @param m
227      *          the {@link Matrix4d matrix} to update <code>this</code> frustum culler's frustum planes from
228      * @param allowTestSpheres
229      *          whether the methods {@link #testSphere(Vector3d, double)}, {@link #testSphere(double, double, double, double)},
230      *          {@link #intersectSphere(Vector3d, double)} or {@link #intersectSphere(double, double, double, double)} will be used.
231      *          If no spheres need to be tested, then <code>false</code> should be used
232      * @return this
233      */
234     ref public FrustumIntersection set(Matrix4d m, bool allowTestSpheres) return {
235         double invl;
236         nxX = m.m03 + m.m00; nxY = m.m13 + m.m10; nxZ = m.m23 + m.m20; nxW = m.m33 + m.m30;
237         if (allowTestSpheres) {
238             invl = Math.invsqrt(nxX * nxX + nxY * nxY + nxZ * nxZ);
239             nxX *= invl; nxY *= invl; nxZ *= invl; nxW *= invl;
240         }
241         planes[0].set(nxX, nxY, nxZ, nxW);
242         pxX = m.m03 - m.m00; pxY = m.m13 - m.m10; pxZ = m.m23 - m.m20; pxW = m.m33 - m.m30;
243         if (allowTestSpheres) {
244             invl = Math.invsqrt(pxX * pxX + pxY * pxY + pxZ * pxZ);
245             pxX *= invl; pxY *= invl; pxZ *= invl; pxW *= invl;
246         }
247         planes[1].set(pxX, pxY, pxZ, pxW);
248         nyX = m.m03 + m.m01; nyY = m.m13 + m.m11; nyZ = m.m23 + m.m21; nyW = m.m33 + m.m31;
249         if (allowTestSpheres) {
250             invl = Math.invsqrt(nyX * nyX + nyY * nyY + nyZ * nyZ);
251             nyX *= invl; nyY *= invl; nyZ *= invl; nyW *= invl;
252         }
253         planes[2].set(nyX, nyY, nyZ, nyW);
254         pyX = m.m03 - m.m01; pyY = m.m13 - m.m11; pyZ = m.m23 - m.m21; pyW = m.m33 - m.m31;
255         if (allowTestSpheres) {
256             invl = Math.invsqrt(pyX * pyX + pyY * pyY + pyZ * pyZ);
257             pyX *= invl; pyY *= invl; pyZ *= invl; pyW *= invl;
258         }
259         planes[3].set(pyX, pyY, pyZ, pyW);
260         nzX = m.m03 + m.m02; nzY = m.m13 + m.m12; nzZ = m.m23 + m.m22; nzW = m.m33 + m.m32;
261         if (allowTestSpheres) {
262             invl = Math.invsqrt(nzX * nzX + nzY * nzY + nzZ * nzZ);
263             nzX *= invl; nzY *= invl; nzZ *= invl; nzW *= invl;
264         }
265         planes[4].set(nzX, nzY, nzZ, nzW);
266         pzX = m.m03 - m.m02; pzY = m.m13 - m.m12; pzZ = m.m23 - m.m22; pzW = m.m33 - m.m32;
267         if (allowTestSpheres) {
268             invl = Math.invsqrt(pzX * pzX + pzY * pzY + pzZ * pzZ);
269             pzX *= invl; pzY *= invl; pzZ *= invl; pzW *= invl;
270         }
271         planes[5].set(pzX, pzY, pzZ, pzW);
272         return this;
273     }
274 
275     /**
276      * Test whether the given point is within the frustum defined by <code>this</code> frustum culler.
277      * 
278      * @param point
279      *          the point to test
280      * @return <code>true</code> if the given point is inside the frustum; <code>false</code> otherwise
281      */
282     public bool testPoint(Vector3d point) {
283         return testPoint(point.x, point.y, point.z);
284     }
285 
286     /**
287      * Test whether the given point <code>(x, y, z)</code> is within the frustum defined by <code>this</code> frustum culler.
288      * 
289      * @param x
290      *          the x-coordinate of the point
291      * @param y
292      *          the y-coordinate of the point
293      * @param z
294      *          the z-coordinate of the point
295      * @return <code>true</code> if the given point is inside the frustum; <code>false</code> otherwise
296      */
297     public bool testPoint(double x, double y, double z) {
298         return nxX * x + nxY * y + nxZ * z + nxW >= 0 &&
299                pxX * x + pxY * y + pxZ * z + pxW >= 0 &&
300                nyX * x + nyY * y + nyZ * z + nyW >= 0 &&
301                pyX * x + pyY * y + pyZ * z + pyW >= 0 &&
302                nzX * x + nzY * y + nzZ * z + nzW >= 0 &&
303                pzX * x + pzY * y + pzZ * z + pzW >= 0;
304     }
305 
306     /**
307      * Test whether the given sphere is partly or completely within or outside of the frustum defined by <code>this</code> frustum culler.
308      * <p>
309      * The algorithm implemented by this method is conservative. This means that in certain circumstances a <i>false positive</i>
310      * can occur, when the method returns <code>true</code> for spheres that do not intersect the frustum.
311      * See <a href="http://iquilezles.org/www/articles/frustumcorrect/frustumcorrect.htm">iquilezles.org</a> for an examination of this problem.
312      * 
313      * @param center
314      *          the sphere's center
315      * @param radius
316      *          the sphere's radius
317      * @return <code>true</code> if the given sphere is partly or completely inside the frustum;
318      *         <code>false</code> otherwise
319      */
320     public bool testSphere(Vector3d center, double radius) {
321         return testSphere(center.x, center.y, center.z, radius);
322     }
323 
324     /**
325      * Test whether the given sphere is partly or completely within or outside of the frustum defined by <code>this</code> frustum culler.
326      * <p>
327      * The algorithm implemented by this method is conservative. This means that in certain circumstances a <i>false positive</i>
328      * can occur, when the method returns <code>true</code> for spheres that do not intersect the frustum.
329      * See <a href="http://iquilezles.org/www/articles/frustumcorrect/frustumcorrect.htm">iquilezles.org</a> for an examination of this problem.
330      * 
331      * @param x
332      *          the x-coordinate of the sphere's center
333      * @param y
334      *          the y-coordinate of the sphere's center
335      * @param z
336      *          the z-coordinate of the sphere's center
337      * @param r
338      *          the sphere's radius
339      * @return <code>true</code> if the given sphere is partly or completely inside the frustum;
340      *         <code>false</code> otherwise
341      */
342     public bool testSphere(double x, double y, double z, double r) {
343         return nxX * x + nxY * y + nxZ * z + nxW >= -r &&
344                pxX * x + pxY * y + pxZ * z + pxW >= -r &&
345                nyX * x + nyY * y + nyZ * z + nyW >= -r &&
346                pyX * x + pyY * y + pyZ * z + pyW >= -r &&
347                nzX * x + nzY * y + nzZ * z + nzW >= -r &&
348                pzX * x + pzY * y + pzZ * z + pzW >= -r;
349     }
350 
351     /**
352      * Determine whether the given sphere is partly or completely within or outside of the frustum defined by <code>this</code> frustum culler.
353      * <p>
354      * The algorithm implemented by this method is conservative. This means that in certain circumstances a <i>false positive</i>
355      * can occur, when the method returns {@link #INTERSECT} for spheres that do not intersect the frustum.
356      * See <a href="http://iquilezles.org/www/articles/frustumcorrect/frustumcorrect.htm">iquilezles.org</a> for an examination of this problem.
357      * 
358      * @param center
359      *          the sphere's center
360      * @param radius
361      *          the sphere's radius
362      * @return {@link #INSIDE} if the given sphere is completely inside the frustum, or {@link #INTERSECT} if the sphere intersects
363      *         the frustum, or {@link #OUTSIDE} if the sphere is outside of the frustum
364      */
365     public int intersectSphere(Vector3d center, double radius) {
366         return intersectSphere(center.x, center.y, center.z, radius);
367     }
368 
369     /**
370      * Determine whether the given sphere is partly or completely within or outside of the frustum defined by <code>this</code> frustum culler.
371      * <p>
372      * The algorithm implemented by this method is conservative. This means that in certain circumstances a <i>false positive</i>
373      * can occur, when the method returns {@link #INTERSECT} for spheres that do not intersect the frustum.
374      * See <a href="http://iquilezles.org/www/articles/frustumcorrect/frustumcorrect.htm">iquilezles.org</a> for an examination of this problem.
375      * 
376      * @param x
377      *          the x-coordinate of the sphere's center
378      * @param y
379      *          the y-coordinate of the sphere's center
380      * @param z
381      *          the z-coordinate of the sphere's center
382      * @param r
383      *          the sphere's radius
384      * @return {@link #INSIDE} if the given sphere is completely inside the frustum, or {@link #INTERSECT} if the sphere intersects
385      *         the frustum, or {@link #OUTSIDE} if the sphere is outside of the frustum
386      */
387     public int intersectSphere(double x, double y, double z, double r) {
388         bool inside = true;
389         double dist;
390         dist = nxX * x + nxY * y + nxZ * z + nxW;
391         if (dist >= -r) {
392             inside &= dist >= r;
393             dist = pxX * x + pxY * y + pxZ * z + pxW;
394             if (dist >= -r) {
395                 inside &= dist >= r;
396                 dist = nyX * x + nyY * y + nyZ * z + nyW;
397                 if (dist >= -r) {
398                     inside &= dist >= r;
399                     dist = pyX * x + pyY * y + pyZ * z + pyW;
400                     if (dist >= -r) {
401                         inside &= dist >= r;
402                         dist = nzX * x + nzY * y + nzZ * z + nzW;
403                         if (dist >= -r) {
404                             inside &= dist >= r;
405                             dist = pzX * x + pzY * y + pzZ * z + pzW;
406                             if (dist >= -r) {
407                                 inside &= dist >= r;
408                                 return inside ? INSIDE : INTERSECT;
409                             }
410                         }
411                     }
412                 }
413             }
414         }
415         return OUTSIDE;
416     }
417 
418     /**
419      * Test whether the given axis-aligned box is partly or completely within or outside of the frustum defined by <code>this</code> frustum culler.
420      * The box is specified via its <code>min</code> and <code>max</code> corner coordinates.
421      * <p>
422      * The algorithm implemented by this method is conservative. This means that in certain circumstances a <i>false positive</i>
423      * can occur, when the method returns <code>true</code> for boxes that do not intersect the frustum.
424      * See <a href="http://iquilezles.org/www/articles/frustumcorrect/frustumcorrect.htm">iquilezles.org</a> for an examination of this problem.
425      * 
426      * @param min
427      *          the minimum corner coordinates of the axis-aligned box
428      * @param max
429      *          the maximum corner coordinates of the axis-aligned box
430      * @return <code>true</code> if the axis-aligned box is completely or partly inside of the frustum; <code>false</code> otherwise
431      */
432     public bool testAab(Vector3d min, Vector3d max) {
433         return testAab(min.x, min.y, min.z, max.x, max.y, max.z);
434     }
435 
436     /**
437      * Test whether the given axis-aligned box is partly or completely within or outside of the frustum defined by <code>this</code> frustum culler.
438      * The box is specified via its min and max corner coordinates.
439      * <p>
440      * The algorithm implemented by this method is conservative. This means that in certain circumstances a <i>false positive</i>
441      * can occur, when the method returns <code>true</code> for boxes that do not intersect the frustum.
442      * See <a href="http://iquilezles.org/www/articles/frustumcorrect/frustumcorrect.htm">iquilezles.org</a> for an examination of this problem.
443      * <p>
444      * Reference: <a href="http://old.cescg.org/CESCG-2002/DSykoraJJelinek/">Efficient View Frustum Culling</a>
445      * 
446      * @param minX
447      *          the x-coordinate of the minimum corner
448      * @param minY
449      *          the y-coordinate of the minimum corner
450      * @param minZ
451      *          the z-coordinate of the minimum corner
452      * @param maxX
453      *          the x-coordinate of the maximum corner
454      * @param maxY
455      *          the y-coordinate of the maximum corner
456      * @param maxZ
457      *          the z-coordinate of the maximum corner
458      * @return <code>true</code> if the axis-aligned box is completely or partly inside of the frustum; <code>false</code> otherwise
459      */
460     public bool testAab(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
461         /*
462          * This is an implementation of the "2.4 Basic intersection test" of the mentioned site.
463          * It does not distinguish between partially inside and fully inside, though, so the test with the 'p' vertex is omitted.
464          */
465         return nxX * (nxX < 0 ? minX : maxX) + nxY * (nxY < 0 ? minY : maxY) + nxZ * (nxZ < 0 ? minZ : maxZ) >= -nxW &&
466                pxX * (pxX < 0 ? minX : maxX) + pxY * (pxY < 0 ? minY : maxY) + pxZ * (pxZ < 0 ? minZ : maxZ) >= -pxW &&
467                nyX * (nyX < 0 ? minX : maxX) + nyY * (nyY < 0 ? minY : maxY) + nyZ * (nyZ < 0 ? minZ : maxZ) >= -nyW &&
468                pyX * (pyX < 0 ? minX : maxX) + pyY * (pyY < 0 ? minY : maxY) + pyZ * (pyZ < 0 ? minZ : maxZ) >= -pyW &&
469                nzX * (nzX < 0 ? minX : maxX) + nzY * (nzY < 0 ? minY : maxY) + nzZ * (nzZ < 0 ? minZ : maxZ) >= -nzW &&
470                pzX * (pzX < 0 ? minX : maxX) + pzY * (pzY < 0 ? minY : maxY) + pzZ * (pzZ < 0 ? minZ : maxZ) >= -pzW;
471     }
472 
473     /**
474      * Test whether the given XY-plane (at <code>Z = 0</code>) is partly or completely within or outside of the frustum defined by <code>this</code> frustum culler.
475      * The plane is specified via its <code>min</code> and <code>max</code> corner coordinates.
476      * <p>
477      * The algorithm implemented by this method is conservative. This means that in certain circumstances a <i>false positive</i>
478      * can occur, when the method returns <code>true</code> for planes that do not intersect the frustum.
479      * See <a href="http://iquilezles.org/www/articles/frustumcorrect/frustumcorrect.htm">iquilezles.org</a> for an examination of this problem.
480      * 
481      * @param min
482      *          the minimum corner coordinates of the XY-plane
483      * @param max
484      *          the maximum corner coordinates of the XY-plane
485      * @return <code>true</code> if the XY-plane is completely or partly inside of the frustum; <code>false</code> otherwise
486      */
487     public bool testPlaneXY(Vector2d min, Vector2d max) {
488         return testPlaneXY(min.x, min.y, max.x, max.y);
489     }
490 
491     /**
492      * Test whether the given XY-plane (at <code>Z = 0</code>) is partly or completely within or outside of the frustum defined by <code>this</code> frustum culler.
493      * The plane is specified via its min and max corner coordinates.
494      * <p>
495      * The algorithm implemented by this method is conservative. This means that in certain circumstances a <i>false positive</i>
496      * can occur, when the method returns <code>true</code> for planes that do not intersect the frustum.
497      * See <a href="http://iquilezles.org/www/articles/frustumcorrect/frustumcorrect.htm">iquilezles.org</a> for an examination of this problem.
498      * <p>
499      * Reference: <a href="http://old.cescg.org/CESCG-2002/DSykoraJJelinek/">Efficient View Frustum Culling</a>
500      * 
501      * @param minX
502      *          the x-coordinate of the minimum corner
503      * @param minY
504      *          the y-coordinate of the minimum corner
505      * @param maxX
506      *          the x-coordinate of the maximum corner
507      * @param maxY
508      *          the y-coordinate of the maximum corner
509      * @return <code>true</code> if the XY-plane is completely or partly inside of the frustum; <code>false</code> otherwise
510      */
511     public bool testPlaneXY(double minX, double minY, double maxX, double maxY) {
512         /*
513          * This is an implementation of the "2.4 Basic intersection test" of the mentioned site.
514          * It does not distinguish between partially inside and fully inside, though, so the test with the 'p' vertex is omitted.
515          */
516         return nxX * (nxX < 0 ? minX : maxX) + nxY * (nxY < 0 ? minY : maxY) >= -nxW &&
517                pxX * (pxX < 0 ? minX : maxX) + pxY * (pxY < 0 ? minY : maxY) >= -pxW &&
518                nyX * (nyX < 0 ? minX : maxX) + nyY * (nyY < 0 ? minY : maxY) >= -nyW &&
519                pyX * (pyX < 0 ? minX : maxX) + pyY * (pyY < 0 ? minY : maxY) >= -pyW &&
520                nzX * (nzX < 0 ? minX : maxX) + nzY * (nzY < 0 ? minY : maxY) >= -nzW &&
521                pzX * (pzX < 0 ? minX : maxX) + pzY * (pzY < 0 ? minY : maxY) >= -pzW;
522     }
523 
524     /**
525      * Test whether the given XZ-plane (at <code>Y = 0</code>) is partly or completely within or outside of the frustum defined by <code>this</code> frustum culler.
526      * The plane is specified via its min and max corner coordinates.
527      * <p>
528      * The algorithm implemented by this method is conservative. This means that in certain circumstances a <i>false positive</i>
529      * can occur, when the method returns <code>true</code> for planes that do not intersect the frustum.
530      * See <a href="http://iquilezles.org/www/articles/frustumcorrect/frustumcorrect.htm">iquilezles.org</a> for an examination of this problem.
531      * <p>
532      * Reference: <a href="http://old.cescg.org/CESCG-2002/DSykoraJJelinek/">Efficient View Frustum Culling</a>
533      * 
534      * @param minX
535      *          the x-coordinate of the minimum corner
536      * @param minZ
537      *          the z-coordinate of the minimum corner
538      * @param maxX
539      *          the x-coordinate of the maximum corner
540      * @param maxZ
541      *          the z-coordinate of the maximum corner
542      * @return <code>true</code> if the XZ-plane is completely or partly inside of the frustum; <code>false</code> otherwise
543      */
544     public bool testPlaneXZ(double minX, double minZ, double maxX, double maxZ) {
545         /*
546          * This is an implementation of the "2.4 Basic intersection test" of the mentioned site.
547          * It does not distinguish between partially inside and fully inside, though, so the test with the 'p' vertex is omitted.
548          */
549         return nxX * (nxX < 0 ? minX : maxX) + nxZ * (nxZ < 0 ? minZ : maxZ) >= -nxW &&
550                pxX * (pxX < 0 ? minX : maxX) + pxZ * (pxZ < 0 ? minZ : maxZ) >= -pxW &&
551                nyX * (nyX < 0 ? minX : maxX) + nyZ * (nyZ < 0 ? minZ : maxZ) >= -nyW &&
552                pyX * (pyX < 0 ? minX : maxX) + pyZ * (pyZ < 0 ? minZ : maxZ) >= -pyW &&
553                nzX * (nzX < 0 ? minX : maxX) + nzZ * (nzZ < 0 ? minZ : maxZ) >= -nzW &&
554                pzX * (pzX < 0 ? minX : maxX) + pzZ * (pzZ < 0 ? minZ : maxZ) >= -pzW;
555     }
556 
557     /**
558      * Determine whether the given axis-aligned box is partly or completely within or outside of the frustum defined by <code>this</code> frustum culler
559      * and, if the box is not inside this frustum, return the index of the plane that culled it.
560      * The box is specified via its <code>min</code> and <code>max</code> corner coordinates.
561      * <p>
562      * The algorithm implemented by this method is conservative. This means that in certain circumstances a <i>false positive</i>
563      * can occur, when the method returns {@link #INTERSECT} for boxes that do not intersect the frustum.
564      * See <a href="http://iquilezles.org/www/articles/frustumcorrect/frustumcorrect.htm">iquilezles.org</a> for an examination of this problem.
565      * 
566      * @param min
567      *          the minimum corner coordinates of the axis-aligned box
568      * @param max
569      *          the maximum corner coordinates of the axis-aligned box
570      * @return the index of the first plane that culled the box, if the box does not intersect the frustum;
571      *         or {@link #INTERSECT} if the box intersects the frustum, or {@link #INSIDE} if the box is fully inside of the frustum.
572      *         The plane index is one of {@link #PLANE_NX}, {@link #PLANE_PX}, {@link #PLANE_NY}, {@link #PLANE_PY}, {@link #PLANE_NZ} and {@link #PLANE_PZ}
573      */
574     public int intersectAab(Vector3d min, Vector3d max) {
575         return intersectAab(min.x, min.y, min.z, max.x, max.y, max.z);
576     }
577 
578     /**
579      * Determine whether the given axis-aligned box is partly or completely within or outside of the frustum defined by <code>this</code> frustum culler
580      * and, if the box is not inside this frustum, return the index of the plane that culled it.
581      * The box is specified via its min and max corner coordinates.
582      * <p>
583      * The algorithm implemented by this method is conservative. This means that in certain circumstances a <i>false positive</i>
584      * can occur, when the method returns {@link #INTERSECT} for boxes that do not intersect the frustum.
585      * See <a href="http://iquilezles.org/www/articles/frustumcorrect/frustumcorrect.htm">iquilezles.org</a> for an examination of this problem.
586      * <p>
587      * Reference: <a href="http://old.cescg.org/CESCG-2002/DSykoraJJelinek/">Efficient View Frustum Culling</a>
588      * 
589      * @param minX
590      *          the x-coordinate of the minimum corner
591      * @param minY
592      *          the y-coordinate of the minimum corner
593      * @param minZ
594      *          the z-coordinate of the minimum corner
595      * @param maxX
596      *          the x-coordinate of the maximum corner
597      * @param maxY
598      *          the y-coordinate of the maximum corner
599      * @param maxZ
600      *          the z-coordinate of the maximum corner
601      * @return the index of the first plane that culled the box, if the box does not intersect the frustum,
602      *         or {@link #INTERSECT} if the box intersects the frustum, or {@link #INSIDE} if the box is fully inside of the frustum.
603      *         The plane index is one of {@link #PLANE_NX}, {@link #PLANE_PX}, {@link #PLANE_NY}, {@link #PLANE_PY}, {@link #PLANE_NZ} and {@link #PLANE_PZ}
604      */
605     public int intersectAab(double minX, double minY, double minZ, double maxX, double maxY, double maxZ) {
606         /*
607          * This is an implementation of the "2.4 Basic intersection test" of the mentioned site.
608          * 
609          * In addition to the algorithm in the paper, this method also returns the index of the first plane that culled the box.
610          */
611         int plane = PLANE_NX;
612         bool inside = true;
613         if (nxX * (nxX < 0 ? minX : maxX) + nxY * (nxY < 0 ? minY : maxY) + nxZ * (nxZ < 0 ? minZ : maxZ) >= -nxW) {
614             plane = PLANE_PX;
615             inside &= nxX * (nxX < 0 ? maxX : minX) + nxY * (nxY < 0 ? maxY : minY) + nxZ * (nxZ < 0 ? maxZ : minZ) >= -nxW;
616             if (pxX * (pxX < 0 ? minX : maxX) + pxY * (pxY < 0 ? minY : maxY) + pxZ * (pxZ < 0 ? minZ : maxZ) >= -pxW) {
617                 plane = PLANE_NY;
618                 inside &= pxX * (pxX < 0 ? maxX : minX) + pxY * (pxY < 0 ? maxY : minY) + pxZ * (pxZ < 0 ? maxZ : minZ) >= -pxW;
619                 if (nyX * (nyX < 0 ? minX : maxX) + nyY * (nyY < 0 ? minY : maxY) + nyZ * (nyZ < 0 ? minZ : maxZ) >= -nyW) {
620                     plane = PLANE_PY;
621                     inside &= nyX * (nyX < 0 ? maxX : minX) + nyY * (nyY < 0 ? maxY : minY) + nyZ * (nyZ < 0 ? maxZ : minZ) >= -nyW;
622                     if (pyX * (pyX < 0 ? minX : maxX) + pyY * (pyY < 0 ? minY : maxY) + pyZ * (pyZ < 0 ? minZ : maxZ) >= -pyW) {
623                         plane = PLANE_NZ;
624                         inside &= pyX * (pyX < 0 ? maxX : minX) + pyY * (pyY < 0 ? maxY : minY) + pyZ * (pyZ < 0 ? maxZ : minZ) >= -pyW;
625                         if (nzX * (nzX < 0 ? minX : maxX) + nzY * (nzY < 0 ? minY : maxY) + nzZ * (nzZ < 0 ? minZ : maxZ) >= -nzW) {
626                             plane = PLANE_PZ;
627                             inside &= nzX * (nzX < 0 ? maxX : minX) + nzY * (nzY < 0 ? maxY : minY) + nzZ * (nzZ < 0 ? maxZ : minZ) >= -nzW;
628                             if (pzX * (pzX < 0 ? minX : maxX) + pzY * (pzY < 0 ? minY : maxY) + pzZ * (pzZ < 0 ? minZ : maxZ) >= -pzW) {
629                                 inside &= pzX * (pzX < 0 ? maxX : minX) + pzY * (pzY < 0 ? maxY : minY) + pzZ * (pzZ < 0 ? maxZ : minZ) >= -pzW;
630                                 return inside ? INSIDE : INTERSECT;
631                             }
632                         }
633                     }
634                 }
635             }
636         }
637         return plane;
638     }
639 
640     /**
641      * Compute the signed distance from the given axis-aligned box to the <code>plane</code>.
642      * 
643      * @param minX
644      *          the x-coordinate of the minimum corner
645      * @param minY
646      *          the y-coordinate of the minimum corner
647      * @param minZ
648      *          the z-coordinate of the minimum corner
649      * @param maxX
650      *          the x-coordinate of the maximum corner
651      * @param maxY
652      *          the y-coordinate of the maximum corner
653      * @param maxZ
654      *          the z-coordinate of the maximum corner
655      * @param plane
656      *          one of 
657      *          {@link #PLANE_NX}, {@link #PLANE_PX},
658      *          {@link #PLANE_NY}, {@link #PLANE_PY}, 
659      *          {@link #PLANE_NZ} and {@link #PLANE_PZ}
660      * @return the signed distance of the axis-aligned box to the plane
661      */
662     public double distanceToPlane(double minX, double minY, double minZ, double maxX, double maxY, double maxZ, int plane) {
663         return planes[plane].x * (planes[plane].x < 0 ? maxX : minX) + planes[plane].y * (planes[plane].y < 0 ? maxY : minY)
664                 + planes[plane].z * (planes[plane].z < 0 ? maxZ : minZ) + planes[plane].w;
665     }
666 
667     /**
668      * Determine whether the given axis-aligned box is partly or completely within or outside of the frustum defined by <code>this</code> frustum culler
669      * and, if the box is not inside this frustum, return the index of the plane that culled it.
670      * The box is specified via its <code>min</code> and <code>max</code> corner coordinates.
671      * <p>
672      * This method differs from {@link #intersectAab(Vector3d, Vector3d)} in that
673      * it allows to mask-off planes that should not be calculated. For example, in order to only test a box against the
674      * left frustum plane, use a mask of {@link #PLANE_MASK_NX}. Or in order to test all planes <i>except</i> the left plane, use 
675      * a mask of <code>(~0 ^ PLANE_MASK_NX)</code>.
676      * <p>
677      * The algorithm implemented by this method is conservative. This means that in certain circumstances a <i>false positive</i>
678      * can occur, when the method returns {@link #INTERSECT} for boxes that do not intersect the frustum.
679      * See <a href="http://iquilezles.org/www/articles/frustumcorrect/frustumcorrect.htm">iquilezles.org</a> for an examination of this problem.
680      * 
681      * @param min
682      *          the minimum corner coordinates of the axis-aligned box
683      * @param max
684      *          the maximum corner coordinates of the axis-aligned box
685      * @param mask
686      *          contains as bitset all the planes that should be tested.
687      *          This value can be any combination of 
688      *          {@link #PLANE_MASK_NX}, {@link #PLANE_MASK_PX},
689      *          {@link #PLANE_MASK_NY}, {@link #PLANE_MASK_PY}, 
690      *          {@link #PLANE_MASK_NZ} and {@link #PLANE_MASK_PZ}
691      * @return the index of the first plane that culled the box, if the box does not intersect the frustum,
692      *         or {@link #INTERSECT} if the box intersects the frustum, or {@link #INSIDE} if the box is fully inside of the frustum.
693      *         The plane index is one of {@link #PLANE_NX}, {@link #PLANE_PX}, {@link #PLANE_NY}, {@link #PLANE_PY}, {@link #PLANE_NZ} and {@link #PLANE_PZ}
694      */
695     public int intersectAab(Vector3d min, Vector3d max, int mask) {
696         return intersectAab(min.x, min.y, min.z, max.x, max.y, max.z, mask);
697     }
698 
699     /**
700      * Determine whether the given axis-aligned box is partly or completely within or outside of the frustum defined by <code>this</code> frustum culler
701      * and, if the box is not inside this frustum, return the index of the plane that culled it.
702      * The box is specified via its min and max corner coordinates.
703      * <p>
704      * This method differs from {@link #intersectAab(double, double, double, double, double, double)} in that
705      * it allows to mask-off planes that should not be calculated. For example, in order to only test a box against the
706      * left frustum plane, use a mask of {@link #PLANE_MASK_NX}. Or in order to test all planes <i>except</i> the left plane, use 
707      * a mask of <code>(~0 ^ PLANE_MASK_NX)</code>.
708      * <p>
709      * The algorithm implemented by this method is conservative. This means that in certain circumstances a <i>false positive</i>
710      * can occur, when the method returns {@link #INTERSECT} for boxes that do not intersect the frustum.
711      * See <a href="http://iquilezles.org/www/articles/frustumcorrect/frustumcorrect.htm">iquilezles.org</a> for an examination of this problem.
712      * <p>
713      * Reference: <a href="http://old.cescg.org/CESCG-2002/DSykoraJJelinek/">Efficient View Frustum Culling</a>
714      * 
715      * @param minX
716      *          the x-coordinate of the minimum corner
717      * @param minY
718      *          the y-coordinate of the minimum corner
719      * @param minZ
720      *          the z-coordinate of the minimum corner
721      * @param maxX
722      *          the x-coordinate of the maximum corner
723      * @param maxY
724      *          the y-coordinate of the maximum corner
725      * @param maxZ
726      *          the z-coordinate of the maximum corner
727      * @param mask
728      *          contains as bitset all the planes that should be tested.
729      *          This value can be any combination of 
730      *          {@link #PLANE_MASK_NX}, {@link #PLANE_MASK_PX},
731      *          {@link #PLANE_MASK_NY}, {@link #PLANE_MASK_PY}, 
732      *          {@link #PLANE_MASK_NZ} and {@link #PLANE_MASK_PZ}
733      * @return the index of the first plane that culled the box, if the box does not intersect the frustum,
734      *         or {@link #INTERSECT} if the box intersects the frustum, or {@link #INSIDE} if the box is fully inside of the frustum.
735      *         The plane index is one of {@link #PLANE_NX}, {@link #PLANE_PX}, {@link #PLANE_NY}, {@link #PLANE_PY}, {@link #PLANE_NZ} and {@link #PLANE_PZ}
736      */
737     public int intersectAab(double minX, double minY, double minZ, double maxX, double maxY, double maxZ, int mask) {
738         /*
739          * This is an implementation of the first algorithm in "2.5 Plane masking and coherency" of the mentioned site.
740          * 
741          * In addition to the algorithm in the paper, this method also returns the index of the first plane that culled the box.
742          */
743         int plane = PLANE_NX;
744         bool inside = true;
745         if ((mask & PLANE_MASK_NX) == 0 || nxX * (nxX < 0 ? minX : maxX) + nxY * (nxY < 0 ? minY : maxY) + nxZ * (nxZ < 0 ? minZ : maxZ) >= -nxW) {
746             plane = PLANE_PX;
747             inside &= nxX * (nxX < 0 ? maxX : minX) + nxY * (nxY < 0 ? maxY : minY) + nxZ * (nxZ < 0 ? maxZ : minZ) >= -nxW;
748             if ((mask & PLANE_MASK_PX) == 0 || pxX * (pxX < 0 ? minX : maxX) + pxY * (pxY < 0 ? minY : maxY) + pxZ * (pxZ < 0 ? minZ : maxZ) >= -pxW) {
749                 plane = PLANE_NY;
750                 inside &= pxX * (pxX < 0 ? maxX : minX) + pxY * (pxY < 0 ? maxY : minY) + pxZ * (pxZ < 0 ? maxZ : minZ) >= -pxW;
751                 if ((mask & PLANE_MASK_NY) == 0 || nyX * (nyX < 0 ? minX : maxX) + nyY * (nyY < 0 ? minY : maxY) + nyZ * (nyZ < 0 ? minZ : maxZ) >= -nyW) {
752                     plane = PLANE_PY;
753                     inside &= nyX * (nyX < 0 ? maxX : minX) + nyY * (nyY < 0 ? maxY : minY) + nyZ * (nyZ < 0 ? maxZ : minZ) >= -nyW;
754                     if ((mask & PLANE_MASK_PY) == 0 || pyX * (pyX < 0 ? minX : maxX) + pyY * (pyY < 0 ? minY : maxY) + pyZ * (pyZ < 0 ? minZ : maxZ) >= -pyW) {
755                         plane = PLANE_NZ;
756                         inside &= pyX * (pyX < 0 ? maxX : minX) + pyY * (pyY < 0 ? maxY : minY) + pyZ * (pyZ < 0 ? maxZ : minZ) >= -pyW;
757                         if ((mask & PLANE_MASK_NZ) == 0 || nzX * (nzX < 0 ? minX : maxX) + nzY * (nzY < 0 ? minY : maxY) + nzZ * (nzZ < 0 ? minZ : maxZ) >= -nzW) {
758                             plane = PLANE_PZ;
759                             inside &= nzX * (nzX < 0 ? maxX : minX) + nzY * (nzY < 0 ? maxY : minY) + nzZ * (nzZ < 0 ? maxZ : minZ) >= -nzW;
760                             if ((mask & PLANE_MASK_PZ) == 0 || pzX * (pzX < 0 ? minX : maxX) + pzY * (pzY < 0 ? minY : maxY) + pzZ * (pzZ < 0 ? minZ : maxZ) >= -pzW) {
761                                 inside &= pzX * (pzX < 0 ? maxX : minX) + pzY * (pzY < 0 ? maxY : minY) + pzZ * (pzZ < 0 ? maxZ : minZ) >= -pzW;
762                                 return inside ? INSIDE : INTERSECT;
763                             }
764                         }
765                     }
766                 }
767             }
768         }
769         return plane;
770     }
771 
772     /**
773      * Determine whether the given axis-aligned box is partly or completely within or outside of the frustum defined by <code>this</code> frustum culler
774      * and, if the box is not inside this frustum, return the index of the plane that culled it.
775      * The box is specified via its <code>min</code> and <code>max</code> corner coordinates.
776      * <p>
777      * This method differs from {@link #intersectAab(Vector3d, Vector3d)} in that
778      * it allows to mask-off planes that should not be calculated. For example, in order to only test a box against the
779      * left frustum plane, use a mask of {@link #PLANE_MASK_NX}. Or in order to test all planes <i>except</i> the left plane, use 
780      * a mask of <code>(~0 ^ PLANE_MASK_NX)</code>.
781      * <p>
782      * In addition, the <code>startPlane</code> denotes the first frustum plane to test the box against. To use this effectively means to store the
783      * plane that previously culled an axis-aligned box (as returned by <code>intersectAab()</code>) and in the next frame use the return value
784      * as the argument to the <code>startPlane</code> parameter of this method. The assumption is that the plane that culled the object previously will also
785      * cull it now (temporal coherency) and the culling computation is likely reduced in that case.
786      * <p>
787      * The algorithm implemented by this method is conservative. This means that in certain circumstances a <i>false positive</i>
788      * can occur, when the method returns {@link #INTERSECT} for boxes that do not intersect the frustum.
789      * See <a href="http://iquilezles.org/www/articles/frustumcorrect/frustumcorrect.htm">iquilezles.org</a> for an examination of this problem.
790      * 
791      * @param min
792      *          the minimum corner coordinates of the axis-aligned box
793      * @param max
794      *          the maximum corner coordinates of the axis-aligned box
795      * @param mask
796      *          contains as bitset all the planes that should be tested.
797      *          This value can be any combination of 
798      *          {@link #PLANE_MASK_NX}, {@link #PLANE_MASK_PX},
799      *          {@link #PLANE_MASK_NY}, {@link #PLANE_MASK_PY}, 
800      *          {@link #PLANE_MASK_NZ} and {@link #PLANE_MASK_PZ}
801      * @param startPlane
802      *          the first frustum plane to test the axis-aligned box against. It is one of
803      *          {@link #PLANE_NX}, {@link #PLANE_PX}, {@link #PLANE_NY}, {@link #PLANE_PY}, {@link #PLANE_NZ} and {@link #PLANE_PZ}
804      * @return the index of the first plane that culled the box, if the box does not intersect the frustum,
805      *         or {@link #INTERSECT} if the box intersects the frustum, or {@link #INSIDE} if the box is fully inside of the frustum.
806      *         The plane index is one of {@link #PLANE_NX}, {@link #PLANE_PX}, {@link #PLANE_NY}, {@link #PLANE_PY}, {@link #PLANE_NZ} and {@link #PLANE_PZ}
807      */
808     public int intersectAab(Vector3d min, Vector3d max, int mask, int startPlane) {
809         return intersectAab(min.x, min.y, min.z, max.x, max.y, max.z, mask, startPlane);
810     }
811 
812     /**
813      * Determine whether the given axis-aligned box is partly or completely within or outside of the frustum defined by <code>this</code> frustum culler
814      * and, if the box is not inside this frustum, return the index of the plane that culled it.
815      * The box is specified via its min and max corner coordinates.
816      * <p>
817      * This method differs from {@link #intersectAab(double, double, double, double, double, double)} in that
818      * it allows to mask-off planes that should not be calculated. For example, in order to only test a box against the
819      * left frustum plane, use a mask of {@link #PLANE_MASK_NX}. Or in order to test all planes <i>except</i> the left plane, use 
820      * a mask of <code>(~0 ^ PLANE_MASK_NX)</code>.
821      * <p>
822      * In addition, the <code>startPlane</code> denotes the first frustum plane to test the box against. To use this effectively means to store the
823      * plane that previously culled an axis-aligned box (as returned by <code>intersectAab()</code>) and in the next frame use the return value
824      * as the argument to the <code>startPlane</code> parameter of this method. The assumption is that the plane that culled the object previously will also
825      * cull it now (temporal coherency) and the culling computation is likely reduced in that case.
826      * <p>
827      * The algorithm implemented by this method is conservative. This means that in certain circumstances a <i>false positive</i>
828      * can occur, when the method returns {@link #INTERSECT} for boxes that do not intersect the frustum.
829      * See <a href="http://iquilezles.org/www/articles/frustumcorrect/frustumcorrect.htm">iquilezles.org</a> for an examination of this problem.
830      * <p>
831      * Reference: <a href="http://old.cescg.org/CESCG-2002/DSykoraJJelinek/">Efficient View Frustum Culling</a>
832      * 
833      * @param minX
834      *          the x-coordinate of the minimum corner
835      * @param minY
836      *          the y-coordinate of the minimum corner
837      * @param minZ
838      *          the z-coordinate of the minimum corner
839      * @param maxX
840      *          the x-coordinate of the maximum corner
841      * @param maxY
842      *          the y-coordinate of the maximum corner
843      * @param maxZ
844      *          the z-coordinate of the maximum corner
845      * @param mask
846      *          contains as bitset all the planes that should be tested.
847      *          This value can be any combination of 
848      *          {@link #PLANE_MASK_NX}, {@link #PLANE_MASK_PX},
849      *          {@link #PLANE_MASK_NY}, {@link #PLANE_MASK_PY}, 
850      *          {@link #PLANE_MASK_NZ} and {@link #PLANE_MASK_PZ}
851      * @param startPlane
852      *          the first frustum plane to test the axis-aligned box against. It is one of
853      *          {@link #PLANE_NX}, {@link #PLANE_PX}, {@link #PLANE_NY}, {@link #PLANE_PY}, {@link #PLANE_NZ} and {@link #PLANE_PZ}
854      * @return the index of the first plane that culled the box, if the box does not intersect the frustum,
855      *         or {@link #INTERSECT} if the box intersects the frustum, or {@link #INSIDE} if the box is fully inside of the frustum.
856      *         The plane index is one of {@link #PLANE_NX}, {@link #PLANE_PX}, {@link #PLANE_NY}, {@link #PLANE_PY}, {@link #PLANE_NZ} and {@link #PLANE_PZ}
857      */
858     public int intersectAab(double minX, double minY, double minZ, double maxX, double maxY, double maxZ, int mask, int startPlane) {
859         /*
860          * This is an implementation of the second algorithm in "2.5 Plane masking and coherency" of the mentioned site.
861          * 
862          * In addition to the algorithm in the paper, this method also returns the index of the first plane that culled the box.
863          */
864         int plane = startPlane;
865         bool inside = true;
866         Vector4d p = planes[startPlane];
867         if ((mask & 1<<startPlane) != 0 && p.x * (p.x < 0 ? minX : maxX) + p.y * (p.y < 0 ? minY : maxY) + p.z * (p.z < 0 ? minZ : maxZ) < -p.w) {
868             return plane;
869         }
870         if ((mask & PLANE_MASK_NX) == 0 || nxX * (nxX < 0 ? minX : maxX) + nxY * (nxY < 0 ? minY : maxY) + nxZ * (nxZ < 0 ? minZ : maxZ) >= -nxW) {
871             plane = PLANE_PX;
872             inside &= nxX * (nxX < 0 ? maxX : minX) + nxY * (nxY < 0 ? maxY : minY) + nxZ * (nxZ < 0 ? maxZ : minZ) >= -nxW;
873             if ((mask & PLANE_MASK_PX) == 0 || pxX * (pxX < 0 ? minX : maxX) + pxY * (pxY < 0 ? minY : maxY) + pxZ * (pxZ < 0 ? minZ : maxZ) >= -pxW) {
874                 plane = PLANE_NY;
875                 inside &= pxX * (pxX < 0 ? maxX : minX) + pxY * (pxY < 0 ? maxY : minY) + pxZ * (pxZ < 0 ? maxZ : minZ) >= -pxW;
876                 if ((mask & PLANE_MASK_NY) == 0 || nyX * (nyX < 0 ? minX : maxX) + nyY * (nyY < 0 ? minY : maxY) + nyZ * (nyZ < 0 ? minZ : maxZ) >= -nyW) {
877                     plane = PLANE_PY;
878                     inside &= nyX * (nyX < 0 ? maxX : minX) + nyY * (nyY < 0 ? maxY : minY) + nyZ * (nyZ < 0 ? maxZ : minZ) >= -nyW;
879                     if ((mask & PLANE_MASK_PY) == 0 || pyX * (pyX < 0 ? minX : maxX) + pyY * (pyY < 0 ? minY : maxY) + pyZ * (pyZ < 0 ? minZ : maxZ) >= -pyW) {
880                         plane = PLANE_NZ;
881                         inside &= pyX * (pyX < 0 ? maxX : minX) + pyY * (pyY < 0 ? maxY : minY) + pyZ * (pyZ < 0 ? maxZ : minZ) >= -pyW;
882                         if ((mask & PLANE_MASK_NZ) == 0 || nzX * (nzX < 0 ? minX : maxX) + nzY * (nzY < 0 ? minY : maxY) + nzZ * (nzZ < 0 ? minZ : maxZ) >= -nzW) {
883                             plane = PLANE_PZ;
884                             inside &= nzX * (nzX < 0 ? maxX : minX) + nzY * (nzY < 0 ? maxY : minY) + nzZ * (nzZ < 0 ? maxZ : minZ) >= -nzW;
885                             if ((mask & PLANE_MASK_PZ) == 0 || pzX * (pzX < 0 ? minX : maxX) + pzY * (pzY < 0 ? minY : maxY) + pzZ * (pzZ < 0 ? minZ : maxZ) >= -pzW) {
886                                 inside &= pzX * (pzX < 0 ? maxX : minX) + pzY * (pzY < 0 ? maxY : minY) + pzZ * (pzZ < 0 ? maxZ : minZ) >= -pzW;
887                                 return inside ? INSIDE : INTERSECT;
888                             }
889                         }
890                     }
891                 }
892             }
893         }
894         return plane;
895     }
896 
897     /**
898      * Test whether the given line segment, defined by the end points <code>a</code> and <code>b</code>, 
899      * is partly or completely within the frustum defined by <code>this</code> frustum culler.
900      * 
901      * @param a
902      *          the line segment's first end point
903      * @param b
904      *          the line segment's second end point
905      * @return <code>true</code> if the given line segment is partly or completely inside the frustum;
906      *         <code>false</code> otherwise
907      */
908     public bool testLineSegment(Vector3d a, Vector3d b) {
909         return testLineSegment(a.x, a.y, a.z, b.x, b.y, b.z);
910     }
911 
912     /**
913      * Test whether the given line segment, defined by the end points <code>(aX, aY, aZ)</code> and <code>(bX, bY, bZ)</code>, 
914      * is partly or completely within the frustum defined by <code>this</code> frustum culler.
915      * 
916      * @param aX
917      *          the x coordinate of the line segment's first end point
918      * @param aY
919      *          the y coordinate of the line segment's first end point
920      * @param aZ
921      *          the z coordinate of the line segment's first end point
922      * @param bX
923      *          the x coordinate of the line segment's second end point
924      * @param bY
925      *          the y coordinate of the line segment's second end point
926      * @param bZ
927      *          the z coordinate of the line segment's second end point
928      * @return <code>true</code> if the given line segment is partly or completely inside the frustum;
929      *         <code>false</code> otherwise
930      */
931     public bool testLineSegment(double aX, double aY, double aZ, double bX, double bY, double bZ) {
932         double da, db;
933         da = Math.fma(nxX, aX, Math.fma(nxY, aY, Math.fma(nxZ, aZ, nxW)));
934         db = Math.fma(nxX, bX, Math.fma(nxY, bY, Math.fma(nxZ, bZ, nxW)));
935         if (da < 0.0f && db < 0.0f)
936             return false;
937         if (da * db < 0.0f) {
938             double p = Math.abs(da) / Math.abs(db - da);
939             double dx = Math.fma(bX - aX, p, aX), dy = Math.fma(bY - aY, p, aY), dz = Math.fma(bZ - aZ, p, aZ);
940             if (da < 0.0f) {
941                 aX = dx; aY = dy; aZ = dz;
942             } else {
943                 bX = dx; bY = dy; bZ = dz;
944             }
945         }
946         da = Math.fma(pxX, aX, Math.fma(pxY, aY, Math.fma(pxZ, aZ, pxW)));
947         db = Math.fma(pxX, bX, Math.fma(pxY, bY, Math.fma(pxZ, bZ, pxW)));
948         if (da < 0.0f && db < 0.0f)
949             return false;
950         if (da * db < 0.0f) {
951             double p = Math.abs(da) / Math.abs(db - da);
952             double dx = Math.fma(bX - aX, p, aX), dy = Math.fma(bY - aY, p, aY), dz = Math.fma(bZ - aZ, p, aZ);
953             if (da < 0.0f) {
954                 aX = dx; aY = dy; aZ = dz;
955             } else {
956                 bX = dx; bY = dy; bZ = dz;
957             }
958         }
959         da = Math.fma(nyX, aX, Math.fma(nyY, aY, Math.fma(nyZ, aZ, nyW)));
960         db = Math.fma(nyX, bX, Math.fma(nyY, bY, Math.fma(nyZ, bZ, nyW)));
961         if (da < 0.0f && db < 0.0f)
962             return false;
963         if (da * db < 0.0f) {
964             double p = Math.abs(da) / Math.abs(db - da);
965             double dx = Math.fma(bX - aX, p, aX), dy = Math.fma(bY - aY, p, aY), dz = Math.fma(bZ - aZ, p, aZ);
966             if (da < 0.0f) {
967                 aX = dx; aY = dy; aZ = dz;
968             } else {
969                 bX = dx; bY = dy; bZ = dz;
970             }
971         }
972         da = Math.fma(pyX, aX, Math.fma(pyY, aY, Math.fma(pyZ, aZ, pyW)));
973         db = Math.fma(pyX, bX, Math.fma(pyY, bY, Math.fma(pyZ, bZ, pyW)));
974         if (da < 0.0f && db < 0.0f)
975             return false;
976         if (da * db < 0.0f) {
977             double p = Math.abs(da) / Math.abs(db - da);
978             double dx = Math.fma(bX - aX, p, aX), dy = Math.fma(bY - aY, p, aY), dz = Math.fma(bZ - aZ, p, aZ);
979             if (da < 0.0f) {
980                 aX = dx; aY = dy; aZ = dz;
981             } else {
982                 bX = dx; bY = dy; bZ = dz;
983             }
984         }
985         da = Math.fma(nzX, aX, Math.fma(nzY, aY, Math.fma(nzZ, aZ, nzW)));
986         db = Math.fma(nzX, bX, Math.fma(nzY, bY, Math.fma(nzZ, bZ, nzW)));
987         if (da < 0.0f && db < 0.0f)
988             return false;
989         if (da * db < 0.0f) {
990             double p = Math.abs(da) / Math.abs(db - da);
991             double dx = Math.fma(bX - aX, p, aX), dy = Math.fma(bY - aY, p, aY), dz = Math.fma(bZ - aZ, p, aZ);
992             if (da < 0.0f) {
993                 aX = dx; aY = dy; aZ = dz;
994             } else {
995                 bX = dx; bY = dy; bZ = dz;
996             }
997         }
998         da = Math.fma(pzX, aX, Math.fma(pzY, aY, Math.fma(pzZ, aZ, pzW)));
999         db = Math.fma(pzX, bX, Math.fma(pzY, bY, Math.fma(pzZ, bZ, pzW)));
1000         return da >= 0.0f || db >= 0.0f;
1001     }
1002 
1003 }