Added C++ implementation of terrafit, which is about 6x faster!!!

Had to move parts of Terra into a library and into its own namespace, due to conflicts with SimGear definitions.

src/Prep/ArrayFit/Makefile.am

@@ -12,3 +12,4 @@ demo_SOURCES = demo.cxx
 demo_LDADD = $(support_LIBS) -lz
 
 INCLUDES = -I$(top_srcdir)/src -I$(top_srcdir)/src/Lib
+

src/Prep/Terra/Array.h

@@ -11,6 +11,7 @@
 #define MIN(a,b) (((a)>(b))?(b):(a))
 #endif
 
+namespace Terra {
 
 template<class T>
 class array {
@@ -121,4 +122,6 @@ inline T& array2<T>::ref(int i, int j)
     return data[j*w + i];
 }
 
+}; // namespace Terra
+
 #endif

src/Prep/Terra/Geom.h

@@ -11,6 +11,8 @@
 #include <assert.h>
 #endif
 
+namespace Terra {
+
 typedef double real;
 #define EPS 1e-6
 #define EPS2 (EPS*EPS)
@@ -20,10 +22,14 @@ typedef int boolean;
 enum Axis {X, Y, Z, W};
 enum Side {Left=-1, On=0, Right=1};
 
+}; // namespace Terra
+
 #include <math.h>
 #include "Vec2.h"
 #include "Vec3.h"
 
+namespace Terra {
+
 #ifndef NULL
 #define NULL 0
 #endif
@@ -174,5 +180,6 @@ inline ostream& operator<<(ostream &out, const Line& l)
     return out << "Line(a=" << l.a << " b=" << l.b << " c=" << l.c << ")";
 }
 
+}; // namespace Terra
 
 #endif

src/Prep/Terra/GreedyInsert.cc

@@ -3,11 +3,14 @@
 #include "GreedyInsert.h"
 
 #include "Mask.h"
-extern ImportMask *MASK;
 
 using std::cerr;
 using std::endl;
 
+namespace Terra {
+
+extern ImportMask *MASK;
+
 
 void TrackedTriangle::update(Subdivision& s)
 {
@@ -288,3 +291,5 @@ real GreedySubdivision::eval(int x,int y)
 
     return z_plane(x,y);
 }
+
+}; // namespace Terra

src/Prep/Terra/GreedyInsert.h

@@ -5,6 +5,8 @@
 #include "Subdivision.h"
 #include "Map.h"
 
+namespace Terra {
+
 class TrackedTriangle : public Triangle
 {
     //
@@ -89,4 +91,6 @@ public:
 #define DATA_POINT_IGNORED 2
 #define DATA_VALUE_UNKNOWN 3
 
+}; // namespace Terra
+
 #endif

src/Prep/Terra/Heap.cc

@@ -5,6 +5,7 @@
 using std::cerr;
 using std::endl;
 
+namespace Terra {
 
 void Heap::swap(int i,int j)
 {
@@ -121,3 +122,6 @@ heap_node *Heap::kill(int i)
 
     return &ref(size);
 }
+
+}; // namespace Terra
+

src/Prep/Terra/Heap.h

@@ -4,6 +4,8 @@
 #include "Geom.h"
 #include "Array.h"
 
+namespace Terra {
+
 #define NOT_IN_HEAP -47
 
 //
@@ -54,6 +56,6 @@ public:
     heap_node *kill(int i);
 };
 
-
+}; // namespace Terra
 
 #endif

src/Prep/Terra/Makefile.am

@@ -1,23 +1,28 @@
 EXTRA_DIST = Makefile.gcc-2.95 Makefile.orig
 
 # don't build xterra by default
-# bin_PROGRAMS = terra xterra
+#bin_PROGRAMS = terra xterra
 bin_PROGRAMS = terra
+noinst_LIBRARIES = libTerra.a
+
+libTerra_a_SOURCES = \
+	Array.h Geom.h \
+        GreedyInsert.cc GreedyInsert.h \
+        Heap.cc Heap.h \
+	Map.cc Map.h \
+        Mask.cc Mask.h \
+        Quadedge.cc Quadedge.h \
+        Subdivision.cc Subdivision.h \
+        Vec2.h Vec3.h
 
 terra_SOURCES = \
-	Array.h Geom.h GreedyInsert.cc GreedyInsert.h Heap.cc Heap.h \
+	cmdline.cc greedy.cc output.cc terra.cc terra.h version.h
-	Map.cc Map.h Mask.cc \
-	Mask.h Quadedge.cc Quadedge.h Subdivision.cc Subdivision.h Vec2.h \
-	Vec3.h cmdline.cc greedy.cc \
-	output.cc terra.cc terra.h version.h
 
-terra_LDADD = $(opengl_LIBS)
+terra_LDADD = libTerra.a
 
-# xterra_SOURCES = \
+
-#	Array.h Geom.h GreedyInsert.cc GreedyInsert.h Heap.cc Heap.h \
+#xterra_SOURCES = \
-#	Map.cc Map.h Mask.cc \
+#        cmdline.cc glHacks.cc glHacks.h greedy.cc gui.cc gui.h \
-#	Mask.h Quadedge.cc Quadedge.h Subdivision.cc Subdivision.h Vec2.h \
-#	Vec3.h cmdline.cc glHacks.cc glHacks.h greedy.cc gui.cc gui.h \
 #	output.cc terra.h version.h xterra.cc
 #
-# xterra_LDADD = $(opengl_LIBS)
+#xterra_LDADD = $(opengl_LIBS) libTerra.a

src/Prep/Terra/Map.cc

@@ -3,6 +3,8 @@
 #include "Geom.h"
 #include "Map.h"
 
+namespace Terra {
+
 using std::cerr;
 using std::endl;
 
@@ -78,3 +80,5 @@ Map *readPGM(istream& in)
 
     return map;
 }
+
+}; // namespace Terra

src/Prep/Terra/Map.h

@@ -8,6 +8,8 @@ using std::istream;
 
 #include "Geom.h"
 
+namespace Terra {
+
 class Map
 {
 public:
@@ -102,8 +104,6 @@ void DirectMap<T>::textRead(istream& in)
 	}
 }
 
-
+}; // namespace Terra
-
-
 
 #endif

src/Prep/Terra/Mask.cc

@@ -8,6 +8,7 @@
 using std::cerr;
 using std::endl;
 
+namespace Terra {
 
 RealMask *readMask(istream& in)
 {
@@ -61,3 +62,5 @@ RealMask *readMask(istream& in)
 
     return mask;
 }
+
+}; // namespace Terra

src/Prep/Terra/Mask.h

@@ -7,6 +7,8 @@
 
 using std::istream;
 
+namespace Terra {
+
 class ImportMask
 {
 
@@ -50,4 +52,6 @@ inline real& RealMask::ref(int i, int j)
 
 extern RealMask *readMask(istream&);
 
+}; // namespace Terra
+
 #endif

src/Prep/Terra/Quadedge.cc

@@ -6,6 +6,8 @@
 using std::cerr;
 using std::endl;
 
+namespace Terra {
+
 Edge::Edge(const Edge&)
 {
     cerr << "Edge: Edge assignments are forbidden." << endl;
@@ -83,3 +85,6 @@ void splice(Edge *a, Edge *b)
     alpha->next = t3;
     beta->next = t4;
 }
+
+}; // namespace Terra
+

src/Prep/Terra/Quadedge.h

@@ -3,6 +3,8 @@
 
 #include "Geom.h"
 
+namespace Terra {
+
 class Triangle;
 
 class Edge : public Labelled {
@@ -77,4 +79,6 @@ inline ostream& operator<<(ostream& out, const Edge *e)
     return out << "{ " << e->Org() << " ---> " << e->Dest() << " }";
 }
 
+}; // namespace Terra
+
 #endif

src/Prep/Terra/Subdivision.cc

@@ -8,6 +8,7 @@
 using std::cerr;
 using std::endl;
 
+namespace Terra {
 
 Edge *Subdivision::makeEdge(Vec2& org, Vec2& dest)
 {
@@ -430,3 +431,5 @@ void Triangle::update(Subdivision&)
 
 
 }
+
+}; // namespace Terra

src/Prep/Terra/Subdivision.h

@@ -3,6 +3,8 @@
 
 #include "Quadedge.h"
 
+namespace Terra {
+
 class Subdivision;
 
 class Triangle : public Labelled {
@@ -89,5 +91,6 @@ inline ostream& operator<<(ostream& out, Triangle& t)
 	       << t.point3() << ")";
 }
 
+}; // namespace Terra
 
 #endif

src/Prep/Terra/Vec2.h

@@ -6,6 +6,8 @@
 using std::ostream;
 using std::istream;
 
+namespace Terra {
+
 class Vec2 {
 protected:
     real elt[2];
@@ -169,5 +171,6 @@ inline istream& operator>>(istream& in, Vec2& v)
     return in >> c >> v[0] >> v[1] >> c;
 }
 
+}; // namespace Terra
 
 #endif

src/Prep/Terra/Vec3.h

@@ -1,6 +1,8 @@
 #ifndef VEC3_INCLUDED // -*- C++ -*-
 #define VEC3_INCLUDED
 
+namespace Terra {
+
 class Vec3 {
 protected:
     real elt[3];
@@ -174,5 +176,6 @@ inline istream& operator>>(istream& in, Vec3& v)
     return in >> c >> v[0] >> v[1] >> v[2] >> c;
 }
 
+}; // namespace Terra
 
 #endif

src/Prep/Terra/cmdline.cc

@@ -9,6 +9,8 @@ using std::cin;
 using std::cerr;
 using std::endl;
 
+namespace Terra {
+
 GreedySubdivision *mesh;
 Map *DEM;
 
@@ -187,3 +189,6 @@ void process_cmdline(int argc, char **argv)
 	in.close();
     }
 }
+
+}; // namespace Terra
+

src/Prep/Terra/greedy.cc

@@ -3,6 +3,8 @@
 using std::cerr;
 using std::endl;
 
+namespace Terra {
+
 void scripted_preinsertion(istream& script)
 {
     char op[4];
@@ -109,3 +111,5 @@ void display_greedy_insertion(void (*callback)())
 
     announce_goal();
 }
+
+}; // namespace Terra

src/Prep/Terra/gui.cc

@@ -10,6 +10,7 @@
 using std::cout;
 using std::endl;
 
+namespace Terra {
 
 int mesh_view;
 int surf_view;
@@ -461,3 +462,5 @@ void xglutKeepAspect(float width, float height)
     }
 #endif
 }
+
+}; // namespace Terra

src/Prep/Terra/gui.h

@@ -1,10 +1,14 @@
 #ifndef GUI_INCLUDED // -*- C++ -*-
 #define GUI_INCLUDED
 
+namespace Terra {
+
 extern int mesh_view;
 extern int surf_view;
 
 extern void gui_init();
 extern void gui_interact();
 
+}; // namespace Terra
+
 #endif

src/Prep/Terra/output.cc

@@ -13,6 +13,7 @@ SG_USING_STD(ostream);
 SG_USING_STD(ofstream);
 SG_USING_STD(streampos);
 
+namespace Terra {
 
 void generate_output(char *filename, FileFormat format)
 {
@@ -202,3 +203,5 @@ void output_dem(ostream& out)
     out << endl;
     cerr << ">> Done." << endl;
 }
+
+}; // namespace Terra

src/Prep/Terra/terra.cc

@@ -2,9 +2,9 @@
 
 main(int argc, char **argv)
 {
-    process_cmdline(argc, argv);
+    Terra::process_cmdline(argc, argv);
 
-    greedy_insertion();
+    Terra::greedy_insertion();
 
-    generate_output();
+    Terra::generate_output();
 }

src/Prep/Terra/terra.h

@@ -5,6 +5,7 @@
 #include "Map.h"
 #include "Mask.h"
 
+namespace Terra {
 
 extern GreedySubdivision *mesh;
 extern Map *DEM;
@@ -34,5 +35,6 @@ extern void output_dem(ostream&);
 
 extern void process_cmdline(int argc, char **argv);
 
+}; // namespace Terra
 
 #endif

src/Prep/Terra/xterra.cc

@@ -8,11 +8,11 @@
 main(int argc, char **argv)
 {
     glutInit(&argc, argv);
-    process_cmdline(argc, argv);
+    Terra::process_cmdline(argc, argv);
 
 
-    gui_init();
+    Terra::gui_init();
 
 
-    gui_interact();
+    Terra::gui_interact();
 }

src/Prep/TerraFit/Makefile.am

@@ -1,6 +1,15 @@
 bin_SCRIPTS = terrafit.py
+bin_PROGRAMS = terrafit
+
+terrafit_SOURCES = terrafit.cc
+terrafit_LDADD = $(top_builddir)/src/Prep/Terra/libTerra.a \
+                 $(top_builddir)/src/Lib/Array/libArray.a \
+                 -lsgbucket -lsgstructure -lsgmisc -lsgdebug -lz
+
 EXTRA_DIST = terrafit.README terrafit.py.in
 
 terrafit.py: $(srcdir)/terrafit.py.in
 	echo ${bindir}
 	sed -e "s%__bindir__%${bindir}%g" $(srcdir)/terrafit.py.in > terrafit.py
+
+INCLUDES = -I$(top_srcdir)/src -I$(top_srcdir)/src/Lib

src/Prep/TerraFit/terrafit.cc

@@ -0,0 +1,235 @@
+/*
+ * terrafit.cc - Use Terra as a faster ArrayFit
+ * (see http://graphics.cs.uiuc.edu/~garland/software/terra.html)
+ *
+ * Written by Ralf Gerlich, started December 2007
+ * Based on terrafit.py by Norman Vine
+ *      and Terra by Michael Garland
+ *
+ * Copyright (C) 2007  Ralf Gerlich - http://www.custom-scenery.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.
+ */
+#include <iostream>
+#include <string>
+
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <dirent.h>
+#include <errno.h>
+#include <string.h>
+#include <unistd.h>
+#include <zlib.h>
+
+#include <simgear/debug/logstream.hxx>
+#include <simgear/bucket/newbucket.hxx>
+#include <simgear/misc/sg_path.hxx>
+#include <simgear/structure/exception.hxx>
+
+#include <Array/array.hxx>
+
+#include <Prep/Terra/GreedyInsert.h>
+#include <Prep/Terra/Map.h>
+#include <Prep/Terra/Mask.h>
+
+class ArrayMap: public Terra::Map {
+public:
+        ArrayMap(TGArray& array): array(array) {
+                width=array.get_cols();
+                height=array.get_rows();
+                depth=32;
+        }
+        
+        virtual Terra::real eval(int i, int j) {
+                return (Terra::real)array.get_array_elev(i,j);
+        }
+        
+        /* No direct reading of .arr.gz files */
+        virtual void rawRead(istream&) {
+        }
+        virtual void textRead(istream&) {
+        }
+protected:
+        TGArray& array;
+};
+
+static Terra::ImportMask default_mask;
+namespace Terra {
+/* GreedyInsertion requires us to declare a mask, even if we
+ * don't need one...
+ */
+Terra::ImportMask *MASK=&default_mask;
+}; // namespace Terra
+
+Terra::real error_threshold=40.0;
+unsigned int min_points=50;
+unsigned int point_limit=1000;
+
+inline int goal_not_met(Terra::GreedySubdivision* mesh)
+{
+    return
+        ( mesh->maxError() > error_threshold &&
+          mesh->pointCount() < point_limit ) ||
+        mesh->pointCount() < min_points;
+        
+}
+
+static void announce_goal(Terra::GreedySubdivision* mesh)
+{
+    SG_LOG(SG_GENERAL, SG_INFO, "Goal conditions met:");
+    SG_LOG(SG_GENERAL, SG_INFO, "     error=" << mesh->maxError() << " [thresh="<< error_threshold << "]");
+    SG_LOG(SG_GENERAL, SG_INFO, "     points=" << mesh->pointCount() << " [limit=" << point_limit << "]");
+}
+
+void greedy_insertion(Terra::GreedySubdivision* mesh)
+{
+
+    while( goal_not_met(mesh) )
+    {
+	if( !mesh->greedyInsert() )
+	    break;
+    }
+
+    announce_goal(mesh);
+}
+
+void fit_file(const std::string& path) {
+        SG_LOG(SG_GENERAL, SG_INFO,"Working on file '" << path << "'");
+        
+        std::string infile,outfile;
+        if (path.compare(path.size()-7,7,".arr.gz")==0) {
+                infile=path.substr(0,path.size()-7);
+        } else if (path.compare(path.size()-4,4,".arr")==0) {
+                infile=path.substr(0,path.size()-4);
+        } else {
+                /* actually should not happen */
+                SG_LOG(SG_GENERAL, SG_ALERT, "unknown input file extension = " << path);
+                throw sg_exception("unknown input file extension!");
+        }
+        
+        outfile=infile+".fit.gz";
+        
+        struct stat statbuf;
+        if (stat(path.c_str(),&statbuf)!=0) {
+                SG_LOG(SG_GENERAL, SG_ALERT ,"ERROR: Unable to stat '" << path << "':" << strerror(errno));
+                return;
+        }
+        
+        time_t src_mtime=statbuf.st_mtime;
+        
+        if (stat(outfile.c_str(),&statbuf)==0 && statbuf.st_mtime>src_mtime) {
+                SG_LOG(SG_GENERAL, SG_INFO ,"Skipping " << outfile << ", source " << path << " is older");
+                return;
+        }
+        
+        SGBucket bucket(0,0); // dummy bucket
+        TGArray inarray(infile);
+        inarray.parse(bucket);
+        
+        ArrayMap *DEM=new ArrayMap(inarray);
+        
+        Terra::GreedySubdivision *mesh;
+        
+        mesh=new Terra::GreedySubdivision(DEM);
+        
+        greedy_insertion(mesh);
+        
+        gzFile fp;
+        if ( (fp = gzopen( outfile.c_str(), "wb9" )) == NULL ) {
+                SG_LOG(SG_GENERAL, SG_ALERT, "ERROR: opening " << outfile << " for writing!");
+                return;
+        }
+        
+        for (int x=0;x<DEM->width;x++) {
+                for (int y=0;y<DEM->height;y++) {
+                        if (mesh->is_used(x,y) != DATA_POINT_USED)
+                                continue;
+                        double vx,vy,vz;
+                        vx=(inarray.get_originx()+x*inarray.get_col_step())/3600.0;
+                        vx=(inarray.get_originy()+y*inarray.get_row_step())/3600.0;
+                        vz=DEM->eval(x,y);
+                        gzprintf(fp,"%+03.8f %+02.8f %0.2f\n",vx,vy,vz);
+                }
+        }
+        
+        gzclose(fp);
+}
+
+void walk_path(const std::string& path) {
+        struct stat statbuf;
+        if (stat(path.c_str(),&statbuf)!=0) {
+                SG_LOG(SG_GENERAL, SG_ALERT ,"ERROR: Unable to stat '" << path << "':" << strerror(errno));
+                return;
+        }
+        
+        if (path.compare(path.size()-7,7,".arr.gz")==0 || path.compare(path.size()-4,4,".arr")==0) {
+                fit_file(path);
+        } else if (S_ISDIR(statbuf.st_mode)) {
+                DIR* dir=opendir(path.c_str());
+                if (!dir) {
+                        SG_LOG(SG_GENERAL, SG_ALERT ,"ERROR: Unable to open directory '" << path << "':" << strerror(errno));
+                        return;
+                }
+                
+                struct dirent* dirent;
+                
+                while ((dirent=readdir(dir))) {
+                        if (!strcmp(dirent->d_name,".") || !strcmp(dirent->d_name,"..")) {
+                                continue; // skip . and ..
+                        }
+                        SGPath subpath(path);
+                        subpath.append(dirent->d_name);
+                        walk_path(subpath.str());
+                }
+                
+                closedir(dir);
+        }
+}
+
+void usage(char* progname, char* msg) {
+        if (msg!=NULL)
+                SG_LOG(SG_GENERAL,SG_ALERT, msg);
+        SG_LOG(SG_GENERAL,SG_INFO, "Usage: " << progname << " [options] <file | path to walk>");
+        SG_LOG(SG_GENERAL,SG_INFO, "\t -h | --help ");
+        SG_LOG(SG_GENERAL,SG_INFO, "\t -m | --minnodes 50");
+        SG_LOG(SG_GENERAL,SG_INFO, "\t -x | --maxnodes 1000");
+        SG_LOG(SG_GENERAL,SG_INFO, "\t -e | --maxerror 40");
+        SG_LOG(SG_GENERAL,SG_INFO, "\t -v | --version");
+        SG_LOG(SG_GENERAL,SG_INFO, "");
+        SG_LOG(SG_GENERAL,SG_INFO, "Algorithm will produce at least <minnodes> fitted nodes, but no");
+        SG_LOG(SG_GENERAL,SG_INFO, "more than <maxnodes>.  Within that range, the algorithm will stop");
+        SG_LOG(SG_GENERAL,SG_INFO, "if the maximum elevation error for any remaining point");
+        SG_LOG(SG_GENERAL,SG_INFO, "drops below <maxerror> meters.");
+        SG_LOG(SG_GENERAL,SG_INFO, "");
+        SG_LOG(SG_GENERAL,SG_INFO, "Increasing the maxnodes value and/or decreasing maxerror");
+        SG_LOG(SG_GENERAL,SG_INFO, "will produce a better surface approximation.");
+        SG_LOG(SG_GENERAL,SG_INFO, "");
+        SG_LOG(SG_GENERAL,SG_INFO, "The input file must be a .arr.gz file such as that produced");
+        SG_LOG(SG_GENERAL,SG_INFO, "by demchop or hgtchop utils.");
+        SG_LOG(SG_GENERAL,SG_INFO, "");
+        SG_LOG(SG_GENERAL,SG_INFO, "**** NOTE ****:");
+        SG_LOG(SG_GENERAL,SG_INFO, "If a directory is input all .arr.gz files in directory will be");
+        SG_LOG(SG_GENERAL,SG_INFO, "processed recursively.");
+        SG_LOG(SG_GENERAL,SG_INFO, "");
+        SG_LOG(SG_GENERAL,SG_INFO, "The output file(s) is/are called .fit.gz and is simply a list of");
+        SG_LOG(SG_GENERAL,SG_INFO, "from the resulting fitted surface nodes.  The user of the");
+        SG_LOG(SG_GENERAL,SG_INFO, ".fit.gz file will need to retriangulate the surface.");
+}
+
+int main(int argc, char** argv) {
+        sglog().setLogLevels( SG_ALL, SG_DEBUG );
+        walk_path(argv[1]); // FIXME: parse options
+}
+