Geometry

00001 #ifndef GRAL_GB_GEOMETRY_TEST_PRIMITIVES3D_C
00002 #define GRAL_GB_GEOMETRY_TEST_PRIMITIVES3D_C
00003 
00004 #include "Geometry/test-primitives3d.h"
00005 #include "Geometry/algebraic-primitives.h"
00006 #include "Geometry/geometric-primitives.h"
00007 #include "Geometry/primitives3d.h"
00008 
00009 #include "Utility/pre-post-conditions.h"
00010 
00011 
00012 template<class P>
00013 void check_intersections(P const& s0, P const& s1,
00014                          P const& t0, P const& t1, P const& t2,
00015                          ostream& out)
00016 {
00017 
00018   typedef segment<P>   segment3;
00019   typedef ray<P>       ray3;
00020   typedef triangle<P>  triangle3;
00021   typedef intersection_segment_triangle<segment3, triangle3> intersection_st;
00022   typedef intersection_ray_triangle<ray3, triangle3>         intersection_rt;
00023 
00024   segment3  S(s0, s1);
00025   triangle3 T(t0, t1, t2);
00026   ray3      R(s0, s1 -s0);
00027   intersection_st  IS(S,T);
00028   intersection_rt  IR(R,T);
00029   
00030   out << "\n" 
00031       << "Segment:  [" << S.p0() << "," << S.p1() << "]\n"
00032       << "Triangle: [" << T.p0() << "," << T.p1() << "," << T.p2() << "]\n"
00033       << "Ray:      [" << R.p0() << "," << R.dir() << "]\n";
00034 
00035   // segment intersection
00036   out << "segment does" 
00037       << (IS.segment_intersects_plane() ? " " : " not ")
00038       << "intersect plane\n"
00039       << "segment does" 
00040       << (IS.segment_intersects_triangle() ? " " : " not ")
00041       << "intersect triangle\n";
00042   if(IS.segment_intersects_plane())
00043     out << "plane intersection: " << IS.intersection() << '\n'; 
00044 
00045   // ray intersection
00046   out << "line does" 
00047       << (IR.ray_intersects_plane() ? " " : " not ")
00048       << "intersect plane\n"
00049       << "ray does" 
00050       << (IR.ray_intersects_triangle() ? " " : " not ")
00051       << "intersect triangle\n";
00052   if(IR.ray_intersects_plane())
00053     out << "plane intersection: " << IR.intersection() << '\n'; 
00054 }
00055 
00056 
00057 template<class POINT>
00058 void test_algebraic_primitives3d<POINT>::do_tests(std::ostream& out)
00059 {
00060   typedef point_traits<POINT>         pt;
00061   typedef dimension_dependent_primitives_3d<POINT> ap;
00062 
00063   POINT b = pt::Origin(3);
00064   int li = pt::LowerIndex(b); // should be the same for all points
00065 
00066   POINT e[3];
00067   for(int i = 0; i < 3; ++i) {
00068     e[i] = pt::Origin(3);
00069     e[i][li + i] = 1;
00070     b[li+i] = 1;
00071   }
00072 
00073   for(int i = 0; i < 3; ++i) {
00074     int j = (i < 2 ? i+1 : 0);
00075     int k = (i > 0 ? i-1 : 2);
00076     out << "ap::det3(e["<<i<<"],e["<<k<<"],e["<<j<<"]) = "
00077         << ap::det3(e[i],e[k],e[j]) << '\n'
00078         << "ap::det3(e["<<i<<"],e["<<j<<"],e["<<k<<"]) = "
00079         << ap::det3(e[i],e[j],e[k]) << '\n'
00080         << "ap::vectorproduct(e["<<i<<"],e["<<j<<"]) = "
00081         << ap::vectorproduct(e[i],e[j]) << '\n'
00082         << "ap::vectorproduct(e["<<i<<"],e["<<k<<"]) = "
00083         << ap::vectorproduct(e[i],e[k]) << '\n';
00084 
00085   }
00086 
00087   POINT x = pt::Origin(3);
00088   ap::solve(e[0],e[1],e[2],x,b);
00089   out << x << '\n';
00090   ap::solve(e[0],e[2],e[1],x,b);
00091   out << x << '\n';
00092 
00093   POINT inv[3];
00094   for(int i = 0; i < 3; ++i) {
00095     inv[i] = pt::Origin(3);
00096   }
00097   out << "ap::inverse(e[0],e[1],e[2]) = ";
00098   ap::inverse(e[0],e[1],e[2],inv[0],inv[1],inv[2]);
00099   out << inv[0] << "  " << inv[1] << "  " << inv[2] << '\n';
00100 
00101   out << "ap::inverse(e[0],e[2],e[1]) = ";
00102   ap::inverse(e[0],e[2],e[1],inv[0],inv[1],inv[2]);
00103   out << inv[0] << "  " << inv[1] << "  " << inv[2] << '\n';
00104 
00105   out<< "ap::condition(e[0],e[1],e[2],ap::Norm_frobenius()) = "
00106      <<  ap::condition(e[0],e[1],e[2],ap::Norm_frobenius()) << "\n";
00107   out<< "ap::condition(e[0],e[1],e[2],ap::Norm_1()) = "
00108      <<  ap::condition(e[0],e[1],e[2],ap::Norm_1()) << "\n";
00109   out<< "ap::condition(e[0],e[1],e[2],ap::Norm_infinity()) = "
00110      <<  ap::condition(e[0],e[1],e[2],ap::Norm_infinity()) << "\n";
00111 
00112   POINT a[3];
00113   for(int i = 0; i < 3; ++i) {
00114     a[i] = pt::Origin(3);
00115     a[i][li+i] = i+1;
00116   }
00117 
00118 
00119   out << "ap::inverse(a[0],a[1],a[2]) = ";
00120   ap::inverse(a[0],a[1],a[2],inv[0],inv[1],inv[2]);
00121   out << inv[0] << "  " << inv[1] << "  " << inv[2] << '\n';
00122   out << "ap::inverse(a[0],a[1],a[2]) = ";
00123   ap::inverse(a[0],a[2],a[1],inv[0],inv[1],inv[2]);
00124   out << inv[0] << "  " << inv[1] << "  " << inv[2] << '\n';
00125 
00126   out << "ap::matrixnorm_frobenius(a[0],a[1],a[2]) = "
00127       << ap::matrixnorm_frobenius(a[0],a[1],a[2]) << '\n'
00128       << "ap::matrixnorm_1(a[0],a[1],a[2]) = "
00129       << ap::matrixnorm_1(a[0],a[1],a[2]) << '\n'
00130       << "ap::matrixnorm_infinity(a[0],a[1],a[2]) = "
00131       << ap::matrixnorm_infinity(a[0],a[1],a[2]) << '\n';
00132 
00133   out<< "ap::condition(a[0],a[1],a[2],ap::Norm_frobenius()) = "
00134      <<  ap::condition(a[0],a[1],a[2],ap::Norm_frobenius()) << "\n";
00135   out<< "ap::condition(a[0],a[1],a[2],ap::Norm_1()) = "
00136      <<  ap::condition(a[0],a[1],a[2],ap::Norm_1()) << "\n";
00137   out<< "ap::condition(a[0],a[1],a[2],ap::Norm_infinity()) = "
00138      <<  ap::condition(a[0],a[1],a[2],ap::Norm_infinity()) << "\n";
00139 
00140 
00141 
00142   typedef segment<POINT>   segment3;
00143   typedef ray<POINT>       ray3;
00144   typedef triangle<POINT>  triangle3;
00145   typedef intersection_segment_triangle<segment3, triangle3> intersection_t;
00146   typedef intersection_ray_triangle<ray3, triangle3>         intersection_rt;
00147 
00148   //------- intersection test #1 -----------
00149 
00150   // S = [ (0.25,0.25,-1), (0.25,0.25, 1)] 
00151   // T = [ (0,0,0), (1,0,0), (0,1,1)]
00152   // intersection (0.25,0.25,0)
00153 
00154   POINT p[2], q[3];
00155   for(int i = 0; i < 2; ++i) {
00156     p[i] = pt::Origin(3);
00157   }  
00158   for(int i = 0; i < 3; ++i) {
00159     q[i] = pt::Origin(3);
00160   }  
00161 
00162   pt::z(p[1]) = 1;  
00163   pt::z(p[0]) = -1;
00164   pt::x(p[0]) = pt::x(p[1]) = pt::y(p[0]) = pt::y(p[1]) = 0.25;
00165 
00166   pt::x(q[1]) = 1;
00167   pt::y(q[2]) = 1;
00168 
00169   check_intersections(p[0],p[1],
00170                       q[0],q[1],q[2], out);
00171 
00172 
00173   //------- intersection test #2 -----------
00174 
00175   // S = [ (0.25,0.25, 0.5), (0.25,0.25, 1)] 
00176   // T = [ (0,0,0), (1,0,0), (0,1,1)]
00177   // intersection does not exist
00178 
00179   // segment: uses const& of POINT
00180   pt::z(p[0]) = 0.5; 
00181   check_intersections(p[0],p[1],
00182                       q[0],q[1],q[2], out);
00183 
00184 
00185 
00186   //------- intersection test #3 -----------
00187   // S = [(0.5,0.5,0), (0.5,0.5,1)]
00188   // T = [(-0.5,-0.5,0), (-0.5, 0.5,1), (-0.5,-0.5,1)
00189   // intersection does not exist
00190 
00191 
00192   pt::x(p[0]) = pt::x(p[1]) = pt::y(p[0]) = pt::y(p[1]) =  0.5;
00193   pt::z(p[0]) = 0; pt::z(p[1]) = 1;
00194   pt::x(q[0]) = pt::x(q[1]) = pt::x(q[2]) = -0.5;
00195   pt::y(q[1]) = 0.5;
00196   pt::y(q[0]) = pt::y(q[2]) = -0.5;
00197   pt::z(q[0]) = 0;
00198   pt::z(q[1]) = pt::z(q[2]) = 1;
00199 
00200   check_intersections(p[0],p[1],
00201                       q[0],q[1],q[2], out);
00202 
00203   //------- intersection test #4 -----------
00204   // segmen/ray parallel to triangle
00205 
00206   {
00207     POINT p0(1,1,2), p1(1,0,2);
00208     POINT q0(0,0,0), q1(1,0,0), q2(0,1,0);
00209     check_intersections(p0,p1,
00210                         q0,q1,q2,out);
00211   }
00212   //------- intersection test #5 -----------
00213   // intersection on edge of triangle
00214   {
00215     POINT p0(.5,.5,-1), p1(.5,.5,1);
00216     POINT q0(0,0,0), q1(1,0,0), q2(0,1,0);
00217     check_intersections(p0,p1,
00218                         q0,q1,q2,out);
00219   }
00220 
00221 
00222 }
00223 
00224 #endif
Geometry/internal/test-primitives3d.C
