1 // Copyright (C) 2006-2009 Mathias Froehlich - Mathias.Froehlich@web.de
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Library General Public
5 // License as published by the Free Software Foundation; either
6 // version 2 of the License, or (at your option) any later version.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Library General Public License for more details.
13 // You should have received a copy of the GNU General Public License
14 // along with this program; if not, write to the Free Software
15 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 #ifndef SGIntersect_HXX
19 #define SGIntersect_HXX
23 intersects(const SGSphere<T>& s1, const SGSphere<T>& s2)
30 T dist = s1.getRadius() + s2.getRadius();
31 return distSqr(s1.getCenter(), s2.getCenter()) <= dist*dist;
35 template<typename T1, typename T2>
37 intersects(const SGBox<T1>& box, const SGSphere<T2>& sphere)
41 // Is more or less trivially included in the next tests
45 if (sphere.getCenter().x() < box.getMin().x() - sphere.getRadius())
47 if (sphere.getCenter().y() < box.getMin().y() - sphere.getRadius())
49 if (sphere.getCenter().z() < box.getMin().z() - sphere.getRadius())
52 if (box.getMax().x() + sphere.getRadius() < sphere.getCenter().x())
54 if (box.getMax().y() + sphere.getRadius() < sphere.getCenter().y())
56 if (box.getMax().z() + sphere.getRadius() < sphere.getCenter().z())
62 template<typename T1, typename T2>
64 intersects(const SGSphere<T1>& sphere, const SGBox<T2>& box)
65 { return intersects(box, sphere); }
68 template<typename T1, typename T2>
70 intersects(const SGVec3<T1>& v, const SGBox<T2>& box)
72 if (v[0] < box.getMin()[0])
74 if (box.getMax()[0] < v[0])
76 if (v[1] < box.getMin()[1])
78 if (box.getMax()[1] < v[1])
80 if (v[2] < box.getMin()[2])
82 if (box.getMax()[2] < v[2])
86 template<typename T1, typename T2>
88 intersects(const SGBox<T1>& box, const SGVec3<T2>& v)
89 { return intersects(v, box); }
94 intersects(const SGRay<T>& ray, const SGPlane<T>& plane)
96 // We compute the intersection point
97 // x = origin + \alpha*direction
98 // from the ray origin and non nomalized direction.
99 // For 0 <= \alpha the ray intersects the infinite plane.
100 // The intersection point x can also be written
102 // where n is the planes normal, dist is the distance of the plane from
103 // the origin in normal direction and y is ana aproriate vector
104 // perpendicular to n.
105 // Equate the x values and take the scalar product with the plane normal n.
106 // dot(n, origin) + \alpha*dot(n, direction) = dist
107 // We can now compute alpha from the above equation.
108 // \alpha = (dist - dot(n, origin))/dot(n, direction)
110 // The negative numerator for the \alpha expression
111 T num = plane.getPositiveDist();
112 num -= dot(plane.getNormal(), ray.getOrigin());
114 // If the numerator is zero, we have the rays origin included in the plane
115 if (fabs(num) <= SGLimits<T>::min())
118 // The denominator for the \alpha expression
119 T den = dot(plane.getNormal(), ray.getDirection());
121 // If we get here, we already know that the rays origin is not included
122 // in the plane. Thus if we have a zero denominator we have
123 // a ray paralell to the plane. That is no intersection.
124 if (fabs(den) <= SGLimits<T>::min())
127 // We would now compute \alpha = num/den and compare with 0 and 1.
128 // But to avoid that expensive division, check equation multiplied by
130 T alphaDen = copysign(1, den)*num;
139 intersects(const SGPlane<T>& plane, const SGRay<T>& ray)
140 { return intersects(ray, plane); }
144 intersects(SGVec3<T>& dst, const SGRay<T>& ray, const SGPlane<T>& plane)
146 // We compute the intersection point
147 // x = origin + \alpha*direction
148 // from the ray origin and non nomalized direction.
149 // For 0 <= \alpha the ray intersects the infinite plane.
150 // The intersection point x can also be written
152 // where n is the planes normal, dist is the distance of the plane from
153 // the origin in normal direction and y is ana aproriate vector
154 // perpendicular to n.
155 // Equate the x values and take the scalar product with the plane normal n.
156 // dot(n, origin) + \alpha*dot(n, direction) = dist
157 // We can now compute alpha from the above equation.
158 // \alpha = (dist - dot(n, origin))/dot(n, direction)
160 // The negative numerator for the \alpha expression
161 T num = plane.getPositiveDist();
162 num -= dot(plane.getNormal(), ray.getOrigin());
164 // If the numerator is zero, we have the rays origin included in the plane
165 if (fabs(num) <= SGLimits<T>::min()) {
166 dst = ray.getOrigin();
170 // The denominator for the \alpha expression
171 T den = dot(plane.getNormal(), ray.getDirection());
173 // If we get here, we already know that the rays origin is not included
174 // in the plane. Thus if we have a zero denominator we have
175 // a ray paralell to the plane. That is no intersection.
176 if (fabs(den) <= SGLimits<T>::min())
179 // We would now compute \alpha = num/den and compare with 0 and 1.
180 // But to avoid that expensive division, check equation multiplied by
186 dst = ray.getOrigin() + alpha*ray.getDirection();
192 intersects(SGVec3<T>& dst, const SGPlane<T>& plane, const SGRay<T>& ray)
193 { return intersects(dst, ray, plane); }
197 intersects(const SGLineSegment<T>& lineSegment, const SGPlane<T>& plane)
199 // We compute the intersection point
200 // x = origin + \alpha*direction
201 // from the line segments origin and non nomalized direction.
202 // For 0 <= \alpha <= 1 the line segment intersects the infinite plane.
203 // The intersection point x can also be written
205 // where n is the planes normal, dist is the distance of the plane from
206 // the origin in normal direction and y is ana aproriate vector
207 // perpendicular to n.
208 // Equate the x values and take the scalar product with the plane normal n.
209 // dot(n, origin) + \alpha*dot(n, direction) = dist
210 // We can now compute alpha from the above equation.
211 // \alpha = (dist - dot(n, origin))/dot(n, direction)
213 // The negative numerator for the \alpha expression
214 T num = plane.getPositiveDist();
215 num -= dot(plane.getNormal(), lineSegment.getOrigin());
217 // If the numerator is zero, we have the lines origin included in the plane
218 if (fabs(num) <= SGLimits<T>::min())
221 // The denominator for the \alpha expression
222 T den = dot(plane.getNormal(), lineSegment.getDirection());
224 // If we get here, we already know that the lines origin is not included
225 // in the plane. Thus if we have a zero denominator we have
226 // a line paralell to the plane. That is no intersection.
227 if (fabs(den) <= SGLimits<T>::min())
230 // We would now compute \alpha = num/den and compare with 0 and 1.
231 // But to avoid that expensive division, compare equations
232 // multiplied by |den|. Note that copysign is usually a compiler intrinsic
233 // that expands in assembler code that not even stalls the cpus pipes.
234 T alphaDen = copysign(1, den)*num;
245 intersects(const SGPlane<T>& plane, const SGLineSegment<T>& lineSegment)
246 { return intersects(lineSegment, plane); }
250 intersects(SGVec3<T>& dst, const SGLineSegment<T>& lineSegment, const SGPlane<T>& plane)
252 // We compute the intersection point
253 // x = origin + \alpha*direction
254 // from the line segments origin and non nomalized direction.
255 // For 0 <= \alpha <= 1 the line segment intersects the infinite plane.
256 // The intersection point x can also be written
258 // where n is the planes normal, dist is the distance of the plane from
259 // the origin in normal direction and y is an aproriate vector
260 // perpendicular to n.
261 // Equate the x values and take the scalar product with the plane normal n.
262 // dot(n, origin) + \alpha*dot(n, direction) = dist
263 // We can now compute alpha from the above equation.
264 // \alpha = (dist - dot(n, origin))/dot(n, direction)
266 // The negative numerator for the \alpha expression
267 T num = plane.getPositiveDist();
268 num -= dot(plane.getNormal(), lineSegment.getOrigin());
270 // If the numerator is zero, we have the lines origin included in the plane
271 if (fabs(num) <= SGLimits<T>::min()) {
272 dst = lineSegment.getOrigin();
276 // The denominator for the \alpha expression
277 T den = dot(plane.getNormal(), lineSegment.getDirection());
279 // If we get here, we already know that the lines origin is not included
280 // in the plane. Thus if we have a zero denominator we have
281 // a line paralell to the plane. That is: no intersection.
282 if (fabs(den) <= SGLimits<T>::min())
285 // We would now compute \alpha = num/den and compare with 0 and 1.
286 // But to avoid that expensive division, check equation multiplied by
287 // the denominator. FIXME: shall we do so? or compute like that?
294 dst = lineSegment.getOrigin() + alpha*lineSegment.getDirection();
300 intersects(SGVec3<T>& dst, const SGPlane<T>& plane, const SGLineSegment<T>& lineSegment)
301 { return intersects(dst, lineSegment, plane); }
304 // Distance of a line segment to a point
307 distSqr(const SGLineSegment<T>& lineSeg, const SGVec3<T>& p)
309 SGVec3<T> ps = p - lineSeg.getStart();
311 T psdotdir = dot(ps, lineSeg.getDirection());
315 SGVec3<T> pe = p - lineSeg.getEnd();
316 if (0 <= dot(pe, lineSeg.getDirection()))
319 return dot(ps, ps) - psdotdir*psdotdir/dot(lineSeg.getDirection(), lineSeg.getDirection());
324 distSqr(const SGVec3<T>& p, const SGLineSegment<T>& lineSeg)
325 { return distSqr(lineSeg, p); }
329 dist(const SGVec3<T>& p, const SGLineSegment<T>& lineSeg)
330 { return sqrt(distSqr(lineSeg, p)); }
333 dist(const SGLineSegment<T>& lineSeg, const SGVec3<T>& p)
334 { return sqrt(distSqr(lineSeg, p)); }
338 intersects(const SGRay<T>& ray, const SGSphere<T>& sphere)
340 // See Tomas Akeniene - Moeller/Eric Haines: Real Time Rendering,
341 // second edition, page 571
342 SGVec3<T> l = sphere.getCenter() - ray.getOrigin();
343 T s = dot(l, ray.getDirection());
346 T r2 = sphere.getRadius2();
347 if (s < 0 && l2 > r2)
350 T d2 = dot(ray.getDirection(), ray.getDirection());
351 // The original test would read
352 // T m2 = l2 - s*s/d2;
355 // but to avoid the expensive division, we multiply by d2
365 intersects(const SGSphere<T>& sphere, const SGRay<T>& ray)
366 { return intersects(ray, sphere); }
370 intersects(const SGLineSegment<T>& lineSegment, const SGSphere<T>& sphere)
372 // See Tomas Akeniene - Moeller/Eric Haines: Real Time Rendering,
373 // second edition, page 571
374 SGVec3<T> l = sphere.getCenter() - lineSegment.getStart();
375 T ld = length(lineSegment.getDirection());
376 T s = dot(l, lineSegment.getDirection())/ld;
379 T r2 = sphere.getRadius2();
380 if (s < 0 && l2 > r2)
397 intersects(const SGSphere<T>& sphere, const SGLineSegment<T>& lineSegment)
398 { return intersects(lineSegment, sphere); }
403 // FIXME do not use that default argument later. Just for development now
404 intersects(SGVec3<T>& x, const SGTriangle<T>& tri, const SGRay<T>& ray, T eps = 0)
406 // See Tomas Akeniene - Moeller/Eric Haines: Real Time Rendering
408 // Method based on the observation that we are looking for a
409 // point x that can be expressed in terms of the triangle points
410 // x = v_0 + u*(v_1 - v_0) + v*(v_2 - v_0)
411 // with 0 <= u, v and u + v <= 1.
412 // OTOH it could be expressed in terms of the ray
414 // Now we can compute u, v and t.
415 SGVec3<T> p = cross(ray.getDirection(), tri.getEdge(1));
417 T denom = dot(p, tri.getEdge(0));
418 T signDenom = copysign(1, denom);
420 SGVec3<T> s = ray.getOrigin() - tri.getBaseVertex();
421 SGVec3<T> q = cross(s, tri.getEdge(0));
423 // t = 1/denom*dot(q, tri.getEdge(1));
424 // To avoid an expensive division we multiply by |denom|
425 T tDenom = signDenom*dot(q, tri.getEdge(1));
428 // For line segment we would test against
431 // with the original t. The multiplied test would read
432 // if (absDenom < tDenom)
435 T absDenom = fabs(denom);
436 T absDenomEps = absDenom*eps;
438 // T u = 1/denom*dot(p, s);
439 T u = signDenom*dot(p, s);
440 if (u < -absDenomEps)
442 // T v = 1/denom*dot(q, d);
445 T v = signDenom*dot(q, ray.getDirection());
446 if (v < -absDenomEps)
449 if (u + v > absDenom + absDenomEps)
452 // return if paralell ??? FIXME what if paralell and in plane?
453 // may be we are ok below than anyway??
454 if (absDenom <= SGLimits<T>::min())
458 // if we have survived here it could only happen with denom == 0
459 // that the point is already in plane. Then return the origin ...
460 if (SGLimitsd::min() < absDenom)
461 x += (tDenom/absDenom)*ray.getDirection();
468 intersects(const SGTriangle<T>& tri, const SGRay<T>& ray, T eps = 0)
470 // FIXME: for now just wrap the other method. When that has prooven
471 // well optimized, implement that special case
473 return intersects(dummy, tri, ray, eps);
478 // FIXME do not use that default argument later. Just for development now
479 intersects(SGVec3<T>& x, const SGTriangle<T>& tri, const SGLineSegment<T>& lineSegment, T eps = 0)
481 // See Tomas Akeniene - Moeller/Eric Haines: Real Time Rendering
483 // Method based on the observation that we are looking for a
484 // point x that can be expressed in terms of the triangle points
485 // x = v_0 + u*(v_1 - v_0) + v*(v_2 - v_0)
486 // with 0 <= u, v and u + v <= 1.
487 // OTOH it could be expressed in terms of the lineSegment
489 // Now we can compute u, v and t.
490 SGVec3<T> p = cross(lineSegment.getDirection(), tri.getEdge(1));
492 T denom = dot(p, tri.getEdge(0));
493 T signDenom = copysign(1, denom);
495 SGVec3<T> s = lineSegment.getStart() - tri.getBaseVertex();
496 SGVec3<T> q = cross(s, tri.getEdge(0));
498 // t = 1/denom*dot(q, tri.getEdge(1));
499 // To avoid an expensive division we multiply by |denom|
500 T tDenom = signDenom*dot(q, tri.getEdge(1));
503 // For line segment we would test against
506 // with the original t. The multiplied test reads
507 T absDenom = fabs(denom);
508 if (absDenom < tDenom)
511 // take the CPU accuracy in account
512 T absDenomEps = absDenom*eps;
514 // T u = 1/denom*dot(p, s);
515 T u = signDenom*dot(p, s);
516 if (u < -absDenomEps)
518 // T v = 1/denom*dot(q, d);
521 T v = signDenom*dot(q, lineSegment.getDirection());
522 if (v < -absDenomEps)
525 if (u + v > absDenom + absDenomEps)
528 // return if paralell ??? FIXME what if paralell and in plane?
529 // may be we are ok below than anyway??
530 if (absDenom <= SGLimits<T>::min())
533 x = lineSegment.getStart();
534 // if we have survived here it could only happen with denom == 0
535 // that the point is already in plane. Then return the origin ...
536 if (SGLimitsd::min() < absDenom)
537 x += (tDenom/absDenom)*lineSegment.getDirection();
544 intersects(const SGTriangle<T>& tri, const SGLineSegment<T>& lineSegment, T eps = 0)
546 // FIXME: for now just wrap the other method. When that has prooven
547 // well optimized, implement that special case
549 return intersects(dummy, tri, lineSegment, eps);
555 closestPoint(const SGTriangle<T>& tri, const SGVec3<T>& p)
557 // This method minimizes the distance function Q(u, v) = || p - x ||
558 // where x is a point in the trialgle given by the vertices v_i
559 // x = v_0 + u*(v_1 - v_0) + v*(v_2 - v_0)
560 // The theoretical analysis is somehow too long for a comment.
561 // May be it is sufficient to see that this code passes all the tests.
563 SGVec3<T> off = tri.getBaseVertex() - p;
564 T a = dot(tri.getEdge(0), tri.getEdge(0));
565 T b = dot(tri.getEdge(0), tri.getEdge(1));
566 T c = dot(tri.getEdge(1), tri.getEdge(1));
567 T d = dot(tri.getEdge(0), off);
568 T e = dot(tri.getEdge(1), off);
593 return tri.getBaseVertex() + tri.getEdge(0);
597 return tri.getBaseVertex() + u*tri.getEdge(0);
603 return tri.getBaseVertex();
604 } else if (c <= -e) {
607 return tri.getBaseVertex() + tri.getEdge(1);
611 return tri.getBaseVertex() + v*tri.getEdge(1);
619 return tri.getBaseVertex();
620 } else if (c <= -e) {
623 return tri.getBaseVertex() + tri.getEdge(1);
627 return tri.getBaseVertex() + v*tri.getEdge(1);
635 return tri.getBaseVertex();
636 } else if (a <= -d) {
639 return tri.getBaseVertex() + tri.getEdge(0);
643 return tri.getBaseVertex() + u*tri.getEdge(0);
647 if (det <= SGLimits<T>::min()) {
650 return tri.getBaseVertex();
655 return tri.getBaseVertex() + u*tri.getEdge(0) + v*tri.getEdge(1);
664 T numer = tmp1 - tmp0;
665 T denom = a - 2*b + c;
666 if (denom <= numer) {
669 return tri.getBaseVertex() + tri.getEdge(0);
673 return tri.getBaseVertex() + u*tri.getEdge(0) + v*tri.getEdge(1);
679 return tri.getBaseVertex() + tri.getEdge(1);
683 return tri.getBaseVertex();
687 return tri.getBaseVertex() + v*tri.getEdge(1);
695 T numer = tmp1 - tmp0;
696 T denom = a - 2*b + c;
697 if (denom <= numer) {
700 return tri.getBaseVertex() + tri.getEdge(1);
704 return tri.getBaseVertex() + u*tri.getEdge(0) + v*tri.getEdge(1);
710 return tri.getBaseVertex() + tri.getEdge(0);
714 return tri.getBaseVertex();
718 return tri.getBaseVertex() + u*tri.getEdge(0);
723 T numer = c + e - b - d;
727 return tri.getVertex(2);
729 T denom = a - 2*b + c;
730 if (denom <= numer) {
733 return tri.getBaseVertex() + tri.getEdge(0);
737 return tri.getBaseVertex() + u*tri.getEdge(0) + v*tri.getEdge(1);
745 closestPoint(const SGVec3<T>& p, const SGTriangle<T>& tri)
746 { return closestPoint(tri, p); }
748 template<typename T, typename T2>
750 intersects(const SGTriangle<T>& tri, const SGSphere<T2>& sphere)
752 // This method minimizes the distance function Q(u, v) = || p - x ||
753 // where x is a point in the trialgle given by the vertices v_i
754 // x = v_0 + u*(v_1 - v_0) + v*(v_2 - v_0)
755 // The theoretical analysis is somehow too long for a comment.
756 // May be it is sufficient to see that this code passes all the tests.
758 SGVec3<T> off = tri.getBaseVertex() - SGVec3<T>(sphere.getCenter());
759 T baseDist2 = dot(off, off);
760 T a = dot(tri.getEdge(0), tri.getEdge(0));
761 T b = dot(tri.getEdge(0), tri.getEdge(1));
762 T c = dot(tri.getEdge(1), tri.getEdge(1));
763 T d = dot(tri.getEdge(0), off);
764 T e = dot(tri.getEdge(1), off);
789 T sqrDist = a + 2*d + baseDist2;
790 return sqrDist <= sphere.getRadius2();
794 T sqrDist = d*u + baseDist2;
795 return sqrDist <= sphere.getRadius2();
801 return baseDist2 <= sphere.getRadius2();
802 } else if (c <= -e) {
805 T sqrDist = c + 2*e + baseDist2;
806 return sqrDist <= sphere.getRadius2();
810 T sqrDist = e*v + baseDist2;
811 return sqrDist <= sphere.getRadius2();
819 return baseDist2 <= sphere.getRadius2();
820 } else if (c <= -e) {
823 T sqrDist = c + 2*e + baseDist2;
824 return sqrDist <= sphere.getRadius2();
828 T sqrDist = e*v + baseDist2;
829 return sqrDist <= sphere.getRadius2();
837 return baseDist2 <= sphere.getRadius2();
838 } else if (a <= -d) {
841 T sqrDist = a + 2*d + baseDist2;
842 return sqrDist <= sphere.getRadius2();
846 T sqrDist = d*u + baseDist2;
847 return sqrDist <= sphere.getRadius2();
851 if (det <= SGLimits<T>::min()) {
852 // sqrDist = baseDist2;
855 return baseDist2 <= sphere.getRadius2();
860 T sqrDist = u*(a*u + b*v + 2*d) + v*(b*u + c*v + 2*e) + baseDist2;
861 return sqrDist <= sphere.getRadius2();
870 T numer = tmp1 - tmp0;
871 T denom = a - 2*b + c;
872 if (denom <= numer) {
875 T sqrDist = a + 2*d + baseDist2;
876 return sqrDist <= sphere.getRadius2();
880 T sqrDist = u*(a*u + b*v + 2*d) + v*(b*u + c*v + 2*e) + baseDist2;
881 return sqrDist <= sphere.getRadius2();
887 T sqrDist = c + 2*e + baseDist2;
888 return sqrDist <= sphere.getRadius2();
892 return baseDist2 <= sphere.getRadius2();
896 T sqrDist = e*v + baseDist2;
897 return sqrDist <= sphere.getRadius2();
905 T numer = tmp1 - tmp0;
906 T denom = a - 2*b + c;
907 if (denom <= numer) {
910 T sqrDist = c + 2*e + baseDist2;
911 return sqrDist <= sphere.getRadius2();
915 T sqrDist = u*(a*u + b*v + 2*d) + v*(b*u + c*v + 2*e)+baseDist2;
916 return sqrDist <= sphere.getRadius2();
922 T sqrDist = a + 2*d + baseDist2;
923 return sqrDist <= sphere.getRadius2();
925 // sqrDist = baseDist2;
928 return baseDist2 <= sphere.getRadius2();
932 T sqrDist = d*u + baseDist2;
933 return sqrDist <= sphere.getRadius2();
938 T numer = c + e - b - d;
942 T sqrDist = c + 2*e + baseDist2;
943 return sqrDist <= sphere.getRadius2();
945 T denom = a - 2*b + c;
946 if (denom <= numer) {
949 T sqrDist = a + 2*d + baseDist2;
950 return sqrDist <= sphere.getRadius2();
954 T sqrDist = u*(a*u + b*v + 2*d) + v*(b*u + c*v + 2*e) + baseDist2;
955 return sqrDist <= sphere.getRadius2();
961 template<typename T1, typename T2>
963 intersects(const SGSphere<T1>& sphere, const SGTriangle<T2>& tri)
964 { return intersects(tri, sphere); }
969 intersects(const SGVec3<T>& v, const SGSphere<T>& sphere)
973 return distSqr(v, sphere.getCenter()) <= sphere.getRadius2();
977 intersects(const SGSphere<T>& sphere, const SGVec3<T>& v)
978 { return intersects(v, sphere); }
983 intersects(const SGBox<T>& box, const SGLineSegment<T>& lineSegment)
985 // See Tomas Akeniene - Moeller/Eric Haines: Real Time Rendering
987 SGVec3<T> c = lineSegment.getCenter() - box.getCenter();
988 SGVec3<T> w = T(0.5)*lineSegment.getDirection();
989 SGVec3<T> v(fabs(w.x()), fabs(w.y()), fabs(w.z()));
990 SGVec3<T> h = T(0.5)*box.getSize();
992 if (fabs(c[0]) > v[0] + h[0])
994 if (fabs(c[1]) > v[1] + h[1])
996 if (fabs(c[2]) > v[2] + h[2])
999 if (fabs(c[1]*w[2] - c[2]*w[1]) > h[1]*v[2] + h[2]*v[1])
1001 if (fabs(c[0]*w[2] - c[2]*w[0]) > h[0]*v[2] + h[2]*v[0])
1003 if (fabs(c[0]*w[1] - c[1]*w[0]) > h[0]*v[1] + h[1]*v[0])
1008 template<typename T>
1010 intersects(const SGLineSegment<T>& lineSegment, const SGBox<T>& box)
1011 { return intersects(box, lineSegment); }
1013 template<typename T>
1015 intersects(const SGBox<T>& box, const SGRay<T>& ray)
1017 // See Tomas Akeniene - Moeller/Eric Haines: Real Time Rendering
1019 for (unsigned i = 0; i < 3; ++i) {
1020 T cMin = box.getMin()[i];
1021 T cMax = box.getMax()[i];
1023 T cOrigin = ray.getOrigin()[i];
1025 T cDir = ray.getDirection()[i];
1026 if (fabs(cDir) <= SGLimits<T>::min()) {
1033 T nearr = - SGLimits<T>::max();
1034 T farr = SGLimits<T>::max();
1036 T T1 = (cMin - cOrigin) / cDir;
1037 T T2 = (cMax - cOrigin) / cDir;
1038 if (T1 > T2) std::swap (T1, T2);/* since T1 intersection with near plane */
1039 if (T1 > nearr) nearr = T1; /* want largest Tnear */
1040 if (T2 < farr) farr = T2; /* want smallest Tfarr */
1041 if (nearr > farr) // farr box is missed
1043 if (farr < 0) // box is behind ray
1049 // make it symmetric
1050 template<typename T>
1052 intersects(const SGRay<T>& ray, const SGBox<T>& box)
1053 { return intersects(box, ray); }
1055 template<typename T1, typename T2>
1057 intersects(const SGBox<T1>& box1, const SGBox<T2>& box2)
1059 if (box2.getMax()[0] < box1.getMin()[0])
1061 if (box1.getMax()[0] < box2.getMin()[0])
1064 if (box2.getMax()[1] < box1.getMin()[1])
1066 if (box1.getMax()[1] < box2.getMin()[1])
1069 if (box2.getMax()[2] < box1.getMin()[2])
1071 if (box1.getMax()[2] < box2.getMin()[2])