strech.c

任意给定三维空间的点集

/*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
 * strech.C -
 *     
\*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*/

#include  <stdlib.h>
#include  <stdio.h>
#include  <assert.h>

#include  <stdlib.h>
#include  <stdio.h>
#include  <iostream>
#include  <ctype.h>
#include  <math.h>
#include  <time.h>
#include  <qapplication.h>
#include  <assert.h>
#include  <GL/glut.h>

//#include  "dmalloc.h"

//#include  "scape/vector3d.h"
#include  "_generic.h"
#include  "sarray.h"
#include  "real.h"
//#include  "line.h"
#include  "point.h"
#include  "bbox.h"
#include  "ftree.h"
//#include  "TreeDiam.h"
//#include  "glwindow.h"
#include  "glbox.h"
#include  "misc.h"

//#include  "terrain.h"
#include  "distance.h"
#include  "algorithm.h"
#include  "3dsload.h"
#include  "glwindow_geometry.h"
#include  "sarray.t"
#include  "strech.h"
//#include  "bbox.h"

/*--- Start of Code ---*/


static void   get_line( FILE  * fl, char  * line )
{
    char  * out;

    out = fgets( line, 1024, fl );
    assert( out != NULL );
}

static void  transform_pnt( const H3dsVert   & v,
                            Point3d  & pnt )
{
    pnt.setCoord( 0, v.x );
    pnt.setCoord( 1, v.y );
    pnt.setCoord( 2, v.z );
}



void  get_obj_ratio( H3dsMeshObj  & obj, double  & ratio ) 
{
    ratio = 1;
}


void  get_scene_ratios( H3dsScene * scene )
{
    //int  verts_num;
    //bool   f_first;
    double  ratio, max_ratio;

    max_ratio = 0;

    printf( "\n\n\n\n\nScene objects: %d\n", scene->meshobjs );
    for  ( int  ind = 0; ind < scene->meshobjs; ind++ ) {
        H3dsMeshObj  & obj( scene->meshobjlist[ ind ] );

        printf( "faces: %d\n", obj.faces );
        get_obj_ratio( obj, ratio );
        if  ( ( ind == 0 )  ||  (ratio > max_ratio ) ) {
            max_ratio = ratio;
        }
    }
}


int        comp_vertices( const  Point3d  & a,
                          const  Point3d  & b )
{
    if  ( a.coords[ 0 ] < b.coords[ 0 ] )
        return  -1;
    if  ( a.coords[ 0 ] > b.coords[ 0 ] )
        return  1;

    if  ( a.coords[ 1 ] < b.coords[ 1 ] )
        return  -1;
    if  ( a.coords[ 1 ] > b.coords[ 1 ] )
        return  1;
    
    if  ( a.coords[ 2 ] < b.coords[ 2 ] )
        return  -1;
    if  ( a.coords[ 2 ] > b.coords[ 2 ] )
        return  1;

    //printf( "bogi\n" );

    return  0;   
}


void   get_all_vertices( H3dsScene   * scene, 
                         ArrPoint3d  & arr_out )
{
    ArrPoint3d  arr;
    Point3d pnt;

    arr.init( 100, 324 );    
    arr_out.init( 100, 366 );

    for  ( int  ind = 0; ind < scene->meshobjs; ind++ ) {
        H3dsMeshObj  & obj( scene->meshobjlist[ ind ] );


        for  ( int  ind = 0; ind < obj.faces; ind++ ) {
            H3dsFace   & face( obj.facelist[ ind ] );

            H3dsVert   & v0( obj.vertlist[ face.p0 ] );
            H3dsVert   & v1( obj.vertlist[ face.p1 ] );
            H3dsVert   & v2( obj.vertlist[ face.p2 ] );
            
            transform_pnt( v0, pnt );                    
            arr.pushx( pnt );

            transform_pnt( v1, pnt );                    
            arr.pushx( pnt );

            transform_pnt( v2, pnt );                    
            arr.pushx( pnt );
        }
    }        

    arr.quickSort( comp_vertices );
    for  ( int  ind = 0; ind < (arr.size() - 1); ind++ ) {
        if   ( ind == 0 ) {
            arr_out.pushx( arr[ ind ] );
            continue;
        }
        if  ( ! arr[ ind ].isEqual( arr[ ind - 1 ] ) ) 
            arr_out.pushx( arr[ ind ] );
    }

    for  ( int  ind = 0; ind < (arr_out.size() - 1); ind++ ) {
        if  ( arr_out[ ind ].isEqual( arr_out[ ind + 1 ] ) ) {
            printf( "Vertices are equal!\nBIG ERROR\n" );
            printf( "ind : %d\n", ind );
            arr_out[ ind ].print();
            arr_out[ ind + 1 ].print();

            exit( -1 );
        } 
    }
    printf( "Array size: %d\n", (int)arr.size() );
}


class  SceneInfo
{
public:
    ArrPoint3d  * p_arr;
    H3dsScene  * scene;
    int  size;
    BBox3d  bbox;
    
    double  * nn_dist;
    double  * longest_edge;   
    double  snap_dist, global_vertex_ratio;

    void  fill_info( H3dsScene   * _scene );
    int   get_pnt_ind( const Point3d  & a );
    void  register_longest_edge( Point3d  & a,
                                 Point3d  & b );
    void  register_longest_edge_all();
    void  computeRatio();

    void  register_len( double  dist, double  & len_min,
                        double  & len_max ) ;
    void  get_longest_shortest_edge( double  & len_min, 
                               double  & len_max );

};


inline void  SceneInfo::register_len( double  dist, double  & len_min,
                                      double  & len_max ) 
{
    if  ( dist < snap_dist ) 
        return;
    if  ( dist < len_min )
        len_min = dist;
    if  ( dist > len_max )
        len_max = dist;
}


void  SceneInfo::get_longest_shortest_edge( double  & len_min, 
                                            double  & len_max )
{
    //double  max_len, len;
    bool   f_first;

    f_first = true;
    for  ( int  ind = 0; ind < scene->meshobjs; ind++ ) {
        H3dsMeshObj  & obj( scene->meshobjlist[ ind ] );

        for  ( int  ind = 0; ind < obj.faces; ind++ ) {
            H3dsFace   & face( obj.facelist[ ind ] );

            H3dsVert   & v0( obj.vertlist[ face.p0 ] );
            H3dsVert   & v1( obj.vertlist[ face.p1 ] );
            H3dsVert   & v2( obj.vertlist[ face.p2 ] );

            Point3d  pnts[ 3 ];
            
            transform_pnt( v0, pnts[ 0 ] );
            transform_pnt( v1, pnts[ 1 ] );
            transform_pnt( v2, pnts[ 2 ] );
            
            if  ( f_first ) {
                len_min = len_max = pnt_distance( pnts[ 0 ], pnts[ 1 ] );
                f_first = false;
            }
            register_len( pnt_distance( pnts[ 0 ], pnts[ 1 ] ),
                          len_min, len_max );
            register_len( pnt_distance( pnts[ 1 ], pnts[ 2 ] ),
                          len_min, len_max );
            register_len( pnt_distance( pnts[ 2 ], pnts[ 0 ] ),
                          len_min, len_max );
        }
    }    
}



int    SceneInfo::get_pnt_ind( const Point3d  & a )
{
    int  pos;
    
    pos = p_arr->binarySearch( a, comp_vertices );
    while  ( pos > 0 ) {
        if  ( pnt_distance( (*p_arr)[ pos ],
                            (*p_arr)[ pos - 1 ] ) < snap_dist ) 
            pos--;
        else
            break;
    }

    return  pos;
}


void   SceneInfo::register_longest_edge( Point3d  & a,
                                         Point3d  & b )
{
    int  a_ind, b_ind;

    a_ind = get_pnt_ind( a );
    b_ind = get_pnt_ind( b );

    assert( a_ind >= 0 );
    if  ( b_ind < 0 ) {
        b.print();
        printf( "\n");
        for  ( int  ind = 0; ind < size; ind++ ) {
            printf( "\n%d  ", ind );
            (*p_arr)[ ind ].print();
        }

    }

    assert( b_ind >= 0 );

    double  len = pnt_distance( a, b );
    if  ( len > longest_edge[ a_ind ] )
        longest_edge[ a_ind ] = len;
    if  ( len > longest_edge[ b_ind ] )
        longest_edge[ b_ind ] = len;
}


void   SceneInfo::register_longest_edge_all()
{
    for  ( int  ind = 0; ind < scene->meshobjs; ind++ ) {
        H3dsMeshObj  & obj( scene->meshobjlist[ ind ] );

        for  ( int  ind = 0; ind < obj.faces; ind++ ) {
            H3dsFace   & face( obj.facelist[ ind ] );

            H3dsVert   & v0( obj.vertlist[ face.p0 ] );
            H3dsVert   & v1( obj.vertlist[ face.p1 ] );
            H3dsVert   & v2( obj.vertlist[ face.p2 ] );

            Point3d  pnts[ 3 ];
            
            transform_pnt( v0, pnts[ 0 ] );
            transform_pnt( v1, pnts[ 1 ] );
            transform_pnt( v2, pnts[ 2 ] );

            register_longest_edge( pnts[ 0 ], pnts[ 1 ] );
            register_longest_edge( pnts[ 1 ], pnts[ 2 ] );
            register_longest_edge( pnts[ 2 ], pnts[ 0 ] );
        }
    }

}


void   SceneInfo::fill_info( H3dsScene   * _scene )
{
    scene = _scene;

    snap_dist = 0.0001;
    p_arr = new  ArrPoint3d;

    get_all_vertices( scene, *(p_arr) );    
    size = (int)p_arr->size();

    // Compute bounding box
    bbox.init();
    for  ( int  ind = 0; ind < size; ind++ ) 
        bbox.bound( (*p_arr)[ ind ] );
    snap_dist = bbox.get_diam() / 1000000.0;

    printf( "snap_dist: %g\n", snap_dist );

    nn_dist = (double *)malloc( sizeof( double ) * size );
    assert( nn_dist != NULL );
    longest_edge = (double *)malloc( sizeof( double ) * size );
    assert( longest_edge != NULL );

    for  ( int  ind = 0; ind < size; ind++ ) {
        nn_dist[ ind ] = 10e50;
        longest_edge[ ind ] = 0;
    }
    
    printf( "computing NN... would take a while.\n" );
    for  ( int  ind = 0; ind < size; ind++ ) {
        if  ( ( ind & 0x7 ) == 0 )
            printf( "ind: %d/%d\n", ind, size );
        for  ( int  jnd = 0; jnd < size; jnd++ ) {
            if  ( ind == jnd ) 
                continue;
            double  len = pnt_distance( (*p_arr)[ ind ], (*p_arr)[ jnd ] );
            if  ( len < snap_dist ) 
                continue;
            if  ( len < nn_dist[ ind ] )
                nn_dist[ ind ] = len;
        }
    }

    //printf( "Computing longest edges for every vertex.\n" );
    register_longest_edge_all();
}


void   SceneInfo::computeRatio()
{
    double  ratio, max_ratio;

    max_ratio = 0;

    for  ( int  ind = 0; ind < size; ind++ ) {
        ratio = longest_edge[ ind ] / nn_dist[ ind ];
        if  ( ind == 0 ) {
            max_ratio = ratio;
            continue;
        }
        if  ( longest_edge[ ind ] == 0 ) 
            continue;
        //printf( "%05d: %10g ", ind, ratio );
        //(*p_arr)[ ind ].print();
        if  ( ratio > 300.0 ) {
            printf( "  XXXXX" );
            printf( "\n" );
            printf( "isEqual? %d\n", (int)(*p_arr)[ ind ].isEqual(  
                                       (*p_arr)[ ind - 1 ] ) );
            printf( "\t%g\n", pnt_distance( (*p_arr)[ ind ], 
                                        (*p_arr)[ ind - 1 ] ) );
            printf( "isEqual? %d\n", (int)(*p_arr)[ ind ].isEqual(  
                                       (*p_arr)[ ind + 1 ] ) );
            printf( "\t%g\n", pnt_distance( (*p_arr)[ ind ], 
                                        (*p_arr)[ ind + 1 ] ) );
            printf( "Longest_edge: %g\n", longest_edge[ ind ] );
            printf( "nn dist     : %g\n", nn_dist[ ind ] );

            printf( "\n\n\n\n" );
        }
        //else
        //    printf( "\n" );
        if  ( ratio > max_ratio )
            max_ratio = ratio;
    }

    printf( "Max ratio: %g\n", max_ratio );
    global_vertex_ratio = max_ratio;
}


SceneInfo   * CompSceneInfo( H3dsScene   * scene )
{
    SceneInfo  * info;