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
25 intersects(const SGSphere<T>& s1, const SGSphere<T>& s2)
32 T dist = s1.getRadius() + s2.getRadius();
33 return distSqr(s1.getCenter(), s2.getCenter()) <= dist*dist;
37 template<typename T1, typename T2>
39 intersects(const SGBox<T1>& box, const SGSphere<T2>& sphere)
46 SGVec3<T1> closest = box.getClosestPoint(sphere.getCenter());
47 return distSqr(closest, SGVec3<T1>(sphere.getCenter())) <= sphere.getRadius2();
50 template<typename T1, typename T2>
52 intersects(const SGSphere<T1>& sphere, const SGBox<T2>& box)
53 { return intersects(box, sphere); }
56 template<typename T1, typename T2>
58 intersects(const SGVec3<T1>& v, const SGBox<T2>& box)
60 if (v[0] < box.getMin()[0])
62 if (box.getMax()[0] < v[0])
64 if (v[1] < box.getMin()[1])
66 if (box.getMax()[1] < v[1])
68 if (v[2] < box.getMin()[2])
70 if (box.getMax()[2] < v[2])
74 template<typename T1, typename T2>
76 intersects(const SGBox<T1>& box, const SGVec3<T2>& v)
77 { return intersects(v, box); }
82 intersects(const SGRay<T>& ray, const SGPlane<T>& plane)
84 // We compute the intersection point
85 // x = origin + \alpha*direction
86 // from the ray origin and non nomalized direction.
87 // For 0 <= \alpha the ray intersects the infinite plane.
88 // The intersection point x can also be written
90 // where n is the planes normal, dist is the distance of the plane from
91 // the origin in normal direction and y is ana aproriate vector
92 // perpendicular to n.
93 // Equate the x values and take the scalar product with the plane normal n.
94 // dot(n, origin) + \alpha*dot(n, direction) = dist
95 // We can now compute alpha from the above equation.
96 // \alpha = (dist - dot(n, origin))/dot(n, direction)
98 // The negative numerator for the \alpha expression
99 T num = plane.getPositiveDist();
100 num -= dot(plane.getNormal(), ray.getOrigin());
102 // If the numerator is zero, we have the rays origin included in the plane
103 if (fabs(num) <= SGLimits<T>::min())
106 // The denominator for the \alpha expression
107 T den = dot(plane.getNormal(), ray.getDirection());
109 // If we get here, we already know that the rays origin is not included
110 // in the plane. Thus if we have a zero denominator we have
111 // a ray paralell to the plane. That is no intersection.
112 if (fabs(den) <= SGLimits<T>::min())
115 // We would now compute \alpha = num/den and compare with 0 and 1.
116 // But to avoid that expensive division, check equation multiplied by
118 T alphaDen = copysign(1, den)*num;
127 intersects(const SGPlane<T>& plane, const SGRay<T>& ray)
128 { return intersects(ray, plane); }
132 intersects(SGVec3<T>& dst, const SGRay<T>& ray, const SGPlane<T>& plane)
134 // We compute the intersection point
135 // x = origin + \alpha*direction
136 // from the ray origin and non nomalized direction.
137 // For 0 <= \alpha the ray intersects the infinite plane.
138 // The intersection point x can also be written
140 // where n is the planes normal, dist is the distance of the plane from
141 // the origin in normal direction and y is ana aproriate vector
142 // perpendicular to n.
143 // Equate the x values and take the scalar product with the plane normal n.
144 // dot(n, origin) + \alpha*dot(n, direction) = dist
145 // We can now compute alpha from the above equation.
146 // \alpha = (dist - dot(n, origin))/dot(n, direction)
148 // The negative numerator for the \alpha expression
149 T num = plane.getPositiveDist();
150 num -= dot(plane.getNormal(), ray.getOrigin());
152 // If the numerator is zero, we have the rays origin included in the plane
153 if (fabs(num) <= SGLimits<T>::min()) {
154 dst = ray.getOrigin();
158 // The denominator for the \alpha expression
159 T den = dot(plane.getNormal(), ray.getDirection());
161 // If we get here, we already know that the rays origin is not included
162 // in the plane. Thus if we have a zero denominator we have
163 // a ray paralell to the plane. That is no intersection.
164 if (fabs(den) <= SGLimits<T>::min())
167 // We would now compute \alpha = num/den and compare with 0 and 1.
168 // But to avoid that expensive division, check equation multiplied by
174 dst = ray.getOrigin() + alpha*ray.getDirection();
180 intersects(SGVec3<T>& dst, const SGPlane<T>& plane, const SGRay<T>& ray)
181 { return intersects(dst, ray, plane); }
185 intersects(const SGLineSegment<T>& lineSegment, const SGPlane<T>& plane)
187 // We compute the intersection point
188 // x = origin + \alpha*direction
189 // from the line segments origin and non nomalized direction.
190 // For 0 <= \alpha <= 1 the line segment intersects the infinite plane.
191 // The intersection point x can also be written
193 // where n is the planes normal, dist is the distance of the plane from
194 // the origin in normal direction and y is ana aproriate vector
195 // perpendicular to n.
196 // Equate the x values and take the scalar product with the plane normal n.
197 // dot(n, origin) + \alpha*dot(n, direction) = dist
198 // We can now compute alpha from the above equation.
199 // \alpha = (dist - dot(n, origin))/dot(n, direction)
201 // The negative numerator for the \alpha expression
202 T num = plane.getPositiveDist();
203 num -= dot(plane.getNormal(), lineSegment.getOrigin());
205 // If the numerator is zero, we have the lines origin included in the plane
206 if (fabs(num) <= SGLimits<T>::min())
209 // The denominator for the \alpha expression
210 T den = dot(plane.getNormal(), lineSegment.getDirection());
212 // If we get here, we already know that the lines origin is not included
213 // in the plane. Thus if we have a zero denominator we have
214 // a line paralell to the plane. That is no intersection.
215 if (fabs(den) <= SGLimits<T>::min())
218 // We would now compute \alpha = num/den and compare with 0 and 1.
219 // But to avoid that expensive division, compare equations
220 // multiplied by |den|. Note that copysign is usually a compiler intrinsic
221 // that expands in assembler code that not even stalls the cpus pipes.
222 T alphaDen = copysign(1, den)*num;
233 intersects(const SGPlane<T>& plane, const SGLineSegment<T>& lineSegment)
234 { return intersects(lineSegment, plane); }
238 intersects(SGVec3<T>& dst, const SGLineSegment<T>& lineSegment, const SGPlane<T>& plane)
240 // We compute the intersection point
241 // x = origin + \alpha*direction
242 // from the line segments origin and non nomalized direction.
243 // For 0 <= \alpha <= 1 the line segment intersects the infinite plane.
244 // The intersection point x can also be written
246 // where n is the planes normal, dist is the distance of the plane from
247 // the origin in normal direction and y is an aproriate vector
248 // perpendicular to n.
249 // Equate the x values and take the scalar product with the plane normal n.
250 // dot(n, origin) + \alpha*dot(n, direction) = dist
251 // We can now compute alpha from the above equation.
252 // \alpha = (dist - dot(n, origin))/dot(n, direction)
254 // The negative numerator for the \alpha expression
255 T num = plane.getPositiveDist();
256 num -= dot(plane.getNormal(), lineSegment.getStart());
258 // If the numerator is zero, we have the lines origin included in the plane
259 if (fabs(num) <= SGLimits<T>::min()) {
260 dst = lineSegment.getStart();
264 // The denominator for the \alpha expression
265 T den = dot(plane.getNormal(), lineSegment.getDirection());
267 // If we get here, we already know that the lines origin is not included
268 // in the plane. Thus if we have a zero denominator we have
269 // a line paralell to the plane. That is: no intersection.
270 if (fabs(den) <= SGLimits<T>::min())
273 // We would now compute \alpha = num/den and compare with 0 and 1.
274 // But to avoid that expensive division, check equation multiplied by
275 // the denominator. FIXME: shall we do so? or compute like that?
282 dst = lineSegment.getStart() + alpha*lineSegment.getDirection();
288 intersects(SGVec3<T>& dst, const SGPlane<T>& plane, const SGLineSegment<T>& lineSegment)
289 { return intersects(dst, lineSegment, plane); }
292 // Distance of a line segment to a point
295 distSqr(const SGLineSegment<T>& lineSeg, const SGVec3<T>& p)
297 SGVec3<T> ps = p - lineSeg.getStart();
299 T psdotdir = dot(ps, lineSeg.getDirection());
303 SGVec3<T> pe = p - lineSeg.getEnd();
304 if (0 <= dot(pe, lineSeg.getDirection()))
307 return dot(ps, ps) - psdotdir*psdotdir/dot(lineSeg.getDirection(), lineSeg.getDirection());
312 distSqr(const SGVec3<T>& p, const SGLineSegment<T>& lineSeg)
313 { return distSqr(lineSeg, p); }
317 dist(const SGVec3<T>& p, const SGLineSegment<T>& lineSeg)
318 { return sqrt(distSqr(lineSeg, p)); }
321 dist(const SGLineSegment<T>& lineSeg, const SGVec3<T>& p)
322 { return sqrt(distSqr(lineSeg, p)); }
326 intersects(const SGRay<T>& ray, const SGSphere<T>& sphere)
328 // See Tomas Akeniene - Moeller/Eric Haines: Real Time Rendering,
329 // second edition, page 571
330 SGVec3<T> l = sphere.getCenter() - ray.getOrigin();
331 T s = dot(l, ray.getDirection());
334 T r2 = sphere.getRadius2();
335 if (s < 0 && l2 > r2)
338 T d2 = dot(ray.getDirection(), ray.getDirection());
339 // The original test would read
340 // T m2 = l2 - s*s/d2;
343 // but to avoid the expensive division, we multiply by d2
353 intersects(const SGSphere<T>& sphere, const SGRay<T>& ray)
354 { return intersects(ray, sphere); }
358 intersects(const SGLineSegment<T>& lineSegment, const SGSphere<T>& sphere)
360 // See Tomas Akeniene - Moeller/Eric Haines: Real Time Rendering,
361 // second edition, page 571
362 SGVec3<T> l = sphere.getCenter() - lineSegment.getStart();
363 T ld = length(lineSegment.getDirection());
364 T s = dot(l, lineSegment.getDirection())/ld;
367 T r2 = sphere.getRadius2();
368 if (s < 0 && l2 > r2)
385 intersects(const SGSphere<T>& sphere, const SGLineSegment<T>& lineSegment)
386 { return intersects(lineSegment, sphere); }
391 // FIXME do not use that default argument later. Just for development now
392 intersects(SGVec3<T>& x, const SGTriangle<T>& tri, const SGRay<T>& ray, T eps = 0)
394 // See Tomas Akeniene - Moeller/Eric Haines: Real Time Rendering
396 // Method based on the observation that we are looking for a
397 // point x that can be expressed in terms of the triangle points
398 // x = v_0 + u*(v_1 - v_0) + v*(v_2 - v_0)
399 // with 0 <= u, v and u + v <= 1.
400 // OTOH it could be expressed in terms of the ray
402 // Now we can compute u, v and t.
403 SGVec3<T> p = cross(ray.getDirection(), tri.getEdge(1));
405 T denom = dot(p, tri.getEdge(0));
406 T signDenom = copysign(1, denom);
408 SGVec3<T> s = ray.getOrigin() - tri.getBaseVertex();
409 SGVec3<T> q = cross(s, tri.getEdge(0));
411 // t = 1/denom*dot(q, tri.getEdge(1));
412 // To avoid an expensive division we multiply by |denom|
413 T tDenom = signDenom*dot(q, tri.getEdge(1));
416 // For line segment we would test against
419 // with the original t. The multiplied test would read
420 // if (absDenom < tDenom)
423 T absDenom = fabs(denom);
424 T absDenomEps = absDenom*eps;
426 // T u = 1/denom*dot(p, s);
427 T u = signDenom*dot(p, s);
428 if (u < -absDenomEps)
430 // T v = 1/denom*dot(q, d);
433 T v = signDenom*dot(q, ray.getDirection());
434 if (v < -absDenomEps)
437 if (u + v > absDenom + absDenomEps)
440 // return if paralell ??? FIXME what if paralell and in plane?
441 // may be we are ok below than anyway??
442 if (absDenom <= SGLimits<T>::min())
446 // if we have survived here it could only happen with denom == 0
447 // that the point is already in plane. Then return the origin ...
448 if (SGLimitsd::min() < absDenom)
449 x += (tDenom/absDenom)*ray.getDirection();
456 intersects(const SGTriangle<T>& tri, const SGRay<T>& ray, T eps = 0)
458 // FIXME: for now just wrap the other method. When that has prooven
459 // well optimized, implement that special case
461 return intersects(dummy, tri, ray, eps);
466 // FIXME do not use that default argument later. Just for development now
467 intersects(SGVec3<T>& x, const SGTriangle<T>& tri, const SGLineSegment<T>& lineSegment, T eps = 0)
469 // See Tomas Akeniene - Moeller/Eric Haines: Real Time Rendering
471 // Method based on the observation that we are looking for a
472 // point x that can be expressed in terms of the triangle points
473 // x = v_0 + u*(v_1 - v_0) + v*(v_2 - v_0)
474 // with 0 <= u, v and u + v <= 1.
475 // OTOH it could be expressed in terms of the lineSegment
477 // Now we can compute u, v and t.
478 SGVec3<T> p = cross(lineSegment.getDirection(), tri.getEdge(1));
480 T denom = dot(p, tri.getEdge(0));
481 T signDenom = copysign(1, denom);
483 SGVec3<T> s = lineSegment.getStart() - tri.getBaseVertex();
484 SGVec3<T> q = cross(s, tri.getEdge(0));
486 // t = 1/denom*dot(q, tri.getEdge(1));
487 // To avoid an expensive division we multiply by |denom|
488 T tDenom = signDenom*dot(q, tri.getEdge(1));
491 // For line segment we would test against
494 // with the original t. The multiplied test reads
495 T absDenom = fabs(denom);
496 if (absDenom < tDenom)
499 // take the CPU accuracy in account
500 T absDenomEps = absDenom*eps;
502 // T u = 1/denom*dot(p, s);
503 T u = signDenom*dot(p, s);
504 if (u < -absDenomEps)
506 // T v = 1/denom*dot(q, d);
509 T v = signDenom*dot(q, lineSegment.getDirection());
510 if (v < -absDenomEps)
513 if (u + v > absDenom + absDenomEps)
516 // return if paralell ??? FIXME what if paralell and in plane?
517 // may be we are ok below than anyway??
518 if (absDenom <= SGLimits<T>::min())
521 x = lineSegment.getStart();
522 // if we have survived here it could only happen with denom == 0
523 // that the point is already in plane. Then return the origin ...
524 if (SGLimitsd::min() < absDenom)
525 x += (tDenom/absDenom)*lineSegment.getDirection();
532 intersects(const SGTriangle<T>& tri, const SGLineSegment<T>& lineSegment, T eps = 0)
534 // FIXME: for now just wrap the other method. When that has prooven
535 // well optimized, implement that special case
537 return intersects(dummy, tri, lineSegment, eps);
543 closestPoint(const SGTriangle<T>& tri, const SGVec3<T>& p)
545 // This method minimizes the distance function Q(u, v) = || p - x ||
546 // where x is a point in the trialgle given by the vertices v_i
547 // x = v_0 + u*(v_1 - v_0) + v*(v_2 - v_0)
548 // The theoretical analysis is somehow too long for a comment.
549 // May be it is sufficient to see that this code passes all the tests.
551 SGVec3<T> off = tri.getBaseVertex() - p;
552 T a = dot(tri.getEdge(0), tri.getEdge(0));
553 T b = dot(tri.getEdge(0), tri.getEdge(1));
554 T c = dot(tri.getEdge(1), tri.getEdge(1));
555 T d = dot(tri.getEdge(0), off);
556 T e = dot(tri.getEdge(1), off);
581 return tri.getBaseVertex() + tri.getEdge(0);
585 return tri.getBaseVertex() + u*tri.getEdge(0);
591 return tri.getBaseVertex();
592 } else if (c <= -e) {
595 return tri.getBaseVertex() + tri.getEdge(1);
599 return tri.getBaseVertex() + v*tri.getEdge(1);
607 return tri.getBaseVertex();
608 } else if (c <= -e) {
611 return tri.getBaseVertex() + tri.getEdge(1);
615 return tri.getBaseVertex() + v*tri.getEdge(1);
623 return tri.getBaseVertex();
624 } else if (a <= -d) {
627 return tri.getBaseVertex() + tri.getEdge(0);
631 return tri.getBaseVertex() + u*tri.getEdge(0);
635 if (det <= SGLimits<T>::min()) {
638 return tri.getBaseVertex();
643 return tri.getBaseVertex() + u*tri.getEdge(0) + v*tri.getEdge(1);
652 T numer = tmp1 - tmp0;
653 T denom = a - 2*b + c;
654 if (denom <= numer) {
657 return tri.getBaseVertex() + tri.getEdge(0);
661 return tri.getBaseVertex() + u*tri.getEdge(0) + v*tri.getEdge(1);
667 return tri.getBaseVertex() + tri.getEdge(1);
671 return tri.getBaseVertex();
675 return tri.getBaseVertex() + v*tri.getEdge(1);
683 T numer = tmp1 - tmp0;
684 T denom = a - 2*b + c;
685 if (denom <= numer) {
688 return tri.getBaseVertex() + tri.getEdge(1);
692 return tri.getBaseVertex() + u*tri.getEdge(0) + v*tri.getEdge(1);
698 return tri.getBaseVertex() + tri.getEdge(0);
702 return tri.getBaseVertex();
706 return tri.getBaseVertex() + u*tri.getEdge(0);
711 T numer = c + e - b - d;
715 return tri.getVertex(2);
717 T denom = a - 2*b + c;
718 if (denom <= numer) {
721 return tri.getBaseVertex() + tri.getEdge(0);
725 return tri.getBaseVertex() + u*tri.getEdge(0) + v*tri.getEdge(1);
733 closestPoint(const SGVec3<T>& p, const SGTriangle<T>& tri)
734 { return closestPoint(tri, p); }
736 template<typename T, typename T2>
738 intersects(const SGTriangle<T>& tri, const SGSphere<T2>& sphere)
740 // This method minimizes the distance function Q(u, v) = || p - x ||
741 // where x is a point in the trialgle given by the vertices v_i
742 // x = v_0 + u*(v_1 - v_0) + v*(v_2 - v_0)
743 // The theoretical analysis is somehow too long for a comment.
744 // May be it is sufficient to see that this code passes all the tests.
746 SGVec3<T> off = tri.getBaseVertex() - SGVec3<T>(sphere.getCenter());
747 T baseDist2 = dot(off, off);
748 T a = dot(tri.getEdge(0), tri.getEdge(0));
749 T b = dot(tri.getEdge(0), tri.getEdge(1));
750 T c = dot(tri.getEdge(1), tri.getEdge(1));
751 T d = dot(tri.getEdge(0), off);
752 T e = dot(tri.getEdge(1), off);
777 T sqrDist = a + 2*d + baseDist2;
778 return sqrDist <= sphere.getRadius2();
782 T sqrDist = d*u + baseDist2;
783 return sqrDist <= sphere.getRadius2();
789 return baseDist2 <= sphere.getRadius2();
790 } else if (c <= -e) {
793 T sqrDist = c + 2*e + baseDist2;
794 return sqrDist <= sphere.getRadius2();
798 T sqrDist = e*v + baseDist2;
799 return sqrDist <= sphere.getRadius2();
807 return baseDist2 <= sphere.getRadius2();
808 } else if (c <= -e) {
811 T sqrDist = c + 2*e + baseDist2;
812 return sqrDist <= sphere.getRadius2();
816 T sqrDist = e*v + baseDist2;
817 return sqrDist <= sphere.getRadius2();
825 return baseDist2 <= sphere.getRadius2();
826 } else if (a <= -d) {
829 T sqrDist = a + 2*d + baseDist2;
830 return sqrDist <= sphere.getRadius2();
834 T sqrDist = d*u + baseDist2;
835 return sqrDist <= sphere.getRadius2();
839 if (det <= SGLimits<T>::min()) {
840 // sqrDist = baseDist2;
843 return baseDist2 <= sphere.getRadius2();
848 T sqrDist = u*(a*u + b*v + 2*d) + v*(b*u + c*v + 2*e) + baseDist2;
849 return sqrDist <= sphere.getRadius2();
858 T numer = tmp1 - tmp0;
859 T denom = a - 2*b + c;
860 if (denom <= numer) {
863 T sqrDist = a + 2*d + baseDist2;
864 return sqrDist <= sphere.getRadius2();
868 T sqrDist = u*(a*u + b*v + 2*d) + v*(b*u + c*v + 2*e) + baseDist2;
869 return sqrDist <= sphere.getRadius2();
875 T sqrDist = c + 2*e + baseDist2;
876 return sqrDist <= sphere.getRadius2();
880 return baseDist2 <= sphere.getRadius2();
884 T sqrDist = e*v + baseDist2;
885 return sqrDist <= sphere.getRadius2();
893 T numer = tmp1 - tmp0;
894 T denom = a - 2*b + c;
895 if (denom <= numer) {
898 T sqrDist = c + 2*e + baseDist2;
899 return sqrDist <= sphere.getRadius2();
903 T sqrDist = u*(a*u + b*v + 2*d) + v*(b*u + c*v + 2*e)+baseDist2;
904 return sqrDist <= sphere.getRadius2();
910 T sqrDist = a + 2*d + baseDist2;
911 return sqrDist <= sphere.getRadius2();
913 // sqrDist = baseDist2;
916 return baseDist2 <= sphere.getRadius2();
920 T sqrDist = d*u + baseDist2;
921 return sqrDist <= sphere.getRadius2();
926 T numer = c + e - b - d;
930 T sqrDist = c + 2*e + baseDist2;
931 return sqrDist <= sphere.getRadius2();
933 T denom = a - 2*b + c;
934 if (denom <= numer) {
937 T sqrDist = a + 2*d + baseDist2;
938 return sqrDist <= sphere.getRadius2();
942 T sqrDist = u*(a*u + b*v + 2*d) + v*(b*u + c*v + 2*e) + baseDist2;
943 return sqrDist <= sphere.getRadius2();
949 template<typename T1, typename T2>
951 intersects(const SGSphere<T1>& sphere, const SGTriangle<T2>& tri)
952 { return intersects(tri, sphere); }
957 intersects(const SGVec3<T>& v, const SGSphere<T>& sphere)
961 return distSqr(v, sphere.getCenter()) <= sphere.getRadius2();
965 intersects(const SGSphere<T>& sphere, const SGVec3<T>& v)
966 { return intersects(v, sphere); }
971 intersects(const SGBox<T>& box, const SGLineSegment<T>& lineSegment)
973 // See Tomas Akeniene - Moeller/Eric Haines: Real Time Rendering
975 SGVec3<T> c = lineSegment.getCenter() - box.getCenter();
976 SGVec3<T> w = T(0.5)*lineSegment.getDirection();
977 SGVec3<T> v(fabs(w.x()), fabs(w.y()), fabs(w.z()));
978 SGVec3<T> h = T(0.5)*box.getSize();
980 if (fabs(c[0]) > v[0] + h[0])
982 if (fabs(c[1]) > v[1] + h[1])
984 if (fabs(c[2]) > v[2] + h[2])
987 if (fabs(c[1]*w[2] - c[2]*w[1]) > h[1]*v[2] + h[2]*v[1])
989 if (fabs(c[0]*w[2] - c[2]*w[0]) > h[0]*v[2] + h[2]*v[0])
991 if (fabs(c[0]*w[1] - c[1]*w[0]) > h[0]*v[1] + h[1]*v[0])
998 intersects(const SGLineSegment<T>& lineSegment, const SGBox<T>& box)
999 { return intersects(box, lineSegment); }
1001 template<typename T>
1003 intersects(const SGBox<T>& box, const SGRay<T>& ray)
1005 // See Tomas Akeniene - Moeller/Eric Haines: Real Time Rendering
1007 for (unsigned i = 0; i < 3; ++i) {
1008 T cMin = box.getMin()[i];
1009 T cMax = box.getMax()[i];
1011 T cOrigin = ray.getOrigin()[i];
1013 T cDir = ray.getDirection()[i];
1014 if (fabs(cDir) <= SGLimits<T>::min()) {
1021 T nearr = - SGLimits<T>::max();
1022 T farr = SGLimits<T>::max();
1024 T T1 = (cMin - cOrigin) / cDir;
1025 T T2 = (cMax - cOrigin) / cDir;
1026 if (T1 > T2) std::swap (T1, T2);/* since T1 intersection with near plane */
1027 if (T1 > nearr) nearr = T1; /* want largest Tnear */
1028 if (T2 < farr) farr = T2; /* want smallest Tfarr */
1029 if (nearr > farr) // farr box is missed
1031 if (farr < 0) // box is behind ray
1037 // make it symmetric
1038 template<typename T>
1040 intersects(const SGRay<T>& ray, const SGBox<T>& box)
1041 { return intersects(box, ray); }
1043 template<typename T1, typename T2>
1045 intersects(const SGBox<T1>& box1, const SGBox<T2>& box2)
1047 if (box2.getMax()[0] < box1.getMin()[0])
1049 if (box1.getMax()[0] < box2.getMin()[0])
1052 if (box2.getMax()[1] < box1.getMin()[1])
1054 if (box1.getMax()[1] < box2.getMin()[1])
1057 if (box2.getMax()[2] < box1.getMin()[2])
1059 if (box1.getMax()[2] < box2.getMin()[2])