From 2a00d8092b2ddb8b64503aa08a89b6f167b99d12 Mon Sep 17 00:00:00 2001 From: Ralf Gerlich Date: Thu, 10 Jan 2008 20:31:18 +0100 Subject: [PATCH] Got rid of the last gts-dependency-leftovers. The gts-checks have been commented out for quite some time already and ArrayFit hasn't been used for years. --- README.gts | 43 ---- configure.ac | 49 ---- src/Lib/Array/Makefile.am | 7 - src/Lib/Array/testgts.cxx | 309 ------------------------ src/Prep/ArrayFit/.cvsignore | 5 - src/Prep/ArrayFit/Makefile.am | 15 -- src/Prep/ArrayFit/arrayfit.cxx | 429 --------------------------------- src/Prep/ArrayFit/demo.cxx | 271 --------------------- src/Prep/Makefile.am | 2 - 9 files changed, 1130 deletions(-) delete mode 100644 README.gts delete mode 100644 src/Lib/Array/testgts.cxx delete mode 100644 src/Prep/ArrayFit/.cvsignore delete mode 100644 src/Prep/ArrayFit/Makefile.am delete mode 100644 src/Prep/ArrayFit/arrayfit.cxx delete mode 100644 src/Prep/ArrayFit/demo.cxx diff --git a/README.gts b/README.gts deleted file mode 100644 index f2c6bd13..00000000 --- a/README.gts +++ /dev/null @@ -1,43 +0,0 @@ -You need the Gnu Triangulated Surfaces library (GTS) installed on your -system to build portions of TerraGear. GTS provides the needed tools -and infrastructure to impliment our terrain simplification approach. -We modeled our approach after the approach outlined in Michael -Garland's paper here: - - http://graphics.cs.uiuc.edu/~garland/software/terra.html - -You can get the latest version of GTS from: - - http://gts.sourceforge.net/ - -More information: - - GTS stands for the GNU Triangulated Surface Library. It is an Open - Source Free Software Library intended to provide a set of useful - functions to deal with 3D surfaces meshed with interconnected - triangles. The source code is available free of charge under the - Free Software LGPL license. - - The code is written entirely in C with an object-oriented approach - based mostly on the design of GTK+. Careful attention is paid to - performance related issues as the initial goal of GTS is to provide - a simple and efficient library to scientists dealing with 3D - computational surface meshes. - - A brief summary of its main features: - - * Simple object-oriented structure giving easy access to - topological properties. - * 2D dynamic Delaunay and constrained Delaunay triangulations. - * Robust geometric predicates (orientation, in circle) using fast - adaptive floating point arithmetic (adapted from the fine work - of Jonathan R. Shewchuk). - * Robust set operations on surfaces (union, intersection, difference). - * Surface refinement and coarsening (multiresolution models). - * Dynamic view-independent continuous level-of-detail. - * Preliminary support for view-dependent level-of-detail. - * Bounding-boxes trees and Kd-trees for efficient point location - and collision/intersection detection. - * Graph operations: traversal, graph partitioning. - * Metric operations (area, volume, curvature ...). - * Triangle strips generation for fast rendering. diff --git a/configure.ac b/configure.ac index 9890e9be..96a6a656 100644 --- a/configure.ac +++ b/configure.ac @@ -49,41 +49,6 @@ 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. -# AC_CHECK_PROG(GLIB, glib-config, yes, no) -# AC_CHECK_PROG(GTS, gts-config, yes, no) -# -# if test "$GLIB" = "no"; then -# echo -# echo "Unable to find glib-config." -# echo -# echo "This program is needed to determine the compiler flags needed for" -# echo "the glib library. Please make sure this linrary is installed and" -# echo "the program is in the search path." -# echo -# echo "Please read README.gts" for more details. -# echo -# exit 1 -# fi -# -# if test "$GTS" = "no"; then -# echo -# echo "Unable to find gts-config." -# echo -# echo "This program is needed to determine the compiler flags needed for" -# echo "the gts library. Please make sure this linrary is installed and" -# echo "the program is in the search path." -# echo -# echo "Please read README.gts" for more details. -# echo -# exit 1 -# -# fi - -# SUPPORT_FLAGS="`gts-config --cflags` `glib-config --cflags`" -# CPPFLAGS="$CPPFLAGS $SUPPORT_FLAGS" - dnl Specify if we want logging (testing build) or not (release build) # set logging default value # with_logging=yes @@ -215,7 +180,6 @@ dnl check for some default libraries AC_SEARCH_LIBS(cos, m) base_LIBS="$LIBS" -# support_LIBS="`gts-config --libs` `glib-config --libs`" AC_SEARCH_LIBS(XCreateWindow, X11) AC_SEARCH_LIBS(XShmCreateImage, Xext) @@ -325,18 +289,6 @@ AC_CHECK_HEADER(newmat/newmat.h) AM_CONDITIONAL(HAVE_NEWMAT, test "x$ac_cv_header_newmat_newmat_h" = "xyes" ) AC_LANG_POP -# AC_CHECK_HEADER(gts.h) -# if test "x$ac_cv_header_gts_h" != "xyes"; then -# echo -# echo "You *must* have the gts library installed on your system to build" -# echo "TerraGear!" -# echo -# echo "Please see README.gts for more details." -# echo -# echo "configure aborted." -# exit -# fi - dnl Check if Generic Polygon Clipping library is installed dnl (from http://www.cs.man.ac.uk/aig/staff/alan/software/) AC_CHECK_HEADERS( gpc.h ) @@ -480,7 +432,6 @@ AC_CONFIG_FILES([ \ src/Lib/TriangleJRS/Makefile \ src/Lib/vpf/Makefile \ src/Prep/Makefile \ - src/Prep/ArrayFit/Makefile \ src/Prep/DemChop/Makefile \ src/Prep/DemInfo/Makefile \ src/Prep/DemRaw2ascii/Makefile \ diff --git a/src/Lib/Array/Makefile.am b/src/Lib/Array/Makefile.am index c1d2de5c..4ba97e7c 100644 --- a/src/Lib/Array/Makefile.am +++ b/src/Lib/Array/Makefile.am @@ -12,11 +12,4 @@ testarray_LDADD = \ -lsgbucket -lsgmath -lsgmisc -lsgdebug -lsgxml \ $(support_LIBS) -lz -# testgts_SOURCES = testgts.cxx - -# testgts_LDADD = \ -# $(top_builddir)/src/Lib/Array/libArray.a \ -# -lsgbucket -lsgmath -lsgmisc -lsgdebug -lsgxml \ -# $(support_LIBS) -lz - INCLUDES = -I$(top_srcdir)/src diff --git a/src/Lib/Array/testgts.cxx b/src/Lib/Array/testgts.cxx deleted file mode 100644 index 373564ac..00000000 --- a/src/Lib/Array/testgts.cxx +++ /dev/null @@ -1,309 +0,0 @@ -// 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. -// -// 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_array_elev( 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.altitude_from_grid( 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.altitude_from_grid( 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.altitude_from_grid( 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.altitude_from_grid( 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; - - // iterate through all remaining points - point_list_iterator current = pending.begin(); - point_list_iterator mark = 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; -} diff --git a/src/Prep/ArrayFit/.cvsignore b/src/Prep/ArrayFit/.cvsignore deleted file mode 100644 index 1b9da297..00000000 --- a/src/Prep/ArrayFit/.cvsignore +++ /dev/null @@ -1,5 +0,0 @@ -.deps -Makefile -Makefile.in -arrayfit -demo diff --git a/src/Prep/ArrayFit/Makefile.am b/src/Prep/ArrayFit/Makefile.am deleted file mode 100644 index f55d40f2..00000000 --- a/src/Prep/ArrayFit/Makefile.am +++ /dev/null @@ -1,15 +0,0 @@ -bin_PROGRAMS = arrayfit demo - -arrayfit_SOURCES = arrayfit.cxx - -arrayfit_LDADD = \ - $(top_builddir)/src/Lib/Array/libArray.a \ - -lsgbucket -lsgmath -lsgmisc -lsgdebug -lsgxml \ - $(support_LIBS) -lz - -demo_SOURCES = demo.cxx - -demo_LDADD = $(support_LIBS) -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 deleted file mode 100644 index 245f6066..00000000 --- a/src/Prep/ArrayFit/arrayfit.cxx +++ /dev/null @@ -1,429 +0,0 @@ -// 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 - http://www.flightgear.org/~curt -// -// 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: arrayfit.cxx,v 1.10 2004-11-19 22:25:51 curt Exp $ - - -#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 triangulate 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 - unsigned int i; - for ( 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] ); - } - - // Strip off .arr or .arr.gz if part of the input file name - string s = infile.str(); - if(s.length() >= 7) { - if(s.substr(s.length() - 7) == ".arr.gz") { - s = s.substr(0, s.length() - 7); - } - } - if(s.length() >= 4) { - if(s.substr(s.length() - 4) == ".arr") { - s = s.substr(0, s.length() - 4); - } - } - infile.set(s); - - SGPath outfile = infile; - 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 ( 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_array_elev( i, j ); - if ( z > -9000 ) { - pending.push_back( Point3D(x, y, z) ); - } else { - // just ignore voids, better just not include - // them, rather than making a stupid guess at a - // value - } - } - } - } - - // 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.altitude_from_grid( 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.altitude_from_grid( 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.altitude_from_grid( 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.altitude_from_grid( 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 && pending.size() > 0 ) { - // 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 = pending.end(); - - // 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 ( 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; -} diff --git a/src/Prep/ArrayFit/demo.cxx b/src/Prep/ArrayFit/demo.cxx deleted file mode 100644 index f3dc4e28..00000000 --- a/src/Prep/ArrayFit/demo.cxx +++ /dev/null @@ -1,271 +0,0 @@ -// demo.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 - http://www.flightgear.org/~curt -// -// 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: demo.cxx,v 1.2 2004-11-19 22:25:51 curt Exp $ - - -#include - -#include - -#include - -SG_USING_STD(cin); -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 ); -} - - -int main( int argc, char **argv ) { - // option defaults - int min_nodes = 50; - int max_nodes = 6000; - double error_threshold = 0.1; - bool verbose = false; - - cout << "Minimum nodes = " << min_nodes << endl; - cout << "Maximum nodes = " << max_nodes << endl; - cout << "Error Threshold = " << error_threshold << endl; - cout << endl; - - // Load the DEM data and make a list of points - point_list pending; - pending.clear(); - - double x, y, z; - - while ( true ) { - cin >> x >> y >> z; - if ( x == 9999 && y == 9999 && z == 9999 ) { - break; - } else { - 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(); - - // Make the corner vertices (enclosing exactly the DEM coverage area) - add_point( fitted, pending[0] ); - GtsVertex *v1 = gts_vertex_new( gts_vertex_class(), - pending[0].x(), - pending[0].y(), - pending[0].z() ); - - add_point( fitted, pending[1] ); - GtsVertex *v2 = gts_vertex_new( gts_vertex_class(), - pending[1].x(), - pending[1].y(), - pending[1].z() ); - - add_point( fitted, pending[2] ); - GtsVertex *v3 = gts_vertex_new( gts_vertex_class(), - pending[2].x(), - pending[2].y(), - pending[2].z() ); - - add_point( fitted, pending[3] ); - GtsVertex *v4 = gts_vertex_new( gts_vertex_class(), - pending[3].x(), - pending[3].y(), - pending[3].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 = pending.end(); - - // 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 ); - 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); - } - - return 0; -} diff --git a/src/Prep/Makefile.am b/src/Prep/Makefile.am index 12b55ce0..8922fa04 100644 --- a/src/Prep/Makefile.am +++ b/src/Prep/Makefile.am @@ -1,5 +1,3 @@ -DISTDIRS = ArrayFit - SUBDIRS = \ DemChop \ DemInfo \