line.h

任意给定三维空间的点集

/*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
 * line.h -
 *     Represents a line in the plane. Adapted from CGAL.
\*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*/

#ifndef  __LINE__H
#define  __LINE__H

#ifndef  ___XGENERIC__H
//#error must include file "generic.h" before this file
#endif  /* ___XGENERIC__H */

template <class NT>
int Msign(const NT &x)
{
    return (x < NT(0)) ? -1: (x > NT(0)) ? 1 : 0;
}

#define  DASSERT( X, Y)   {if  ( !(X) ) { fprintf( stderr, Y ); assert(X); }}
typedef  double  real;
typedef  Vector3d   point;

#define  DOUBLE_EPS  1e-10

class  segment;
class  line;
//class  dpoint;

/* compute inteersection of two lines */
point  intersect_lines( const line  & l1, const line  & l2, 
                        bool  & f_intersect );
point * intersect_lines_ptr( const line  & l1, const line  & l2, 
                             bool  & f_intersect );


class  line {
private:
    real   arr[ 3 ];
    double   d_arr[ 3 ];
    int  _id, _true_id;
    bool   fSplitter;

    inline void   init( const real  & a, const real  & b, const real  & c,
                        int  id = -1 ) 
    { 
#ifdef  EXACT_ARITHMETIC
        /*
        Integer   lc = a.denominator() * b.denominator() * c.denominator();

        Integer   an( a.numerator() * (lc / a.denominator()) );
        Integer   bn( b.numerator() * (lc / b.denominator()) );
        Integer   cn( c.numerator() * (lc / c.denominator()) );
        
        // We can even do better, by dividing in the gcd of an, bn, cn 
        Integer   gc( gcd( gcd( an, bn ), cn ) );

        if  ( gc == 0 )
            gc = 1;

        arr[ 0 ] = real( an / gc );
        arr[ 1 ] = real( bn / gc );
        arr[ 2 ] = real( cn / gc );
        */
        arr[ 0 ] = a ;
        arr[ 1 ] = b ;
        arr[ 2 ] = c ;
        d_arr[ 0 ] = ( arr[ 0 ] );
        d_arr[ 1 ] = ( arr[ 1 ] );
        d_arr[ 2 ] = ( arr[ 2 ] );
#else           
        arr[ 0 ] = a;
        arr[ 1 ] = b;
        arr[ 2 ] = c;
#endif

        fSplitter = false;
        _id = id;
        _true_id = -1;
    }

public:
    const real  & a() const
    {
        return  arr[ 0 ];
    }
    const real  & b() const
    {
        return  arr[ 1 ];
    }
    const real  & c() const
    {
        return  arr[ 2 ];
    }

    const double  & d_a() const
    {
        return  d_arr[ 0 ];
    }
    const double  & d_b() const
    {
        return  d_arr[ 1 ];
    }
    const double  & d_c() const
    {
        return  d_arr[ 2 ];
    }

    inline void   compress()  const
    {
#ifdef  EXACT_ARITHMETIC
        Integer   lc = a().denominator() * b().denominator() * c().denominator();

        Integer   an( a().numerator() * (lc / a().denominator()) );
        Integer   bn( b().numerator() * (lc / b().denominator()) );
        Integer   cn( c().numerator() * (lc / c().denominator()) );
        
        // We can even do better, by dividing in the gcd of an, bn, cn 
        Integer   gc( gcd( gcd( an, bn ), cn ) );

        if  ( gc == 0 )
            gc = 1;

        *(real *)&(arr[ 0 ]) = real( an / gc );
        *(real *)&(arr[ 1 ]) = real( bn / gc );
        *(real *)&(arr[ 2 ]) = real( cn / gc );
#endif //  EXACT_ARITHMETIC
    }


    line() 
    {
        init( 0, 0, 0 );
    }

    line( const point  & p, const point  & q, int  id = -1 ) 
    {
        init( p.y - q.y,
              q.x - p.x,
              p.x*q.y - p.y*q.x, id );
    }

    
    line( const real  & _a, const real  & _b, const real  & _c, int  id = -1 )
    {
        init( _a, _b, _c, id );
    }

    void  init( const line  & l)
    {
        arr[ 0 ] = l.arr[ 0 ];
        arr[ 1 ] = l.arr[ 1 ];
        arr[ 2 ] = l.arr[ 2 ];        
        d_arr[ 0 ] = l.d_arr[ 0 ];
        d_arr[ 1 ] = l.d_arr[ 1 ];
        d_arr[ 2 ] = l.d_arr[ 2 ];
        _id = l._id;
        _true_id = l._true_id;
        fSplitter = l.fSplitter;
    }
    line( const line  & l ) 
    {
        init( l );
    }

    line  & operator=(const line  &l)
    {
        arr[ 0 ] = l.arr[ 0 ];
        arr[ 1 ] = l.arr[ 1 ];
        arr[ 2 ] = l.arr[ 2 ];        
        d_arr[ 0 ] = l.d_arr[ 0 ];
        d_arr[ 1 ] = l.d_arr[ 1 ];
        d_arr[ 2 ] = l.d_arr[ 2 ];
        _id = l._id;
        fSplitter = l.fSplitter;
        _true_id = l._true_id;
        
        return  *this;
    }

    int  id() const
    {
        return  _id;
    }

    int  true_id() const  {
        return  _true_id;
    }
    void  set_true_id( int  val ) {
        _true_id = val;
    }

    void   _setID( int  id_num ) {
        _id = id_num;
    }
    bool operator==(const line  &l) const
    {
        int  sc, slc;

        sc = slc = 0;
        //return ( ( a() == l.a() )  &&  ( b() == l.b() )
        //      &&  ( c() == l.c() ) );
        if  ( ( l._id > 0 )  &&  ( l._id == _id ) ) {
            printf( "bogi!\n" );
            return  true;
        }
        
        if  ( ( a() * l.c() != l.a() * c() )
              ||  ( b() * l.c() != l.b() * c() ) )
            return false;
        sc  = Msign(c());
        slc = Msign(l.c());
        if ( sc == slc ) {
            return (sc == 0) ?  ( a()*l.b() ==  b()*l.a() )
            && (Msign(a() ) == Msign( l.a() ))
                 && (Msign(b() ) == Msign( l.b() ))
            : true;
        } else { 
            return false;
        } 
    }

    bool isRealEqual(const line  &l) const
    {
        int  sc, slc;
        double  res;

        sc = slc = 0;
        //return ( ( a() == l.a() )  &&  ( b() == l.b() )
        //      &&  ( c() == l.c() ) );
        if  ( ( l._id > 0 )  &&  ( l._id == _id ) ) 
            return  true;
        
        res = d_a() * l.d_c() - l.d_a() * d_c();
        if  ( res < -DOUBLE_EPS )
            return  false;
        if  ( res > DOUBLE_EPS )
            return  false;

        res = ( d_b() * l.d_c() - l.d_b() * d_c() );
        if  ( res < -DOUBLE_EPS )
            return  false;
        if  ( res > DOUBLE_EPS )
            return  false;

        return  *this == l;
    }

    bool  isSameUpToSign( const line  & l ) {
        return  ( ( ( l.a() == a() )  &&  ( l.b() == b() )  
                    &&  ( l.c() == c() ) )
                  ||  ( ( l.a() == -a() )  &&  ( l.b() == -b() )
                        &&  ( l.c() == -c() ) ) );
    }

    bool  isSame( const  line  & l ) const {
        if (  (a() * l.c() != l.a() * c())
              ||(b() * l.c() != l.b() * c()) )
            {
                return false;
            }
        int sc  = Msign(c());
        int slc = Msign(l.c());
        if ( sc == slc ) {
            return (sc == 0) ?  ( a()*l.b() ==  b()*l.a() )
                : true;
        } else {
            if  ( ( sc == -slc )  &&  ( sc != 0 ) )
                return  true;
            else
                return  false;
        }
    }

 
    bool operator!=( const line  & l ) const
    {
        return !(*this == l);
    }
    
    void  setSplitter( bool  _fSplitter ) {
        fSplitter = _fSplitter;
    }

    bool  isSplitter() const {
        return  fSplitter;
    }

void   print( int  ind )  const
{
    //print_spaces( ind );
    printf( "          " );
	printf( "l:(a:%g, b:%g, c:%g) id: %d  rid: %d\n",
                ( a() ), ( b() ),
                ( c() ),
                _id, _true_id                  
                );
}


    real  x_at_y( const real  &y ) const
    {
        
        DASSERT( (a() != real(0)),
           "Line::x_at_y(const FT &y) is undefined for horizontal line" );
        return ( -b()*y - c() ) / a();
    }


    real  y_at_x( const real  &x ) const
    {
        if  ( b() == 0 ) {
            print( 20 );
            printf( "what?\n" );
            assert( false );
    	}

        DASSERT( (b() != real(0)),
                 "Line::yat_x(const FT &y) is undefined for vertical line");
        return ( ( (-a()*x) - c() ) / b()) ;
    }

    double  d_y_at_x( const double  &x ) const
    {
        assert( ( d_b() != 0 ) );
        return ( ( (-d_a()*x) - d_c() ) / d_b()) ;
    }

    line  perpendicular( const point  &p ) const
    {
        return line( -b() , a() , b() * p.x - a() * p.y  );
    }
    

    line  opposite() const
    {        
        return line( -a(), -b(), -c() );
    }

    point  direction( void ) const
    {
        return  point( b(), -a(), 0 );
    }

    point  point_on(int i) const
    {
        if (i == 0){
            return is_vertical() ? point( (-b()-c())/a(), 1, 0 )
                : point( 1, -(a()+c())/b(), 0);
        }
        if (i == 1) {
            return is_vertical() ? point( (-b()-c())/a() + b(),
                                          1 - a(), 0 )
                       : point( 1 + b(),
                                -(a()+c())/b() - a(), 0 );
        }
        // we add i times the direction
        if  ( is_vertical() ) {
            return point( (-b()-c())/a() + i*b(),
                          real(1) - i*a(), 0 );
        }
        return point( 1 + real(i)*b(),
                      -(a()+c())/b() - i*a(), 0 );
    }

    point point_on() const
    {
        return is_vertical() ? point( (-b()-c())/a(), real(1), 0 )
            : point( real(1), -(a()+c())/b(), 0 );
    }

    point projection(const point &p) const
    {
        if (is_horizontal()) {
            return point(p.x, -c()/b(), 0 );
        }
        
        if (is_vertical()) {
            return point( -c()/a(), p.y, 0 );
        }
        real ab = a()/b(), ba = b()/a(), ca = c()/a();
        real y = ( -p.x + ab*p.y - ca ) / ( ba + ab );
     
        return point(-ba * y - ca, y, 0);
    }


    #define  ON_ORIENTED_BOUNDARY   0
    #define  ON_POSITIVE_SIDE  1
    #define  ON_NEGATIVE_SIDE - 1
    int     oriented_side( const point &p) const
    {
        real n = a()*p.x + b()*p.y + c();
        return (n == real(0)) ? ON_ORIENTED_BOUNDARY
            : ( (n > real(0)) ? ON_POSITIVE_SIDE
                : ON_NEGATIVE_SIDE);
    }

    point  intersectx( const line  & l, bool  & f_intersect ) const 
    {
        return  intersect_lines( (*this), l, f_intersect );
    }

    //int  get_point_signx( const point  & p ) const;

    int  get_point_sign( const point  & p ) const;
    //int  get_point_sign( const float_point  & p ) const;
    //int  get_dpoint_sign( const dpoint  & p ) const;

    real  get_point_val( const point  & p ) const
    {
        return ( a()*p.x + b()*p.y + c() );
    }

    bool has_on_boundary(const point &p) const
    {
        return (a()*p.x + b()*p.y + c()) == real(0);
    }

    inline bool has_on_positive_side( const point &p) const
    {
        return (a()*p.x + b()*p.y + c()) >  real(0);
    }
    
    inline bool has_on_negative_side( const point &p) const
    {
        return (a()*p.x + b()*p.y + c()) <  real(0);
    }

    inline bool isSteeper( const line  & l ) const
    {
        if  ( b() == 0 ) 
            return  true;
        if  ( l.b() == 0 )
            return  false;

        return  ( (-a()*l.b()) > (-l.a()*b()) );
    }

    
    bool is_horizontal() const
    {
        
        return a() == real(0) ;
    }
    
    bool is_vertical() const
    {
        return b() == real(0) ;
    }

    bool is_degenerate() const
    {
        return (a() == real(0)) && (b() == real(0)) ;
    }

    bool  isContain( const point  & pnt ) const 
    {
        real  val =  (a() * pnt.x + b() * pnt.y + c());
        return  ( real_eq( val, 0.0 ) );        
    }
};




inline  void  print_spaces( int   count ) 
{
    for   ( int  ind = 0; ind < count; ind++ )
        printf( " " );
}


inline  void   print( int  ind, const line  & l ) 
{
    print_spaces( ind );

    printf( "l:(a:%g, b:%g, c:%g)\n", ( l.a() ), ( l.b() ), 
            ( l.c() ) );
}

void  line_cache_init();


#else   /* __LINE__H */
#error  Header file line.h included twice
#endif  /* __LINE__H */

/*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
 *     
 * line.h - End of File
\*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*/