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)
44 SGVec3<T1> closest = box.getClosestPoint(sphere.getCenter());
45 return distSqr(closest, SGVec3<T1>(sphere.getCenter())) <= sphere.getRadius2();
48 template<typename T1, typename T2>
50 intersects(const SGSphere<T1>& sphere, const SGBox<T2>& box)
51 { return intersects(box, sphere); }
54 template<typename T1, typename T2>
56 intersects(const SGVec3<T1>& v, const SGBox<T2>& box)
58 if (v[0] < box.getMin()[0])
60 if (box.getMax()[0] < v[0])
62 if (v[1] < box.getMin()[1])
64 if (box.getMax()[1] < v[1])
66 if (v[2] < box.getMin()[2])
68 if (box.getMax()[2] < v[2])
72 template<typename T1, typename T2>
74 intersects(const SGBox<T1>& box, const SGVec3<T2>& v)
75 { return intersects(v, box); }
80 intersects(const SGRay<T>& ray, const SGPlane<T>& plane)
82 // We compute the intersection point
83 // x = origin + \alpha*direction
84 // from the ray origin and non nomalized direction.
85 // For 0 <= \alpha the ray intersects the infinite plane.
86 // The intersection point x can also be written
88 // where n is the planes normal, dist is the distance of the plane from
89 // the origin in normal direction and y is ana aproriate vector
90 // perpendicular to n.
91 // Equate the x values and take the scalar product with the plane normal n.
92 // dot(n, origin) + \alpha*dot(n, direction) = dist
93 // We can now compute alpha from the above equation.
94 // \alpha = (dist - dot(n, origin))/dot(n, direction)
96 // The negative numerator for the \alpha expression
97 T num = plane.getPositiveDist();
98 num -= dot(plane.getNormal(), ray.getOrigin());
100 // If the numerator is zero, we have the rays origin included in the plane
101 if (fabs(num) <= SGLimits<T>::min())
104 // The denominator for the \alpha expression
105 T den = dot(plane.getNormal(), ray.getDirection());
107 // If we get here, we already know that the rays origin is not included
108 // in the plane. Thus if we have a zero denominator we have
109 // a ray paralell to the plane. That is no intersection.
110 if (fabs(den) <= SGLimits<T>::min())
113 // We would now compute \alpha = num/den and compare with 0 and 1.
114 // But to avoid that expensive division, check equation multiplied by
116 T alphaDen = copysign(1, den)*num;
125 intersects(const SGPlane<T>& plane, const SGRay<T>& ray)
126 { return intersects(ray, plane); }
130 intersects(SGVec3<T>& dst, const SGRay<T>& ray, const SGPlane<T>& plane)
132 // We compute the intersection point
133 // x = origin + \alpha*direction
134 // from the ray origin and non nomalized direction.
135 // For 0 <= \alpha the ray intersects the infinite plane.
136 // The intersection point x can also be written
138 // where n is the planes normal, dist is the distance of the plane from
139 // the origin in normal direction and y is ana aproriate vector
140 // perpendicular to n.
141 // Equate the x values and take the scalar product with the plane normal n.
142 // dot(n, origin) + \alpha*dot(n, direction) = dist
143 // We can now compute alpha from the above equation.
144 // \alpha = (dist - dot(n, origin))/dot(n, direction)
146 // The negative numerator for the \alpha expression
147 T num = plane.getPositiveDist();
148 num -= dot(plane.getNormal(), ray.getOrigin());
150 // If the numerator is zero, we have the rays origin included in the plane
151 if (fabs(num) <= SGLimits<T>::min()) {
152 dst = ray.getOrigin();
156 // The denominator for the \alpha expression
157 T den = dot(plane.getNormal(), ray.getDirection());
159 // If we get here, we already know that the rays origin is not included
160 // in the plane. Thus if we have a zero denominator we have
161 // a ray paralell to the plane. That is no intersection.
162 if (fabs(den) <= SGLimits<T>::min())
165 // We would now compute \alpha = num/den and compare with 0 and 1.
166 // But to avoid that expensive division, check equation multiplied by
172 dst = ray.getOrigin() + alpha*ray.getDirection();
178 intersects(SGVec3<T>& dst, const SGPlane<T>& plane, const SGRay<T>& ray)
179 { return intersects(dst, ray, plane); }
183 intersects(const SGLineSegment<T>& lineSegment, const SGPlane<T>& plane)
185 // We compute the intersection point
186 // x = origin + \alpha*direction
187 // from the line segments origin and non nomalized direction.
188 // For 0 <= \alpha <= 1 the line segment intersects the infinite plane.
189 // The intersection point x can also be written
191 // where n is the planes normal, dist is the distance of the plane from
192 // the origin in normal direction and y is ana aproriate vector
193 // perpendicular to n.
194 // Equate the x values and take the scalar product with the plane normal n.
195 // dot(n, origin) + \alpha*dot(n, direction) = dist
196 // We can now compute alpha from the above equation.
197 // \alpha = (dist - dot(n, origin))/dot(n, direction)
199 // The negative numerator for the \alpha expression
200 T num = plane.getPositiveDist();
201 num -= dot(plane.getNormal(), lineSegment.getOrigin());
203 // If the numerator is zero, we have the lines origin included in the plane
204 if (fabs(num) <= SGLimits<T>::min())
207 // The denominator for the \alpha expression
208 T den = dot(plane.getNormal(), lineSegment.getDirection());
210 // If we get here, we already know that the lines origin is not included
211 // in the plane. Thus if we have a zero denominator we have
212 // a line paralell to the plane. That is no intersection.
213 if (fabs(den) <= SGLimits<T>::min())
216 // We would now compute \alpha = num/den and compare with 0 and 1.
217 // But to avoid that expensive division, compare equations
218 // multiplied by |den|. Note that copysign is usually a compiler intrinsic
219 // that expands in assembler code that not even stalls the cpus pipes.
220 T alphaDen = copysign(1, den)*num;
231 intersects(const SGPlane<T>& plane, const SGLineSegment<T>& lineSegment)
232 { return intersects(lineSegment, plane); }
236 intersects(SGVec3<T>& dst, const SGLineSegment<T>& lineSegment, const SGPlane<T>& plane)
238 // We compute the intersection point
239 // x = origin + \alpha*direction
240 // from the line segments origin and non nomalized direction.
241 // For 0 <= \alpha <= 1 the line segment intersects the infinite plane.
242 // The intersection point x can also be written
244 // where n is the planes normal, dist is the distance of the plane from
245 // the origin in normal direction and y is an aproriate vector
246 // perpendicular to n.
247 // Equate the x values and take the scalar product with the plane normal n.
248 // dot(n, origin) + \alpha*dot(n, direction) = dist
249 // We can now compute alpha from the above equation.
250 // \alpha = (dist - dot(n, origin))/dot(n, direction)
252 // The negative numerator for the \alpha expression
253 T num = plane.getPositiveDist();
254 num -= dot(plane.getNormal(), lineSegment.getOrigin());
256 // If the numerator is zero, we have the lines origin included in the plane
257 if (fabs(num) <= SGLimits<T>::min()) {
258 dst = lineSegment.getOrigin();
262 // The denominator for the \alpha expression
263 T den = dot(plane.getNormal(), lineSegment.getDirection());
265 // If we get here, we already know that the lines origin is not included
266 // in the plane. Thus if we have a zero denominator we have
267 // a line paralell to the plane. That is: no intersection.
268 if (fabs(den) <= SGLimits<T>::min())
271 // We would now compute \alpha = num/den and compare with 0 and 1.
272 // But to avoid that expensive division, check equation multiplied by
273 // the denominator. FIXME: shall we do so? or compute like that?
280 dst = lineSegment.getOrigin() + alpha*lineSegment.getDirection();
286 intersects(SGVec3<T>& dst, const SGPlane<T>& plane, const SGLineSegment<T>& lineSegment)
287 { return intersects(dst, lineSegment, plane); }
290 // Distance of a line segment to a point
293 distSqr(const SGLineSegment<T>& lineSeg, const SGVec3<T>& p)
295 SGVec3<T> ps = p - lineSeg.getStart();
297 T psdotdir = dot(ps, lineSeg.getDirection());
301 SGVec3<T> pe = p - lineSeg.getEnd();
302 if (0 <= dot(pe, lineSeg.getDirection()))
305 return dot(ps, ps) - psdotdir*psdotdir/dot(lineSeg.getDirection(), lineSeg.getDirection());
310 distSqr(const SGVec3<T>& p, const SGLineSegment<T>& lineSeg)
311 { return distSqr(lineSeg, p); }
315 dist(const SGVec3<T>& p, const SGLineSegment<T>& lineSeg)
316 { return sqrt(distSqr(lineSeg, p)); }
319 dist(const SGLineSegment<T>& lineSeg, const SGVec3<T>& p)
320 { return sqrt(distSqr(lineSeg, p)); }
324 intersects(const SGRay<T>& ray, const SGSphere<T>& sphere)
326 // See Tomas Akeniene - Moeller/Eric Haines: Real Time Rendering,
327 // second edition, page 571
328 SGVec3<T> l = sphere.getCenter() - ray.getOrigin();
329 T s = dot(l, ray.getDirection());
332 T r2 = sphere.getRadius2();
333 if (s < 0 && l2 > r2)
336 T d2 = dot(ray.getDirection(), ray.getDirection());
337 // The original test would read
338 // T m2 = l2 - s*s/d2;
341 // but to avoid the expensive division, we multiply by d2
351 intersects(const SGSphere<T>& sphere, const SGRay<T>& ray)
352 { return intersects(ray, sphere); }
356 intersects(const SGLineSegment<T>& lineSegment, const SGSphere<T>& sphere)
358 // See Tomas Akeniene - Moeller/Eric Haines: Real Time Rendering,
359 // second edition, page 571
360 SGVec3<T> l = sphere.getCenter() - lineSegment.getStart();
361 T ld = length(lineSegment.getDirection());
362 T s = dot(l, lineSegment.getDirection())/ld;
365 T r2 = sphere.getRadius2();
366 if (s < 0 && l2 > r2)
383 intersects(const SGSphere<T>& sphere, const SGLineSegment<T>& lineSegment)
384 { return intersects(lineSegment, sphere); }
389 // FIXME do not use that default argument later. Just for development now
390 intersects(SGVec3<T>& x, const SGTriangle<T>& tri, const SGRay<T>& ray, T eps = 0)
392 // See Tomas Akeniene - Moeller/Eric Haines: Real Time Rendering
394 // Method based on the observation that we are looking for a
395 // point x that can be expressed in terms of the triangle points
396 // x = v_0 + u*(v_1 - v_0) + v*(v_2 - v_0)
397 // with 0 <= u, v and u + v <= 1.
398 // OTOH it could be expressed in terms of the ray
400 // Now we can compute u, v and t.
401 SGVec3<T> p = cross(ray.getDirection(), tri.getEdge(1));
403 T denom = dot(p, tri.getEdge(0));
404 T signDenom = copysign(1, denom);
406 SGVec3<T> s = ray.getOrigin() - tri.getBaseVertex();
407 SGVec3<T> q = cross(s, tri.getEdge(0));
409 // t = 1/denom*dot(q, tri.getEdge(1));
410 // To avoid an expensive division we multiply by |denom|
411 T tDenom = signDenom*dot(q, tri.getEdge(1));
414 // For line segment we would test against
417 // with the original t. The multiplied test would read
418 // if (absDenom < tDenom)
421 T absDenom = fabs(denom);
422 T absDenomEps = absDenom*eps;
424 // T u = 1/denom*dot(p, s);
425 T u = signDenom*dot(p, s);
426 if (u < -absDenomEps)
428 // T v = 1/denom*dot(q, d);
431 T v = signDenom*dot(q, ray.getDirection());
432 if (v < -absDenomEps)
435 if (u + v > absDenom + absDenomEps)
438 // return if paralell ??? FIXME what if paralell and in plane?
439 // may be we are ok below than anyway??
440 if (absDenom <= SGLimits<T>::min())
444 // if we have survived here it could only happen with denom == 0
445 // that the point is already in plane. Then return the origin ...
446 if (SGLimitsd::min() < absDenom)
447 x += (tDenom/absDenom)*ray.getDirection();
454 intersects(const SGTriangle<T>& tri, const SGRay<T>& ray, T eps = 0)
456 // FIXME: for now just wrap the other method. When that has prooven
457 // well optimized, implement that special case
459 return intersects(dummy, tri, ray, eps);
464 // FIXME do not use that default argument later. Just for development now
465 intersects(SGVec3<T>& x, const SGTriangle<T>& tri, const SGLineSegment<T>& lineSegment, T eps = 0)
467 // See Tomas Akeniene - Moeller/Eric Haines: Real Time Rendering
469 // Method based on the observation that we are looking for a
470 // point x that can be expressed in terms of the triangle points
471 // x = v_0 + u*(v_1 - v_0) + v*(v_2 - v_0)
472 // with 0 <= u, v and u + v <= 1.
473 // OTOH it could be expressed in terms of the lineSegment
475 // Now we can compute u, v and t.
476 SGVec3<T> p = cross(lineSegment.getDirection(), tri.getEdge(1));
478 T denom = dot(p, tri.getEdge(0));
479 T signDenom = copysign(1, denom);
481 SGVec3<T> s = lineSegment.getStart() - tri.getBaseVertex();
482 SGVec3<T> q = cross(s, tri.getEdge(0));
484 // t = 1/denom*dot(q, tri.getEdge(1));
485 // To avoid an expensive division we multiply by |denom|
486 T tDenom = signDenom*dot(q, tri.getEdge(1));
489 // For line segment we would test against
492 // with the original t. The multiplied test reads
493 T absDenom = fabs(denom);
494 if (absDenom < tDenom)
497 // take the CPU accuracy in account
498 T absDenomEps = absDenom*eps;
500 // T u = 1/denom*dot(p, s);
501 T u = signDenom*dot(p, s);
502 if (u < -absDenomEps)
504 // T v = 1/denom*dot(q, d);
507 T v = signDenom*dot(q, lineSegment.getDirection());
508 if (v < -absDenomEps)
511 if (u + v > absDenom + absDenomEps)
514 // return if paralell ??? FIXME what if paralell and in plane?
515 // may be we are ok below than anyway??
516 if (absDenom <= SGLimits<T>::min())
519 x = lineSegment.getStart();
520 // if we have survived here it could only happen with denom == 0
521 // that the point is already in plane. Then return the origin ...
522 if (SGLimitsd::min() < absDenom)
523 x += (tDenom/absDenom)*lineSegment.getDirection();
530 intersects(const SGTriangle<T>& tri, const SGLineSegment<T>& lineSegment, T eps = 0)
532 // FIXME: for now just wrap the other method. When that has prooven
533 // well optimized, implement that special case
535 return intersects(dummy, tri, lineSegment, eps);
541 closestPoint(const SGTriangle<T>& tri, const SGVec3<T>& p)
543 // This method minimizes the distance function Q(u, v) = || p - x ||
544 // where x is a point in the trialgle given by the vertices v_i
545 // x = v_0 + u*(v_1 - v_0) + v*(v_2 - v_0)
546 // The theoretical analysis is somehow too long for a comment.
547 // May be it is sufficient to see that this code passes all the tests.
549 SGVec3<T> off = tri.getBaseVertex() - p;
550 T a = dot(tri.getEdge(0), tri.getEdge(0));
551 T b = dot(tri.getEdge(0), tri.getEdge(1));
552 T c = dot(tri.getEdge(1), tri.getEdge(1));
553 T d = dot(tri.getEdge(0), off);
554 T e = dot(tri.getEdge(1), off);
579 return tri.getBaseVertex() + tri.getEdge(0);
583 return tri.getBaseVertex() + u*tri.getEdge(0);
589 return tri.getBaseVertex();
590 } else if (c <= -e) {
593 return tri.getBaseVertex() + tri.getEdge(1);
597 return tri.getBaseVertex() + v*tri.getEdge(1);
605 return tri.getBaseVertex();
606 } else if (c <= -e) {
609 return tri.getBaseVertex() + tri.getEdge(1);
613 return tri.getBaseVertex() + v*tri.getEdge(1);
621 return tri.getBaseVertex();
622 } else if (a <= -d) {
625 return tri.getBaseVertex() + tri.getEdge(0);
629 return tri.getBaseVertex() + u*tri.getEdge(0);
633 if (det <= SGLimits<T>::min()) {
636 return tri.getBaseVertex();
641 return tri.getBaseVertex() + u*tri.getEdge(0) + v*tri.getEdge(1);
650 T numer = tmp1 - tmp0;
651 T denom = a - 2*b + c;
652 if (denom <= numer) {
655 return tri.getBaseVertex() + tri.getEdge(0);
659 return tri.getBaseVertex() + u*tri.getEdge(0) + v*tri.getEdge(1);
665 return tri.getBaseVertex() + tri.getEdge(1);
669 return tri.getBaseVertex();
673 return tri.getBaseVertex() + v*tri.getEdge(1);
681 T numer = tmp1 - tmp0;
682 T denom = a - 2*b + c;
683 if (denom <= numer) {
686 return tri.getBaseVertex() + tri.getEdge(1);
690 return tri.getBaseVertex() + u*tri.getEdge(0) + v*tri.getEdge(1);
696 return tri.getBaseVertex() + tri.getEdge(0);
700 return tri.getBaseVertex();
704 return tri.getBaseVertex() + u*tri.getEdge(0);
709 T numer = c + e - b - d;
713 return tri.getVertex(2);
715 T denom = a - 2*b + c;
716 if (denom <= numer) {
719 return tri.getBaseVertex() + tri.getEdge(0);
723 return tri.getBaseVertex() + u*tri.getEdge(0) + v*tri.getEdge(1);
731 closestPoint(const SGVec3<T>& p, const SGTriangle<T>& tri)
732 { return closestPoint(tri, p); }
734 template<typename T, typename T2>
736 intersects(const SGTriangle<T>& tri, const SGSphere<T2>& sphere)
738 // This method minimizes the distance function Q(u, v) = || p - x ||
739 // where x is a point in the trialgle given by the vertices v_i
740 // x = v_0 + u*(v_1 - v_0) + v*(v_2 - v_0)
741 // The theoretical analysis is somehow too long for a comment.
742 // May be it is sufficient to see that this code passes all the tests.
744 SGVec3<T> off = tri.getBaseVertex() - SGVec3<T>(sphere.getCenter());
745 T baseDist2 = dot(off, off);
746 T a = dot(tri.getEdge(0), tri.getEdge(0));
747 T b = dot(tri.getEdge(0), tri.getEdge(1));
748 T c = dot(tri.getEdge(1), tri.getEdge(1));
749 T d = dot(tri.getEdge(0), off);
750 T e = dot(tri.getEdge(1), off);
775 T sqrDist = a + 2*d + baseDist2;
776 return sqrDist <= sphere.getRadius2();
780 T sqrDist = d*u + baseDist2;
781 return sqrDist <= sphere.getRadius2();
787 return baseDist2 <= sphere.getRadius2();
788 } else if (c <= -e) {
791 T sqrDist = c + 2*e + baseDist2;
792 return sqrDist <= sphere.getRadius2();
796 T sqrDist = e*v + baseDist2;
797 return sqrDist <= sphere.getRadius2();
805 return baseDist2 <= sphere.getRadius2();
806 } else if (c <= -e) {
809 T sqrDist = c + 2*e + baseDist2;
810 return sqrDist <= sphere.getRadius2();
814 T sqrDist = e*v + baseDist2;
815 return sqrDist <= sphere.getRadius2();
823 return baseDist2 <= sphere.getRadius2();
824 } else if (a <= -d) {
827 T sqrDist = a + 2*d + baseDist2;
828 return sqrDist <= sphere.getRadius2();
832 T sqrDist = d*u + baseDist2;
833 return sqrDist <= sphere.getRadius2();
837 if (det <= SGLimits<T>::min()) {
838 // sqrDist = baseDist2;
841 return baseDist2 <= sphere.getRadius2();
846 T sqrDist = u*(a*u + b*v + 2*d) + v*(b*u + c*v + 2*e) + baseDist2;
847 return sqrDist <= sphere.getRadius2();
856 T numer = tmp1 - tmp0;
857 T denom = a - 2*b + c;
858 if (denom <= numer) {
861 T sqrDist = a + 2*d + baseDist2;
862 return sqrDist <= sphere.getRadius2();
866 T sqrDist = u*(a*u + b*v + 2*d) + v*(b*u + c*v + 2*e) + baseDist2;
867 return sqrDist <= sphere.getRadius2();
873 T sqrDist = c + 2*e + baseDist2;
874 return sqrDist <= sphere.getRadius2();
878 return baseDist2 <= sphere.getRadius2();
882 T sqrDist = e*v + baseDist2;
883 return sqrDist <= sphere.getRadius2();
891 T numer = tmp1 - tmp0;
892 T denom = a - 2*b + c;
893 if (denom <= numer) {
896 T sqrDist = c + 2*e + baseDist2;
897 return sqrDist <= sphere.getRadius2();
901 T sqrDist = u*(a*u + b*v + 2*d) + v*(b*u + c*v + 2*e)+baseDist2;
902 return sqrDist <= sphere.getRadius2();
908 T sqrDist = a + 2*d + baseDist2;
909 return sqrDist <= sphere.getRadius2();
911 // sqrDist = baseDist2;
914 return baseDist2 <= sphere.getRadius2();
918 T sqrDist = d*u + baseDist2;
919 return sqrDist <= sphere.getRadius2();
924 T numer = c + e - b - d;
928 T sqrDist = c + 2*e + baseDist2;
929 return sqrDist <= sphere.getRadius2();
931 T denom = a - 2*b + c;
932 if (denom <= numer) {
935 T sqrDist = a + 2*d + baseDist2;
936 return sqrDist <= sphere.getRadius2();
940 T sqrDist = u*(a*u + b*v + 2*d) + v*(b*u + c*v + 2*e) + baseDist2;
941 return sqrDist <= sphere.getRadius2();
947 template<typename T1, typename T2>
949 intersects(const SGSphere<T1>& sphere, const SGTriangle<T2>& tri)
950 { return intersects(tri, sphere); }
955 intersects(const SGVec3<T>& v, const SGSphere<T>& sphere)
959 return distSqr(v, sphere.getCenter()) <= sphere.getRadius2();
963 intersects(const SGSphere<T>& sphere, const SGVec3<T>& v)
964 { return intersects(v, sphere); }
969 intersects(const SGBox<T>& box, const SGLineSegment<T>& lineSegment)
971 // See Tomas Akeniene - Moeller/Eric Haines: Real Time Rendering
973 SGVec3<T> c = lineSegment.getCenter() - box.getCenter();
974 SGVec3<T> w = T(0.5)*lineSegment.getDirection();
975 SGVec3<T> v(fabs(w.x()), fabs(w.y()), fabs(w.z()));
976 SGVec3<T> h = T(0.5)*box.getSize();
978 if (fabs(c[0]) > v[0] + h[0])
980 if (fabs(c[1]) > v[1] + h[1])
982 if (fabs(c[2]) > v[2] + h[2])
985 if (fabs(c[1]*w[2] - c[2]*w[1]) > h[1]*v[2] + h[2]*v[1])
987 if (fabs(c[0]*w[2] - c[2]*w[0]) > h[0]*v[2] + h[2]*v[0])
989 if (fabs(c[0]*w[1] - c[1]*w[0]) > h[0]*v[1] + h[1]*v[0])
996 intersects(const SGLineSegment<T>& lineSegment, const SGBox<T>& box)
997 { return intersects(box, lineSegment); }
1001 intersects(const SGBox<T>& box, const SGRay<T>& ray)
1003 // See Tomas Akeniene - Moeller/Eric Haines: Real Time Rendering
1005 for (unsigned i = 0; i < 3; ++i) {
1006 T cMin = box.getMin()[i];
1007 T cMax = box.getMax()[i];
1009 T cOrigin = ray.getOrigin()[i];
1011 T cDir = ray.getDirection()[i];
1012 if (fabs(cDir) <= SGLimits<T>::min()) {
1019 T nearr = - SGLimits<T>::max();
1020 T farr = SGLimits<T>::max();
1022 T T1 = (cMin - cOrigin) / cDir;
1023 T T2 = (cMax - cOrigin) / cDir;
1024 if (T1 > T2) std::swap (T1, T2);/* since T1 intersection with near plane */
1025 if (T1 > nearr) nearr = T1; /* want largest Tnear */
1026 if (T2 < farr) farr = T2; /* want smallest Tfarr */
1027 if (nearr > farr) // farr box is missed
1029 if (farr < 0) // box is behind ray
1035 // make it symmetric
1036 template<typename T>
1038 intersects(const SGRay<T>& ray, const SGBox<T>& box)
1039 { return intersects(box, ray); }
1041 template<typename T1, typename T2>
1043 intersects(const SGBox<T1>& box1, const SGBox<T2>& box2)
1045 if (box2.getMax()[0] < box1.getMin()[0])
1047 if (box1.getMax()[0] < box2.getMin()[0])
1050 if (box2.getMax()[1] < box1.getMin()[1])
1052 if (box1.getMax()[1] < box2.getMin()[1])
1055 if (box2.getMax()[2] < box1.getMin()[2])
1057 if (box1.getMax()[2] < box2.getMin()[2])