From cde416c9f3355b4f4bb8cc0db9be5fb4740b4e9d Mon Sep 17 00:00:00 2001 From: curt Date: Sun, 16 Mar 2003 21:52:10 +0000 Subject: [PATCH] Initial revision of a tool that uses Michael Garland's algorithm to fit a simpler TIN surface to a regular array of height values. The algorithm is described here: http://graphics.cs.uiuc.edu/~garland/software/terra.html --- src/Prep/ArrayFit/Makefile.am | 9 + src/Prep/ArrayFit/arrayfit.cxx | 413 +++++++++++++++++++++++++++++++++ 2 files changed, 422 insertions(+) create mode 100644 src/Prep/ArrayFit/Makefile.am create mode 100644 src/Prep/ArrayFit/arrayfit.cxx diff --git a/src/Prep/ArrayFit/Makefile.am b/src/Prep/ArrayFit/Makefile.am new file mode 100644 index 00000000..71565e54 --- /dev/null +++ b/src/Prep/ArrayFit/Makefile.am @@ -0,0 +1,9 @@ +noinst_PROGRAMS = arrayfit + +arrayfit_SOURCES = arrayfit.cxx + +arrayfit_LDADD = \ + $(top_builddir)/src/Lib/Array/libArray.a \ + -lsgbucket -lsgmath -lsgmisc -lsgdebug -lsgxml -lz + +INCLUDES = -I$(top_srcdir)/src -I$(top_srcdir)/src/Lib diff --git a/src/Prep/ArrayFit/arrayfit.cxx b/src/Prep/ArrayFit/arrayfit.cxx new file mode 100644 index 00000000..6633d2cc --- /dev/null +++ b/src/Prep/ArrayFit/arrayfit.cxx @@ -0,0 +1,413 @@ +// arrayfit.cxx +// +// 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. Lather, rinse, repeat. +// +// The resulting fitted set of nodes is written out to a file so the +// main tile builder can later load these nodes and incorporate them +// into the tile surface. +// +// Written by Curtis Olson, started March 2003. +// +// Copyright (C) 2003 Curtis L. Olson - curt@flightgear.org +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. +// +// $Id$ + + +#include +#include + +#include + +#include + +SG_USING_STD(cout); +SG_USING_STD(endl); + + +// transform point to lat/lon degree coordinates and append to list +static void add_point( point_list &list, Point3D p ) { + Point3D tp( p.x() / 3600.0, p.y() / 3600.0, p.z() ); + list.push_back( tp ); +} + + +static void pick_first_face( GtsFace *f, GtsFace **first ) { + if ( *first == NULL ) + *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 ); +} + + +static void usage( char *prog ) { + cout << "Usage: " << prog << " [ --options ]" << endl; + cout << "\t--input=file.arr" << endl; + cout << "\t--minnodes=50" << endl; + cout << "\t--maxnodes=600" << endl; + cout << "\t--maxerror=50" << endl; + cout << "\t--verbose" << endl; + cout << endl; + cout << "Algorithm will produce at least 50 fitted nodes, but no" << endl; + cout << "more than 600. Within that range, the algorithm will stop"<< endl; + cout << "if the maximum elevation error for any remaining point" << endl; + cout << "drops below 50 meters." << endl; + cout << endl; + cout << "Increasing the maxnodes value and/or decreasing maxerror" << endl; + cout << "will produce a better surface approximation." << endl; + cout << endl; + cout << "The input file must be a .arr file such as that produced" << endl; + cout << "by demchop or hgtchop utils." << endl; + cout << endl; + cout << "The output file is called .fit and is simply a list of" << endl; + cout << "from the resulting fitted surface nodes. The user of the" << endl; + cout << ".fit file will need to retriangulate the surface." << endl; + exit(-1); +} + + +int main( int argc, char **argv ) { + // option defaults + SGPath infile; + int min_nodes = 50; + int max_nodes = 600; + double error_threshold = 50.0; + bool verbose = false; + + // Parse command line arguments + for ( int i = 1; i < argc; ++i ) { + string arg = argv[i]; + + if ( arg.find("--input=") == 0 ) { + infile.set( arg.substr(8) ); + } else if ( arg.find("--minnodes=") == 0 ) { + min_nodes = atoi( arg.substr(11).c_str() ); + } else if ( arg.find("--maxnodes=") == 0 ) { + max_nodes = atoi( arg.substr(11).c_str() ); + } else if ( arg.find("--maxerror=") == 0 ) { + error_threshold = atof( arg.substr(11).c_str() ); + } else if ( arg.find("--verbose") == 0 ) { + verbose = true; + } else { + usage( argv[0] ); + } + } + + if ( ! infile.str().length() ) { + usage( argv[0] ); + } + + SGPath outfile = infile.base(); + cout << outfile.str() << endl; + while ( outfile.extension().length() ) { + outfile = outfile.base(); + cout << outfile.str() << endl; + } + outfile.concat( ".fit.gz" ); + + cout << "Input file = " << infile.str() << endl; + cout << "Outfile file = " << outfile.str() << endl; + cout << "Minimum nodes = " << min_nodes << endl; + cout << "Maximum nodes = " << max_nodes << endl; + cout << "Error Threshold = " << error_threshold << endl; + cout << endl; + + SGBucket b(0.0, 0.0); // build a dummy bucket + TGArray a( infile.str() ); + a.parse( b ); + + // 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) ); + } + } + } + + // Create the empty fitted list (this is the list we are working + // so hard to create.) :-) + point_list fitted; + fitted.clear(); + double_list errors; + errors.clear(); + + // 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; + add_point( fitted, Point3D( x, y, z) ); + errors.push_back( 13000.0 ); + 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; + add_point( fitted, Point3D( x, y, z) ); + errors.push_back( 13000.0 ); + 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; + add_point( fitted, Point3D( x, y, z) ); + errors.push_back( 13000.0 ); + 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; + add_point( fitted, Point3D( x, y, z) ); + errors.push_back( 13000.0 ); + 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 + // terrain data + 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 ); + + if ( verbose ) { + gts_surface_print_stats( surface, stdout ); + } + + // add the remaining points incrementally (from the pending list) + + bool done = false; + int count = 4; + GtsPoint *p = 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 + << " fitted list size = " << fitted.size() << 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() ); + + guess = gts_point_locate( p, surface, guess ); + + double error = calc_error( (GtsTriangle *)guess, p ); + if ( error > max_error ) { + max_error = error; + mark = current; + found = guess; + } + } + + // check stop conditions + if ( (max_error < error_threshold && (int)fitted.size() >= min_nodes) || + (int)fitted.size() >= max_nodes ) + { + // we are done + done = true; + } else { + // add the next point and keep going + cout << "adding " << *mark << " (" + << max_error << ")" << endl; + add_point( fitted, *mark ); + errors.push_back( max_error ); + 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; + // } + } + + if ( verbose ) { + FILE *fp = fopen( "surface.gts", "w" ); + gts_surface_write( surface, fp ); + fclose(fp); + } + + if ( verbose ) { + cout << endl; + } + } + + if ( verbose ) { + FILE *fp = fopen( "surface.gts", "w" ); + gts_surface_write( surface, fp ); + fclose(fp); + } + + gzFile gzfp; + if ( (gzfp = gzopen( outfile.c_str(), "wb9" )) == NULL ) { + cout << "ERROR: cannot open " << outfile.str() << " for writing!" + << endl; + exit(-1); + } + + gzprintf( gzfp, "%d\n", fitted.size() ); + for ( unsigned int i = 0; i < fitted.size(); ++i ) { + gzprintf( gzfp, "%.15f %.15f %.2f %.1f\n", + fitted[i].x(), fitted[i].y(), fitted[i].z(), errors[i] ); + } + gzclose( gzfp ); + +#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; +}