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