gdiam.c

任意给定三维空间的点集

    //printf( "area: %g\n", area );
    return  ( ( area < (ldouble)0.0 )? -1: 
              ( (area > (ldouble)0.0)? 1 : 0 ) );
}

inline  bool    Left(  const point2d  & a, 
                         const point2d  & b,
                         const point2d  & c )
{
    return  AreaSign( a, b, c ) > 0;
}

inline  bool    Left_colinear(  const point2d  & a, 
                                const point2d  & b,
                                const point2d  & c )
{
    return  AreaSign( a, b, c ) >= 0;
}


struct CompareByAngle {  
public:    
    point2d  base;

    bool operator()(const point2d_ptr  & a, const point2d_ptr   & b ) {
	int  sgn;
	gdiam_real  len1, len2;

        if  ( a->equal( *b ) )
            return  false;
        assert( a != NULL );
        assert( b != NULL );
        if  ( a->equal( base ) ) {
            assert( false );
            return  true;
        }
        if  ( b->equal( base ) ) {
            assert( false );
            return  false;
        }

        //printf( "before: %p %p...\n", a, b );
        //fflush( stdout );
        /*
        if  ( a->equal( base ) )
            return  true;
        if  ( b->equal( base ) )
            return  false;
        */
        //printf( "a: (%g, %g)\n", a->x, a->y );
        //printf( "b: (%g, %g)\n", b->x, b->y );
        //fflush( stdout );

	sgn = AreaSign( base, *a, *b );
        if  ( sgn != 0 )
            return  (sgn > 0);
        len1 = base.dist( *a );
        len2 = base.dist( *b );
        
        //printf( "bogi!\n" );
        //fflush( stdout );

        return (len1 > len2);
    }
};


point2d_ptr  get_min_point( vec_point_2d  & in,
                            int  & index ) 
{
    point2d_ptr  min_pnt = in[ 0 ];

    index = 0;

    for  ( int  ind = 0; ind < (int)in.size(); ind++ ) { 
        if  ( in[ ind ]->y < min_pnt->y ) {
            min_pnt = in[ ind ];
            index = ind;
        } else
            if  ( ( in[ ind ]->y == min_pnt->y ) 
                  &&  ( in[ ind ]->x < min_pnt->x ) ) {
                min_pnt = in[ ind ];
                index = ind;
            }
    }

    return  min_pnt;        
}


const void  dump( vec_point_2d   & vec ) 
{
    for  ( int  ind = 0; ind < (int)vec.size(); ind++ ) {
        printf( "-- %11d (%-11g, %-11g)\n",                
                ind, vec[ ind ]->x,
                vec[ ind ]->y );
    }
    printf( "\n\n" );
}


static void  print_pnt( point2d_ptr  pnt ) 
{
    printf( "(%g, %g)\n", pnt->x, pnt->y );
}


void   verify_convex_hull( vec_point_2d  & in, vec_point_2d  & ch )
{
    ldouble  area;

    //dump( ch );
    //fflush( stdout );
    //fflush( stderr );
    for  ( int  ind = 0; ind < (int)( ch.size() - 2 ); ind++ )
        assert( Left( *( ch[ ind ] ), *( ch[ ind + 1 ] ),
                      *( ch[ ind + 2 ] ) ) );
    assert( Left( *( ch[ ch.size() - 2 ] ), *( ch[ ch.size() - 1 ] ),
                  *( ch[ 0 ] ) ) );
    assert( Left( *( ch[ ch.size() - 1 ] ), *( ch[ 0 ] ),
                  *( ch[ 1 ] ) ) );

    for  ( int  ind = 0; ind < (int)( in.size() - 2 ); ind++ ) {
         for  ( int  jnd = 0; jnd < (int)( ch.size() - 1 ); jnd++ ) {
             if  ( ch[ jnd ]->equal_real( *( in[ ind ] ) ) )
                 continue;
             if  ( ch[ jnd + 1 ]->equal( *( in[ ind ] ) ) )
                 continue;
             
             if  ( ! Left_colinear( *( ch[ jnd ] ), *( ch[ jnd + 1 ] ),
                                    *( in[ ind ] ) ) ) {
                 area = Area( *( ch[ jnd ] ), *( ch[ jnd + 1 ] ),
                              *( in[ ind ] ) );
                 if  ( fabs( area ) < 1e-12 )
                     continue;
                 printf( "Failure in progress!\n\n" );
                 print_pnt( ch[ jnd ] );
                 print_pnt( ch[ jnd + 1 ] );
                 print_pnt( in[ ind ] );

                 //dump( ch );
                 printf( "ch[ jnd ]: (%g, %g)\n", 
                         ch[ jnd ]->x, ch[ jnd ]->y );
                 printf( "ch[ jnd + 1 ]: (%g, %g)\n", 
                         ch[ jnd + 1 ]->x, ch[ jnd + 1 ]->y );
                 printf( "ch[ ind ]: (%g, %g)\n", 
                         in[ ind ]->x, in[ ind ]->y );
                 printf( "Area: %g\n", (double)Area( *( ch[ jnd ] ), 
                                             *( ch[ jnd + 1 ] ),
                                             *( in[ ind ] ) ) );
                 printf( "jnd: %d, jnd + 1: %d, ind: %d\n",
                         jnd, jnd+1, ind );
                 fflush( stdout );
                 fflush( stderr );
                 
                 assert( Left( *( ch[ jnd ] ), *( ch[ jnd + 1 ] ),
                               *( in[ ind ] ) ) );
             }
         }
         if  ( ch[ 0 ]->equal( *( in[ ind ] ) ) )
             continue;
         if  ( ch[ ch.size() - 1 ]->equal( *( in[ ind ] ) ) )
             continue;
         assert( Left( *( ch[ ch.size() - 1 ] ), 
                       *( ch[ 0 ] ),
                       *( in[ ind ] ) ) );
    }
    printf( "Convex hull seems to be ok!\n" );
}



static  void   remove_duplicate( vec_point_2d  & in, 
                                 point2d_ptr  ptr,
                                 int  start )
{
    int  dest;

    dest = start;
    while  ( start < (int)in.size() ) {
        if  ( in[ start ]->equal_real( *ptr ) ) {
            start++;
            continue;
        }
        in[ dest++ ] = in[ start++ ];
    }
    while  ( (int)in.size() > dest ) 
        in.pop_back();
}


static  void   remove_consecutive_dup( vec_point_2d  & in )
{
    int  dest, pos;

    if  ( in.size() < 1 ) 
        return;
    
    dest = pos = 0;

    // copy first entry
    in[ dest++ ] = in[ pos++ ];
    while  ( pos < ( (int)in.size() - 1 ) ) {
        if  ( in[ pos - 1 ]->equal_real( *(in[ pos ]) ) ) {
            pos++;
            continue;
        }
        in[ dest++ ] = in[ pos++ ];        
    }
    in[ dest++ ] = in[ pos++ ];

    while  ( (int)in.size() > dest ) 
        in.pop_back();
}


void  convex_hull( vec_point_2d  & in, vec_point_2d  & out )
{
    
    CompareByAngle  comp;
    int  ind, position;
    
    assert( in.size() > 1 );
    
    comp.base = *( get_min_point( in, position ));
    std::swap( in[ 0 ], in[ position ] );
    
    remove_duplicate( in, in[ 0 ], 1 );
    
    /*
    printf( "comp.base: (%g, %g)\n",
            comp.base.x,
            comp.base.y );
    */
    for  ( int  ind = 0; ind < (int)in.size(); ind++ ) {
        assert( in[ ind ] != NULL );
    }
        
    int  size;

    size = in.size();

    /*
    if  ( in.size() == 24 ) {
        dump( in );
        fflush( stdout );
        fflush( stderr );
        }*/

    //printf( "sort( %d, %d, comp )\n", 1, in.size() );
    sort( in.begin() + 1, in.end(), comp );
    remove_consecutive_dup( in );
    
    //dump( in );
    /*
    for  ( int  ind = 0; ind < (int)in.size(); ind++ ) {
        double  x_delta, y_delta;

        x_delta = in[ ind ]->x - comp.base.x;
        y_delta = in[ ind ]->y - comp.base.y;

        printf( "-- %11d (%-11g, %-11g)  atan2: %-11g\n",                
                ind, in[ ind ]->x,
                in[ ind ]->y,
                atan2( y_delta, x_delta ) );
                }*/

    //fflush( stdout );
    //printf( "pushing...\n" );
    //fflush( stdout );
    out.push_back( in[ in.size() - 1 ] );
    out.push_back( in[ 0 ] );

    // perform the graham scan
    ind = 1;
    int  last;

    while ( ind < (int)in.size() ) {
        if  ( out[ out.size() -1 ]->equal_real( *( in[ ind ] ) ) ) {
            ind++;
            continue;
        }
        last = out.size();
        assert( last > 1 );

        if ( Left( *(out[ last - 2 ]),
                   *(out[ last - 1 ]),
                   *(in[ ind ]) ) ) {
            if  ( ! out[ last - 1 ]->equal( *( in[ ind ] ) ) )
                out.push_back( in[ ind ] );
            ind++;
        } else {
            if  ( out.size() < 3 ) {
                dump( out );
                printf( "Input:\n" );
                dump( in );
                printf( "ind: %d\n", ind );
                fflush( stdout );
                assert( out.size() > 2 );
            }

            out.pop_back();
        }
    }

    // we pushed in[ in.size() -1 ] twice in the output
    out.pop_back();

    //verify_convex_hull( in, out );
}


struct bbox_2d_info
{
    gdiam_point_2d_t  vertices[ 4 ];
    gdiam_real  area;

    void  get_dir( int  ind, gdiam_point_2d  out ) {
        out[ 0 ] = vertices[ ind ][ 0 ]
            - vertices[ (ind + 1) % 4 ][ 0 ];
        out[ 1 ] = vertices[ ind ][ 1 ]
            - vertices[ (ind + 1) % 4 ][ 1 ];
    }

    void  dump() {
        printf( "--- bbox 2d ----------------------------\n" );
        for  ( int  ind = 0; ind < 4; ind++ )
            printf( "%d: (%g, %g)\n", ind, vertices[ ind ][0],
                    vertices[ ind ][1] );
        for  ( int  ind = 0; ind < 4; ind++ ) {
            gdiam_point_2d_t  dir;

            get_dir( ind, dir );

            printf( "dir %d: (%g, %g)\n", ind, dir[0],
                    dir[1] );
        }
        printf( "\\----------------------------------\n" );
    }

};

#define  EPS  1e-6
inline  bool   eq_real( gdiam_real  a, gdiam_real  b ) {
    return  ( ( ( b - EPS ) < a )  
              &&  ( a < (b + EPS) ) );
}

class  MinAreaRectangle
{
private:
    vec_point_2d   ch;
    gdiam_real  * angles;

public:
    void  dump_ch() {
        for  ( int  ind = 0; ind < (int)ch.size(); ind++ ) {
            printf( "%d: (%g, %g)\n", ind, ch[ ind ]->x,
                    ch[ ind ]->y );
        }
    }

    void  compute_edge_dir( int  u ) {
        int  u2;
        double  ang; 
        double  x1, y1, x2, y2;

        u2 = (u + 1) % ch.size();

        x1 = ch[ u ]->x;
        y1 = ch[ u ]->y;
        x2 = ch[ u2 ]->x;
        y2 = ch[ u2 ]->y;

        ang = atan2( y2 - y1, x2 - x1 );
        if  ( ang < 0 )
            ang += 2.0 * PI;
        angles[ u ] = ang; 
       
        // floating point nonsence. A huck to solve this...
        if  ( ( u > 0 )  
              &&  ( angles[ u ] < angles[ u - 1 ] ) 
              &&   eq_real( angles[ u ], angles[ u - 1 ] ) )
            angles[ u ] = angles[ u - 1 ];

        /*
        printf( "angles[ %4d ]: %.20f ", u, ang );
        if  ( u > 0 ) 
            printf( "%2d", (angles[ u ] >= angles[ u - 1 ] ) );
        if  ( u > 1 ) {
            printf( " area: %.20f ", (double)Area( *( ch[ u - 2 ] ),
                                                  *( ch[ u - 1 ] ),
                                                  *( ch[ u ] ) ) );
        }
        printf( "\n" );
        */
    }

    /* Finding the first vertex whose edge is over a given angle */
    int  find_vertex( int  start_vertex,
                      double  trg_angle ) {
        double  prev_angle, curr_angle;
        bool  f_found;
        int  ver, prev_vertex;

        f_found = false;

        prev_vertex = start_vertex - 1;
        if  ( prev_vertex < 0 ) 
            prev_vertex = ch.size() - 1;

        prev_angle = angles[ prev_vertex ];
        ver = start_vertex;
        while  ( ! f_found ) {
            curr_angle = angles[ ver ];
            if ( prev_angle <= curr_angle ) 
                f_found = ( ( prev_angle < trg_angle)
                            &&  ( trg_angle <= curr_angle ) ); 
            else 
                f_found = ( ( prev_angle < trg_angle )
                            ||  ( trg_angle <= curr_angle )); 
            
            if ( f_found) 
                break; 
            else 
                prev_angle = curr_angle; 
            
            ver = ( ver + 1 ) % ch.size();
        } 

        return  ver;
    }


    /* Computing the intersection point of two lines, each given by a
   point and slope */
    void  compute_crossing( int  a_ind, double  a_angle,
                            int  b_ind, double  b_angle,
                            gdiam_real  & x, 
                            gdiam_real  & y ) {
        double  a_slope, b_slope, x_a, x_b, y_a, y_b;

        if ( eq_real( a_angle, 0.0 )  ||  eq_real( a_angle, PI ) ) {
            x = ch[ b_ind ]->x;
            y = ch[ a_ind ]->y;
            return;
        }
        if ( eq_real( a_angle, PI/2.0 )  ||  eq_real( a_angle, PI * 1.5 ) ) {
            x = ch[ a_ind ]->x;
            y = ch[ b_ind ]->y;
            return;
        }
        
        a_slope = tan( a_angle ); 
        b_slope = tan( b_angle ); 
        x_a = ch[ a_ind ]->x; 
        y_a = ch[ a_ind ]->y; 
        x_b = ch[ b_ind ]->x; 
        y_b = ch[ b_ind ]->y; 

        double  a_const, b_const, x_out, y_out;

        // Let see:
        //  l_a = y = a_slope * x + a_const ->