Geometry

00001 #ifndef NMWR_GB_ALGEBRAIC_PRIMITIVES_2D_H
00002 #define NMWR_GB_ALGEBRAIC_PRIMITIVES_2D_H
00003 
00004 
00005 // $LICENSE
00006 
00007 //----------------------------------------------------------------
00008 //
00009 // specialize dimension_dependent_primitives<POINT>
00010 // to dimension_dependent_primitives_2d<POINT> for 2D points.
00011 //
00012 //----------------------------------------------------------------
00013 
00014 #include <algorithm> // STL: swap
00015 
00016 
00017 template<class POINT>
00018 struct dimension_dependent_primitives_2d 
00019   : public virtual basic_algebraic_primitives<POINT> 
00020 {
00021 
00022   static double det2(const POINT& p1, const POINT& p2)
00023     {
00024       int l = point_traits<POINT>::LowerIndex(p1);
00025       // IOMgr::Info() << "lower index: " << l << endl;
00026       // IOMgr::Info() << p1 << " " << p2 << endl;
00027       return (p1[l]*p2[l+1]-p1[l+1]*p2[l]);
00028     }
00029 
00030   static double signed_triangle_area(const POINT& p1, const POINT& p2, const POINT& p3)
00031     { 
00032       // IOMgr::Info() << "signed triangle area: " 
00033       //                    << p1 << " " << p2 << " " << p3 << endl;
00034       return (0.5 * det2(p1-p2,p1-p3)); }
00035 
00036   static double triangle_area(const POINT& p1, const POINT& p2, const POINT& p3)
00037     {
00038       double a(signed_triangle_area(p1,p2,p3));
00039       return(a >= 0 ? a : -a);
00040     }
00041 
00042   static POINT normal_with_same_length(const POINT& p)
00043     {
00044       int l = LowerIndex(p);
00045       return POINT(-p[l+1],p[l]);
00046     }
00047 
00048   static POINT normed_normal(const POINT& p)
00049     {
00050       POINT q(normal_with_same_length(p));
00051       return (q/norm_2(q));
00052     }
00053 
00054   static void transpose(POINT & p, POINT & q) {
00055     std::swap(y(p),x(q));
00056   }
00057 
00058   // solve (A1,A2)*X = b
00059   static void solve2(POINT const& A1, POINT const& A2,
00060                      POINT      & X,  POINT const& b)
00061   {
00062     scalar d = det2(A1,A2);
00063     x(X) = ( y(A2) * x(b) - x(A2) * y(b))/d;
00064     y(X) = (-y(A1) * x(b) + x(A1) * y(b))/d;
00065   }
00066 
00067   static void solveA(POINT  A1, POINT  A2,
00068                      POINT& X , POINT  b)
00069   {
00070     if(fabs(x(A1)) < fabs(y(A1))) {
00071       std::swap(x(A1),y(A1));
00072       std::swap(x(A2),y(A2));
00073       std::swap(x(b) ,y(b));
00074     }
00075     // Lu = b, l_21 = a21/a11 = y(A1)/x(A1)
00076     scalar u1 = x(b);
00077     scalar u2 = y(b)-(y(A1)/x(A1))*u1;
00078     // Rx = u, r11 = x(A1), r12 = x(A2), r22 = y(A2) - y(A1)x(A1)^-1x(A2)
00079     y(X) = u2/(y(A2) - (y(A1)/x(A1))*x(A2));
00080     x(X) = (u1 - x(A2)*y(X))/x(A1);
00081   }
00082 };
00083 
00084 #endif
Geometry/primitives2d.h
