diff --git a/configure.ac b/configure.ac index c3f535b0..e31c6454 100644 --- a/configure.ac +++ b/configure.ac @@ -50,6 +50,11 @@ esac AC_SUBST(AR) AC_SUBST(ARFLAGS) +dnl Add gts/glib includes, this will probably need to be made more +dnl flexible in the future. +CFLAGS="$CFLAGS `gts-config --cflags`" +CXXFLAGS="$CXXFLAGS `gts-config --cflags`" + dnl Specify if we want logging (testing build) or not (release build) # set logging default value # with_logging=yes @@ -95,7 +100,7 @@ null_LIBS="$LIBS" AC_CHECK_LIB(m, cos) -base_LIBS="$LIBS" +base_LIBS="$LIBS `gts-config --libs`" AC_CHECK_LIB(socket, socket) AC_CHECK_LIB(X11, XCreateWindow) diff --git a/src/Lib/Array/Makefile.am b/src/Lib/Array/Makefile.am index 9bed153e..9506e2e7 100644 --- a/src/Lib/Array/Makefile.am +++ b/src/Lib/Array/Makefile.am @@ -2,7 +2,7 @@ noinst_LIBRARIES = libArray.a libArray_a_SOURCES = array.cxx array.hxx -noinst_PROGRAMS = testarray +noinst_PROGRAMS = testarray testgts testarray_SOURCES = testarray.cxx @@ -10,4 +10,10 @@ testarray_LDADD = \ $(top_builddir)/src/Lib/Array/libArray.a \ -lsgbucket -lsgmath -lsgmisc -lsgdebug -lsgxml -lz +testgts_SOURCES = testgts.cxx + +testgts_LDADD = \ + $(top_builddir)/src/Lib/Array/libArray.a \ + -lsgbucket -lsgmath -lsgmisc -lsgdebug -lsgxml -lz + INCLUDES = -I$(top_srcdir)/src diff --git a/src/Lib/Array/testgts.cxx b/src/Lib/Array/testgts.cxx new file mode 100644 index 00000000..8a3eb118 --- /dev/null +++ b/src/Lib/Array/testgts.cxx @@ -0,0 +1,309 @@ +// Loads a .arr file (chopped intermediate form of DEM) and leverages +// portions of gts to impliment the terrain simplification algorithm +// in Michael Garlands paper located here: +// +// http://graphics.cs.uiuc.edu/~garland/software/terra.html +// +// Essentially start with two triangles forming the bounding surface. +// Then add the point that has the greatest error. Retriangulate. +// Recalcuate errors for each remaining point, add the one with the +// greatest error. +// +// Outputs to a file that can be visualized with gtsview (available +// from http://gts.sf.net) + +#include + +#include + +#include "array.hxx" + +SG_USING_STD(cout); +SG_USING_STD(endl); + +static void pick_first_face( GtsFace *f, GtsFace **first ) { + if ( *first == NULL ) + *first = f; +} + + +static GtsPoint *global_pt; + +static void find_enclosing_face( GtsFace *f, GtsFace **first ) { + if ( gts_point_is_in_triangle( global_pt, (GtsTriangle *)f ) != GTS_OUT ) { + *first = f; + } +} + + +// if p lies inside plane (in terms of x,y position) return the +// distance from the point to the triangle in the z direction. +// Otherwise return 0. +double calc_error( GtsTriangle *t, GtsPoint *p ) { + if ( gts_point_is_in_triangle( p, t ) == GTS_OUT ) { + // point outside triangle, bail + return 0; + } + + double a, b, c, d; + gts_triangle_normal( t, &a, &b, &c ); + GtsVertex *v1, *v2, *v3; + gts_triangle_vertices( t, &v1, &v2, &v3 ); + GtsPoint *v = (GtsPoint *)v1; + d = a * v->x + b * v->y + c * v->z; + + // cout << "p = " << Point3D( p->x, p->y, p->z ) << endl; + // cout << "coeff = " << Point3D( a, b, c ) << endl; + if ( c < 0.00000001 ) { + cout << "Really small C coefficient" << endl; + exit(-1); + } + + double e = ( d - a * p->x - b * p->y ) / c; + + return fabs( e - p->z ); +} + + +int main( int argc, char **argv ) { + bool verbose = false; + double error_threshold = 20; + + if ( argc != 2 ) { + cout << "Usage: " << argv[0] << " work_dir" << endl; + exit(-1); + } + + string work_dir = argv[1]; + + double lon, lat; + lon = -146.248360; lat = 61.133950; // PAVD (Valdez, AK) + lon = -110.664244; lat = 33.352890; // P13 + lon = -122.374843; lat = 37.619002; // KSFO + + SGBucket b( lon, lat ); + string base = b.gen_base_path(); + string path = work_dir + "/" + base; + + string arrayfile = path + "/" + b.gen_index_str() + ".arr"; + cout << "arrayfile = " << arrayfile << endl; + + TGArray a(arrayfile); + a.parse( b ); + + // Old fit + cout << "Old fit(200) = " << a.fit(200) << endl; + cout << "Old fit(100) = " << a.fit(100) << endl; + cout << "Old fit(50) = " << a.fit(50) << endl; + cout << "Old fit(25) = " << a.fit(25) << endl; + + // Test libgts (gnu triangulation library functions) + + // Load the DEM data and make a list of points + point_list pending; + pending.clear(); + + double x, y, z; + double basex = a.get_originx(); + double basey = a.get_originy(); + double dx = a.get_col_step(); + double dy = a.get_row_step(); + + for ( int i = 0; i < a.get_cols(); ++i ) { + for ( int j = 0; j < a.get_rows(); ++j ) { + if ( (i == 0 && j == 0) || + (i == a.get_cols() - 1 && j == 0 ) || + (i == a.get_cols() - 1 && j == a.get_rows() - 1 ) || + (i == 0 && j == a.get_rows() - 1 ) ) + { + // skip corners since they will be added seperately + } else { + x = basex + i * dx; + y = basey + j * dy; + z = a.get_point( i, j ); + pending.push_back( Point3D(x, y, z) ); + } + } + } + + // Make an (empty) surface + + // Make the corner vertices (enclosing exactly the DEM coverage area) + x = basex; + y = basey; + z = a.interpolate_altitude( x, y ); + cout << "adding = " << Point3D( x, y, z) << endl; + GtsVertex *v1 = gts_vertex_new( gts_vertex_class(), x, y, z ); + + x = basex + dx * (a.get_cols() - 1); + y = basey; + z = a.interpolate_altitude( x, y ); + cout << "adding = " << Point3D( x, y, z) << endl; + GtsVertex *v2 = gts_vertex_new( gts_vertex_class(), x, y, z ); + + x = basex + dx * (a.get_cols() - 1); + y = basey + dy * (a.get_rows() - 1); + z = a.interpolate_altitude( x, y ); + cout << "adding = " << Point3D( x, y, z) << endl; + GtsVertex *v3 = gts_vertex_new( gts_vertex_class(), x, y, z ); + + x = basex; + y = basey + dy * (a.get_rows() - 1); + z = a.interpolate_altitude( x, y ); + cout << "adding = " << Point3D( x, y, z) << endl; + GtsVertex *v4 = gts_vertex_new( gts_vertex_class(), x, y, z ); + + GSList *list = NULL; + list = g_slist_prepend( list, v1 ); + list = g_slist_prepend( list, v2 ); + list = g_slist_prepend( list, v3 ); + list = g_slist_prepend( list, v4 ); + + // make a triangle the completely encloses the 4 corners of our + // DEM + GtsTriangle *t = gts_triangle_enclosing( gts_triangle_class(), + list, 2.0 ); + + // Make the (empty) surface + GtsSurface *surface = gts_surface_new( gts_surface_class(), + gts_face_class(), + gts_edge_class(), + gts_vertex_class() ); + + // add the enclosing surface + gts_surface_add_face( surface, gts_face_new(gts_face_class(), + t->e1, t->e2, t->e3) ); + + // Add the four corners + GtsVertex *result; + result = gts_delaunay_add_vertex( surface, v1, NULL ); + result = gts_delaunay_add_vertex( surface, v2, NULL ); + result = gts_delaunay_add_vertex( surface, v3, NULL ); + result = gts_delaunay_add_vertex( surface, v4, NULL ); + + gts_surface_print_stats( surface, stdout ); + + // add points incrementally from the pending list + + bool done = false; + int count = 4; + GtsPoint *p = gts_point_new( gts_point_class(), 0, 0, 0 ); + global_pt = gts_point_new( gts_point_class(), 0, 0, 0 ); + + while ( !done ) { + // iterate through all the surface faces + + if ( verbose ) { + gts_surface_print_stats( surface, stdout ); + } + cout << "points left = " << pending.size() + << " points added = " << count << endl; + GtsFace *first = NULL; + GtsFace *guess = NULL; + GtsFace *found = NULL; + gts_surface_foreach_face( surface, (GtsFunc)pick_first_face, &first ); + + double max_error = 0; + point_list_iterator mark = NULL; + + // iterate through all remaining points + point_list_iterator current = pending.begin(); + const_point_list_iterator last = pending.end(); + for ( ; current != last; ++current ) { + // cout << *current << endl; + gts_point_set( p, current->x(), current->y(), + current->z() ); + gts_point_set( global_pt, current->x(), current->y(), + current->z() ); + + guess = gts_point_locate( p, surface, guess ); + // gts_surface_foreach_face( surface, (GtsFunc)find_enclosing_face, + // &guess ); + + double error = calc_error( (GtsTriangle *)guess, p ); + if ( error > max_error ) { + max_error = error; + mark = current; + found = guess; + } + } + + if ( max_error > error_threshold ) { + cout << "adding " << *mark << " (" + << max_error << ")" << endl; + GtsVertex *v = gts_vertex_new( gts_vertex_class(), + mark->x(), + mark->y(), + mark->z() ); + GtsVertex *result = gts_delaunay_add_vertex( surface, v, guess ); + if ( result != NULL ) { + cout << " error adding vertex! " << *mark << endl; + } else { + ++count; + } + pending.erase( mark ); + + // GtsFace *f = gts_delaunay_check( surface ); + // if ( f == NULL ) { + // cout << "valid delauney triangulation" << endl; + // } else { + // cout << "NOT VALID DELAUNEY TRIANGULATION" << endl; + // } + } else { + done = true; + } + + FILE *fp = fopen( "surface.gts", "w" ); + gts_surface_write( surface, fp ); + fclose(fp); + + cout << endl; + } + + FILE *fp = fopen( "surface.gts", "w" ); + gts_surface_write( surface, fp ); + fclose(fp); + +#if 0 + + // Make the corner vertices + GtsVertex *v1 = gts_vertex_new( gts_vertex_class(), 0, 0, 0 ); + GtsVertex *v2 = gts_vertex_new( gts_vertex_class(), 10, 0, 5 ); + GtsVertex *v3 = gts_vertex_new( gts_vertex_class(), 10, 10, 10 ); + GtsVertex *v4 = gts_vertex_new( gts_vertex_class(), 0, 10, 5 ); + + // Make the 5 edges + GtsEdge *e1 = gts_edge_new( gts_edge_class(), v1, v2 ); + GtsEdge *e2 = gts_edge_new( gts_edge_class(), v2, v3 ); + GtsEdge *e3 = gts_edge_new( gts_edge_class(), v3, v4 ); + GtsEdge *e4 = gts_edge_new( gts_edge_class(), v4, v1 ); + GtsEdge *e5 = gts_edge_new( gts_edge_class(), v2, v4 ); + + // Make the two faces + GtsFace *f1 = gts_face_new( gts_face_class(), e1, e5, e4 ); + GtsFace *f2 = gts_face_new( gts_face_class(), e2, e3, e5 ); + + // Make the (empty) surface + GtsSurface *surface = gts_surface_new( gts_surface_class(), + gts_face_class(), + gts_edge_class(), + gts_vertex_class() ); + + // Add the two faces + gts_surface_add_face( surface, f1 ); + gts_surface_add_face( surface, f2 ); + + // Add some vertices + for ( int j = 0; j <= 10; ++j ) { + for ( int i = 0; i <= 10; ++i ) { + GtsVertex *v = gts_vertex_new( gts_vertex_class(), i, j, (i - j) ); + GtsVertex *result = gts_delaunay_add_vertex (surface, v, NULL); + if ( result != NULL ) { + cout << " error adding vertex! " << i << " " << j << endl; + } + } + } +#endif + + return 0; +}