planeh3.h

很多二维 三维几何计算算法 C++ 类库

PlaneH3<R>::d() const
{ return get(base).e3; }

template < class R >
CGAL_KERNEL_INLINE
typename PlaneH3<R>::Line_3
PlaneH3<R>::perpendicular_line(const typename PlaneH3<R>::Point_3& p) const
{ return Line_3( p, orthogonal_direction() ); }

template < class R >
CGAL_KERNEL_INLINE
typename PlaneH3<R>::Plane_3
PlaneH3<R>::opposite() const
{ return PlaneH3<R>(-a(), -b(), -c(), -d() ); }

template < class R >
CGAL_KERNEL_INLINE
typename PlaneH3<R>::Point_3
PlaneH3<R>::projection(const typename PlaneH3<R>::Point_3& p) const
{ return _projection( p, *this ); }

template < class R >
CGAL_KERNEL_INLINE
typename PlaneH3<R>::Point_3
PlaneH3<R>::point() const
{
  const RT RT0(0);
  if ( a() != RT0 )
  {
      return Point_3( -d(), RT0, RT0, a() );
  }
  if ( b() != RT0 )
  {
      return Point_3( RT0, -d(), RT0, b() );
  }
  CGAL_kernel_assertion ( c() != RT0);
  return Point_3( RT0, RT0, -d(), c() );
}

template < class R >
CGAL_KERNEL_INLINE
typename PlaneH3<R>::Vector_3
PlaneH3<R>::base1() const
{
 // point():
 // a() != RT0 : Point_3( -d(), RT0, RT0, a() );
 // b() != RT0 : Point_3( RT0, -d(), RT0, b() );
 //            : Point_3( RT0, RT0, -d(), c() );
 // point1():
 // a() != RT0 : Point_3( -b()-d(), a(), RT0, a() );
 // b() != RT0 : Point_3( RT0, -c()-d(), b(), b() );
 //            : Point_3( c(), RT0, -a()-d(), c() );

  const RT RT0(0);
  if ( a() != RT0 )
  {
      return Vector_3( -b(), a(), RT0, a() );
  }
  if ( b() != RT0 )
  {
      return Vector_3( RT0, -c(), b(), b() );
  }
  CGAL_kernel_assertion ( c() != RT(0) );
  return Vector_3( c(), RT0, -a(), c() );
}

template < class R >
inline
typename PlaneH3<R>::Vector_3
PlaneH3<R>::base2() const
{
  Vector_3 a = orthogonal_vector();
  Vector_3 b = base1();
  return Vector_3(a.hy()*b.hz() - a.hz()*b.hy(),
                         a.hz()*b.hx() - a.hx()*b.hz(),
                         a.hx()*b.hy() - a.hy()*b.hx(),
                         a.hw()*b.hw() );
}
// Actually, the following should work, but bcc doesn't like it:
// { return cross_product( orthogonal_vector(), base1() ); }


template < class R >
inline
typename PlaneH3<R>::Point_3
PlaneH3<R>::point1() const
{ return point() + base1(); }

template < class R >
inline
typename PlaneH3<R>::Point_3
PlaneH3<R>::point2() const
{ return point() + base2(); }

template < class R >
inline
typename PlaneH3<R>::Direction_3
PlaneH3<R>::orthogonal_direction() const
{ return Direction_3(a(), b(), c() ); }

template < class R >
inline
typename PlaneH3<R>::Vector_3
PlaneH3<R>::orthogonal_vector() const
{ return Vector_3(a(), b(), c() ); }

template < class R >
bool
PlaneH3<R>::is_degenerate() const
{
 const RT RT0(0);
 return ( (a() == RT0 ) && (b() == RT0 ) && (c() == RT0 ) );
}

template < class R >
bool
PlaneH3<R>::has_on_positive_side( const typename PlaneH3<R>::Point_3& p) const
{
 return (a()*p.hx() + b()*p.hy() + c()*p.hz() + d()*p.hw() > RT(0) );
}

template < class R >
bool
PlaneH3<R>::has_on_negative_side( const typename PlaneH3<R>::Point_3& p) const
{
 return (a()*p.hx() + b()*p.hy() + c()*p.hz() + d()*p.hw() < RT(0) );
}


template < class R >
bool
PlaneH3<R>::has_on( const typename PlaneH3<R>::Point_3& p) const
{
 return (a()*p.hx() + b()*p.hy() + c()*p.hz() + d()*p.hw() == RT(0) );
}

template < class R >
bool
PlaneH3<R>::has_on( const typename PlaneH3<R>::Line_3& l) const
{
 Point_3   p   = l.point();
 Vector_3  ld  = l.direction().to_vector();
 Vector_3  ov  = orthogonal_vector();

 return (  ( a()*p.hx() + b()*p.hy() + c()*p.hz() + d()*p.hw()   == RT(0) )
         &&( ld.hx()*ov.hx() + ld.hy()*ov.hy() + ld.hz()*ov.hz() == RT(0) ) );
}

template < class R >
Oriented_side
PlaneH3<R>::oriented_side( const typename PlaneH3<R>::Point_3& p) const
{
 RT value = a()*p.hx() + b()*p.hy() + c()*p.hz() + d()*p.hw() ;
 if (value > RT(0) )
 {
    return ON_POSITIVE_SIDE;
 }
 else
 {
    return
    (value < RT(0) ) ? ON_NEGATIVE_SIDE : ON_ORIENTED_BOUNDARY;
 }
}


template < class R >
bool
PlaneH3<R>::operator==(const PlaneH3<R>& l) const
{
 if (  (a() * l.d() != l.a() * d() )
     ||(b() * l.d() != l.b() * d() )
     ||(c() * l.d() != l.c() * d() ) )
 {
    return false;
 }
 int sd  = static_cast<int>(CGAL_NTS sign(d()));
 int sld = static_cast<int>(CGAL_NTS sign(l.d()));
 if ( sd == sld )
 {
    if (sd == 0)
    {
        return (  (a()*l.b() == b()*l.a() )
                &&(a()*l.c() == c()*l.a() )
                &&(b()*l.c() == c()*l.b() )
                &&(CGAL_NTS sign(a() )== CGAL_NTS sign( l.a() ))
                &&(CGAL_NTS sign(b() )== CGAL_NTS sign( l.b() ))
                &&(CGAL_NTS sign(c() )== CGAL_NTS sign( l.c() )) );
    }
    else
    {
        return true;
    }
 }
 else
 {
    return false;
 }
}

template < class R >
typename PlaneH3<R>::Aff_transformation_3
PlaneH3<R>::transform_to_2d() const
{
  const RT  RT0(0);
  const RT  RT1(1);
  Vector_3 nov = orthogonal_vector();
  Vector_3 e1v = point1()-point() ;
  Vector_3 e2v = point2()-point() ;
  RT orthohx = nov.hx();
  RT orthohy = nov.hy();
  RT orthohz = nov.hz();
  RT e1phx   = e1v.hx();
  RT e1phy   = e1v.hy();
  RT e1phz   = e1v.hz();
  RT e2phx   = e2v.hx();
  RT e2phy   = e2v.hy();
  RT e2phz   = e2v.hz();

  RT t11 =  -( orthohy*e2phz - orthohz*e2phy );
  RT t12 =   ( orthohx*e2phz - orthohz*e2phx );
  RT t13 =  -( orthohx*e2phy - orthohy*e2phx );

  RT t21 =   ( orthohy*e1phz - orthohz*e1phy );
  RT t22 =  -( orthohx*e1phz - orthohz*e1phx );
  RT t23 =   ( orthohx*e1phy - orthohy*e1phx );

  RT t31 =   ( e1phy*e2phz - e1phz*e2phy );
  RT t32 =  -( e1phx*e2phz - e1phz*e2phx );
  RT t33 =   ( e1phx*e2phy - e1phy*e2phx );

  RT scale = det3x3_by_formula( orthohx, orthohy, orthohz,
                                     e1phx,   e1phy,   e1phz,
                                     e2phx,   e2phy,   e2phz );

  Aff_transformation_3
     point_to_origin(TRANSLATION,  - ( point() - ORIGIN ) );
  Aff_transformation_3
     rotate_and_more( t11,    t12,   t13,   RT0,
                      t21,    t22,   t23,   RT0,
                      t31,    t32,   t33,   RT0,
                                            scale);

  Point_3 ortho( orthohx, orthohy, orthohz );
  Point_3 e1p( e1phx, e1phy, e1phz );
  Point_3 e2p( e2phx, e2phy, e2phz );
  CGAL_kernel_assertion((   ortho.transform(rotate_and_more)
        == Point_3( RT(0), RT(0), RT(1)) ));
  CGAL_kernel_assertion((   e1p.transform(rotate_and_more)
        == Point_3( RT(1), RT(0), RT(0)) ));
  CGAL_kernel_assertion((   e2p.transform(rotate_and_more)
        == Point_3( RT(0), RT(1), RT(0)) ));

  return  rotate_and_more * point_to_origin;
}

template < class R >
CGAL_KERNEL_INLINE
typename PlaneH3<R>::Point_2
PlaneH3<R>::to_2d(const typename PlaneH3<R>::Point_3& p) const
{
  Point_3 tp = p.transform( transform_to_2d() );
  return Point_2( tp.hx(), tp.hy(), tp.hw());
}


template < class R >
CGAL_KERNEL_INLINE
typename PlaneH3<R>::Point_3
PlaneH3<R>::to_3d(const typename PlaneH3<R>::Point_2& p)  const
{
  Point_3 hp( p.hx(), p.hy(), RT(0.0), p.hw());
  return hp.transform( transform_to_2d().inverse() );
}

CGAL_END_NAMESPACE

#endif  // CGAL_PLANEH3_H