Original version of Michael Garlands terra program version 0.7 (public domain)

src/Prep/Terra/Array.h

@@ -0,0 +1,124 @@
+#ifndef ARRAY_INCLUDED // -*- C++ -*-
+#define ARRAY_INCLUDED
+
+//
+// Array classes
+//
+// Taken from gfxTools.h 1.2
+
+#ifndef MAX
+#define MAX(a,b) (((a)>(b))?(a):(b))
+#define MIN(a,b) (((a)>(b))?(b):(a))
+#endif
+
+
+template<class T>
+class array {
+protected:
+    T *data;
+    int len;
+public:
+    array() { data=NULL; len=0; }
+    array(int l) { init(l); }
+    ~array() { free(); }
+
+
+    inline void init(int l);
+    inline void free();
+    inline void resize(int l);
+
+    inline T& ref(int i);
+    inline T& operator[](int i) { return data[i]; }
+    inline T& operator()(int i) { return ref(i); }
+    inline int length() { return len; }
+    inline int maxLength() { return len; }
+};
+
+template<class T>
+inline void array<T>::init(int l)
+{
+    data = new T[l];
+    len = l;
+}
+
+template<class T>
+inline void array<T>::free()
+{
+    if( data )
+    {
+        delete[] data;
+        data = NULL;
+    }
+}
+
+template<class T>
+inline T& array<T>::ref(int i)
+{
+#ifdef SAFETY
+    assert( data );
+    assert( i>=0 && i<len );
+#endif
+    return data[i];
+}
+
+template<class T>
+inline void array<T>::resize(int l)
+{
+    T *old = data;
+    data = new T[l];
+    data = (T *)memcpy(data,old,MIN(len,l)*sizeof(T));
+    len = l;
+    delete[] old;
+}
+
+
+
+template<class T>
+class array2 {
+protected:
+    T *data;
+    int w, h;
+public:
+    array2() { data=NULL; w=h=0; }
+    array2(int w, int h) { init(w,h); }
+    ~array2() { free(); }
+
+    inline void init(int w, int h);
+    inline void free();
+
+    inline T& ref(int i, int j);
+    inline T& operator()(int i,int j) { return ref(i,j); }
+    inline int width() { return w; }
+    inline int height() { return h; }
+};
+
+template<class T>
+inline void array2<T>::init(int width,int height)
+{
+    w = width;
+    h = height;
+    data = new T[w*h];
+}
+
+template<class T>
+inline void array2<T>::free()
+{
+    if( data )
+    {
+        delete[] data;
+        data = NULL;
+    }
+}
+
+template<class T>
+inline T& array2<T>::ref(int i, int j)
+{
+#ifdef SAFETY
+    assert( data );
+    assert( i>=0 && i<w );
+    assert( j>=0 && j<h );
+#endif
+    return data[j*w + i];
+}
+
+#endif

src/Prep/Terra/Geom.h

@@ -0,0 +1,182 @@
+#ifndef GEOM_INCLUDED // -*- C++ -*-
+#define GEOM_INCLUDED
+
+////////////////////////////////////////////////////////////////////////
+//
+// Define some basic types and values
+//
+////////////////////////////////////////////////////////////////////////
+
+#ifdef SAFETY
+#include <assert.h>
+#endif
+
+typedef double real;
+#define EPS 1e-6
+#define EPS2 (EPS*EPS)
+
+typedef int boolean;
+
+enum Axis {X, Y, Z, W};
+enum Side {Left=-1, On=0, Right=1};
+
+#include <math.h>
+#include "Vec2.h"
+#include "Vec3.h"
+
+#ifndef NULL
+#define NULL 0
+#endif
+
+#ifndef True
+#define True 1
+#define False 0
+#endif
+
+class Labelled {
+public:
+    int token;
+};
+
+
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// Here we define some useful geometric functions
+//
+////////////////////////////////////////////////////////////////////////
+
+//
+// triArea returns TWICE the area of the oriented triangle ABC.
+// The area is positive when ABC is oriented counterclockwise.
+inline real triArea(const Vec2& a, const Vec2& b, const Vec2& c)
+{
+    return (b[X] - a[X])*(c[Y] - a[Y]) - (b[Y] - a[Y])*(c[X] - a[X]);
+}
+
+inline boolean ccw(const Vec2& a, const Vec2& b, const Vec2& c)
+{
+    return triArea(a, b, c) > 0;
+}
+
+inline boolean rightOf(const Vec2& x, const Vec2& org, const Vec2& dest)
+{
+    return ccw(x, dest, org);
+}
+
+inline boolean leftOf(const Vec2& x, const Vec2& org, const Vec2& dest)
+{
+    return ccw(x, org, dest);
+}
+
+// Returns True if the point d is inside the circle defined by the
+// points a, b, c. See Guibas and Stolfi (1985) p.107.
+//
+inline boolean inCircle(const Vec2& a, const Vec2& b, const Vec2& c,
+			const Vec2& d)
+{
+    return (a[0]*a[0] + a[1]*a[1]) * triArea(b, c, d) -
+	   (b[0]*b[0] + b[1]*b[1]) * triArea(a, c, d) +
+ 	   (c[0]*c[0] + c[1]*c[1]) * triArea(a, b, d) -
+	   (d[0]*d[0] + d[1]*d[1]) * triArea(a, b, c) > EPS;
+}
+
+
+
+class Plane {
+public:
+
+    real a, b, c;
+
+    Plane() { }
+    Plane(const Vec3& p, const Vec3& q, const Vec3& r) { init(p,q,r); }
+
+    inline void init(const Vec3& p, const Vec3& q, const Vec3& r);
+
+    real operator()(real x,real y) { return a*x + b*y + c; }
+    real operator()(int x, int y)  { return a*x + b*y + c; }
+};
+
+inline void Plane::init(const Vec3& p, const Vec3& q, const Vec3& r)
+// find the plane z=ax+by+c passing through three points p,q,r
+{
+    // We explicitly declare these (rather than putting them in a
+    // Vector) so that they can be allocated into registers.
+    real ux = q[X]-p[X], uy = q[Y]-p[Y], uz = q[Z]-p[Z];
+    real vx = r[X]-p[X], vy = r[Y]-p[Y], vz = r[Z]-p[Z];
+    real den = ux*vy-uy*vx;
+
+    a = (uz*vy - uy*vz)/den;
+    b = (ux*vz - uz*vx)/den;
+    c = p[Z] - a*p[X] - b*p[Y];
+}
+
+
+
+class Line {
+
+private:
+    real a, b, c;
+
+public:
+    Line(const Vec2& p, const Vec2& q)
+    {
+	Vec2 t = q - p;
+	real l = t.length();
+#ifdef SAFETY
+	assert(l!=0);
+#endif
+	a =   t[Y] / l;
+	b = - t[X] / l;
+	c = -(a*p[X] + b*p[Y]);
+    }
+
+    inline real eval(const Vec2& p) const
+    {
+	return (a*p[X] + b*p[Y] + c);
+    }
+
+    inline Side classify(const Vec2& p) const
+    {
+	real d = eval(p);
+
+	if( d < -EPS )
+	    return Left;
+	else if( d > EPS )
+	    return Right;
+	else
+	    return On;
+    }
+
+    inline Vec2 intersect(const Line& l) const
+    {
+	Vec2 p;
+	intersect(l, p);
+	return p;
+    }
+
+    inline void intersect(const Line& l, Vec2& p) const
+    {
+	real den = a*l.b - b*l.a;
+#ifdef SAFETY
+	assert(den!=0);
+#endif
+	p[X] = (b*l.c - c*l.b)/den;
+	p[Y] = (c*l.a - a*l.c)/den;
+    }
+
+#ifdef IOSTREAMH
+    friend ostream& operator<<(ostream&, const Line&);
+#endif
+};
+
+#ifdef IOSTREAMH
+inline ostream& operator<<(ostream &out, const Line& l)
+{
+    return out << "Line(a=" << l.a << " b=" << l.b << " c=" << l.c << ")";
+}
+#endif
+
+
+#endif

src/Prep/Terra/GreedyInsert.cc

@@ -0,0 +1,269 @@
+#include <iostream.h>
+#include "GreedyInsert.h"
+
+#include "Mask.h"
+extern ImportMask *MASK;
+
+void TrackedTriangle::update(Subdivision& s)
+{
+    GreedySubdivision& gs = (GreedySubdivision&)s;
+    gs.scanTriangle(*this);
+}
+
+
+
+
+
+GreedySubdivision::GreedySubdivision(Map *map)
+{
+    H = map;
+    heap = new Heap(128);
+
+    int w = H->width;
+    int h = H->height;
+
+    is_used.init(w, h);
+    int x,y;
+    for(x=0;x<w;x++)
+	for(y=0;y<h;y++)
+	    is_used(x,y) = DATA_POINT_UNUSED;
+
+
+    initMesh(Vec2(0,0),
+	     Vec2(0, h-1),
+	     Vec2(w-1, h-1),
+	     Vec2(w-1, 0));
+
+    is_used(0, 0)     = DATA_POINT_USED;
+    is_used(0, h-1)   = DATA_POINT_USED;
+    is_used(w-1, h-1) = DATA_POINT_USED;
+    is_used(w-1, 0)   = DATA_POINT_USED;
+
+    count = 4;
+}
+
+
+
+
+
+Triangle *GreedySubdivision::allocFace(Edge *e)
+{
+    Triangle *t = new TrackedTriangle(e);
+
+    heap->insert(t, -1.0);
+
+    return t;
+}
+
+
+
+
+void GreedySubdivision::compute_plane(Plane& plane,
+				      Triangle& T,
+				      Map& map)
+{
+    const Vec2& p1 = T.point1();
+    const Vec2& p2 = T.point2();
+    const Vec2& p3 = T.point3();
+
+    Vec3 v1(p1, map(p1[X], p1[Y]));
+    Vec3 v2(p2, map(p2[X], p2[Y]));
+    Vec3 v3(p3, map(p3[X], p3[Y]));
+
+    plane.init(v1, v2, v3);
+}
+
+///////////////////////////
+//
+// This is indeed an ugly hack.
+// It should be replaced
+//
+static int vec2_y_compar(const void *a,const void *b)
+{
+    Vec2 &p1=*(Vec2 *)a,
+	 &p2=*(Vec2 *)b;
+
+    return (p1[Y]==p2[Y]) ? 0 : (p1[Y] < p2[Y] ? -1 : 1);
+}
+
+static void order_triangle_points(Vec2 *by_y,
+				  const Vec2& p1,
+				  const Vec2& p2,
+				  const Vec2& p3)
+{
+    by_y[0] = p1;
+    by_y[1] = p2;
+    by_y[2] = p3;
+
+    qsort(by_y,3,sizeof(Vec2),vec2_y_compar);
+}
+
+
+
+void GreedySubdivision::scan_triangle_line(Plane& plane,
+					   int y,
+					   real x1, real x2,
+					   Candidate& candidate)
+{
+    int startx = (int)ceil(MIN(x1,x2));
+    int endx   = (int)floor(MAX(x1,x2));
+
+    if( startx > endx ) return;
+
+    real z0 = plane(startx, y);
+    real dz = plane.a;
+    real z, diff;
+
+    for(int x=startx;x<=endx;x++)
+    {
+	if( !is_used(x,y) )
+	{
+	    z = H->eval(x,y);
+	    diff = fabs(z - z0);
+
+	    candidate.consider(x, y, MASK->apply(x, y, diff));
+	}
+
+	z0 += dz;
+    }
+}
+
+
+void GreedySubdivision::scanTriangle(TrackedTriangle& T)
+{
+    Plane z_plane;
+    compute_plane(z_plane, T, *H);
+
+    Vec2 by_y[3];
+    order_triangle_points(by_y,T.point1(),T.point2(),T.point3());
+    Vec2& v0 = by_y[0];
+    Vec2& v1 = by_y[1];
+    Vec2& v2 = by_y[2];
+
+
+    int y;
+    int starty, endy;
+    Candidate candidate;
+
+    real dx1 = (v1[X] - v0[X]) / (v1[Y] - v0[Y]);
+    real dx2 = (v2[X] - v0[X]) / (v2[Y] - v0[Y]);
+
+    real x1 = v0[X];
+    real x2 = v0[X];
+
+    starty = (int)v0[Y];
+    endy   = (int)v1[Y];
+    for(y=starty;y<endy;y++) {
+	scan_triangle_line(z_plane, y, x1, x2, candidate);
+
+        x1 += dx1;
+        x2 += dx2;
+    }
+    
+    /////////////////////////////
+
+    dx1 = (v2[X] - v1[X]) / (v2[Y] - v1[Y]);
+    x1 = v1[X];
+
+    starty = (int)v1[Y];
+    endy   = (int)v2[Y];
+    for(y=starty;y<=endy;y++) {
+	scan_triangle_line(z_plane, y, x1, x2, candidate);
+
+        x1 += dx1;
+        x2 += dx2;
+    }
+
+    /////////////////////////////////
+    //
+    // We have now found the appropriate candidate point.
+    //
+    if( candidate.import < 1e-4 )
+    {
+	if( T.token != NOT_IN_HEAP )
+	    heap->kill(T.token);
+
+#ifdef SAFETY
+	T.setCandidate(-69, -69, 0.0);
+#endif
+    }
+    else
+    {
+	assert( !is_used(candidate.x, candidate.y) );
+
+	T.setCandidate(candidate.x, candidate.y, candidate.import);
+	if( T.token == NOT_IN_HEAP )
+	    heap->insert(&T, candidate.import);
+	else
+	    heap->update(&T, candidate.import);
+    }
+}
+
+Edge *GreedySubdivision::select(int sx, int sy, Triangle *t)
+{
+    if( is_used(sx, sy) )
+    {
+	cerr << "   WARNING: Tried to reinsert point: " << sx<<" "<<sy<<endl;
+	return NULL;
+    }
+
+    is_used(sx,sy) = DATA_POINT_USED;
+    count++;
+    return insert(Vec2(sx,sy), t);
+}
+
+
+int GreedySubdivision::greedyInsert()
+{
+    heap_node *node = heap->extract();
+
+    if( !node ) return False;
+
+    TrackedTriangle &T = *(TrackedTriangle *)node->obj;
+    int sx, sy;
+    T.getCandidate(&sx, &sy);
+
+    select(sx, sy, &T);
+
+    return True;
+}
+
+real GreedySubdivision::maxError()
+{
+    heap_node *node = heap->top();
+
+    if( !node )
+	return 0.0;
+
+    return node->import;
+}
+
+real GreedySubdivision::rmsError()
+{
+    real err = 0.0;
+    int width = H->width;
+    int height = H->height;
+
+
+
+    for(int i=0; i<width; i++)
+	for(int j=0; j<height; j++)
+	{
+	    real diff = eval(i, j) - H->eval(i, j);
+	    err += diff * diff;
+	}
+
+    return sqrt(err / (width * height));
+}
+
+
+real GreedySubdivision::eval(int x,int y) 
+{
+    Vec2 p(x,y);
+    Triangle *T = locate(p)->Lface();
+
+    Plane z_plane;
+    compute_plane(z_plane, *T, *H);
+
+    return z_plane(x,y);
+}

src/Prep/Terra/GreedyInsert.h

@@ -0,0 +1,92 @@
+#ifndef GREEDYINSERT_INCLUDED // -*- C++ -*-
+#define GREEDYINSERT_INCLUDED
+
+#include "Heap.h"
+#include "Subdivision.h"
+#include "Map.h"
+
+class TrackedTriangle : public Triangle
+{
+    //
+    // candidate position
+    int sx, sy;
+
+
+public:
+    TrackedTriangle(Edge *e, int t=NOT_IN_HEAP)
+	: Triangle(e, t)
+    {
+    }
+
+    void update(Subdivision&);
+
+
+    void setCandidate(int x,int y, real) { sx=x; sy=y; }
+    void getCandidate(int *x, int *y) { *x=sx; *y=sy; }
+};
+
+
+class Candidate
+{
+public:
+
+    int x, y;
+    real import;
+
+    Candidate() { import = -HUGE; }
+
+    void consider(int sx, int sy, real i)
+    {
+	if( i > import )
+	{
+	    x = sx;
+	    y = sy;
+	    import = i;
+	}
+    }
+};
+
+
+class GreedySubdivision : public Subdivision
+{
+    Heap *heap;
+    unsigned int count;
+
+protected:
+
+    Map *H;
+
+    Triangle *allocFace(Edge *e);
+
+    void compute_plane(Plane&, Triangle&, Map&);
+
+    void scan_triangle_line(Plane& plane,
+			    int y, real x1, real x2,
+			    Candidate& candidate);
+
+public:
+    GreedySubdivision(Map *map);
+
+    array2<char> is_used;
+
+    Edge *select(int sx, int sy, Triangle *t=NULL);
+
+    Map& getData() { return *H; }
+
+    void scanTriangle(TrackedTriangle& t);
+    int greedyInsert();
+
+    unsigned int pointCount() { return count; }
+    real maxError();
+    real rmsError();
+    real eval(int x,int y);
+};
+
+//
+// These are the possible values of is_used(x,y):
+#define DATA_POINT_UNUSED  0
+#define DATA_POINT_USED    1
+#define DATA_POINT_IGNORED 2
+#define DATA_VALUE_UNKNOWN 3
+
+#endif

src/Prep/Terra/Heap.cc

@@ -0,0 +1,120 @@
+#include <assert.h>
+#include <iostream.h>
+#include "Heap.h"
+
+
+void Heap::swap(int i,int j)
+{
+    heap_node tmp = ref(i);
+
+    ref(i) = ref(j);
+    ref(j) = tmp;
+
+    ref(i).obj->token = i;
+    ref(j).obj->token = j;
+}
+
+void Heap::upheap(int i)
+{
+    if( i==0 ) return;
+
+    if( ref(i).import > ref(parent(i)).import ) {
+        swap(i,parent(i));
+        upheap(parent(i));
+    }
+}
+
+void Heap::downheap(int i)
+{
+    if (i>=size) return;        // perhaps just extracted the last
+
+    int largest = i,
+        l = left(i),
+        r = right(i);
+
+    if( l<size && ref(l).import > ref(largest).import ) largest = l;
+    if( r<size && ref(r).import > ref(largest).import ) largest = r;
+
+    if( largest != i ) {
+        swap(i,largest);
+        downheap(largest);
+    }
+}
+
+
+
+void Heap::insert(Labelled *t,real v)
+{
+    if( size == maxLength() )
+    {
+	cerr << "NOTE: Growing heap from " << size << " to " << 2*size << endl;
+	resize(2*size);
+    }
+
+    int i = size++;
+
+    ref(i).obj = t;
+    ref(i).import = v;
+
+    ref(i).obj->token = i;
+
+    upheap(i);
+}
+
+void Heap::update(Labelled *t,real v)
+{
+    int i = t->token;
+
+    if( i >= size )
+    {
+	cerr << "WARNING: Attempting to update past end of heap!" << endl;
+	return;
+    }
+    else if( i == NOT_IN_HEAP )
+    {
+	cerr << "WARNING: Attempting to update object not in heap!" << endl;
+	return;
+    }
+
+    real old=ref(i).import;
+    ref(i).import = v;
+
+    if( v<old )
+        downheap(i);
+    else
+        upheap(i);
+}
+
+
+
+heap_node *Heap::extract()
+{
+    if( size<1 ) return 0;
+
+    swap(0,size-1);
+    size--;
+
+    downheap(0);
+
+    ref(size).obj->token = NOT_IN_HEAP;
+
+    return &ref(size);
+}
+
+heap_node *Heap::kill(int i)
+{
+    if( i>=size )
+	cerr << "WARNING: Attempt to delete invalid heap node." << endl;
+
+    swap(i, size-1);
+    size--;
+    ref(size).obj->token = NOT_IN_HEAP;
+
+    if( ref(i).import < ref(size).import )
+	downheap(i);
+    else
+	upheap(i);
+
+
+    return &ref(size);
+}

src/Prep/Terra/Heap.h

@@ -0,0 +1,59 @@
+#ifndef HEAP_INCLUDED // -*- C++ -*-
+#define HEAP_INCLUDED
+
+#include "Geom.h"
+#include "Array.h"
+
+#define NOT_IN_HEAP -47
+
+//
+//
+// This file extracted from ~/anim/lab/mlab
+//
+//
+
+class heap_node {
+public:
+    real import;
+    Labelled *obj;
+
+    heap_node() { obj=NULL; import=0.0; }
+    heap_node(Labelled *t, double i=0.0) { obj=t; import=i; }
+    heap_node(const heap_node& h) { import=h.import; obj=h.obj; }
+};
+
+
+
+class Heap : public array<heap_node> {
+
+    //
+    // The actual size of the heap.  array::length()
+    // simply returns the amount of allocated space
+    int size;
+
+    void swap(int i, int j);
+
+    int parent(int i) { return (i-1)/2; }
+    int left(int i) { return 2*i+1; }
+    int right(int i) { return 2*i+2; }
+
+    void upheap(int i);
+    void downheap(int i);
+
+public:
+
+    Heap() { size=0; }
+    Heap(int s) : array<heap_node>(s) { size=0; }
+
+
+    void insert(Labelled *, real);
+    void update(Labelled *, real);
+
+    heap_node *extract();
+    heap_node *top() { return size<1 ? (heap_node *)NULL : &ref(0); }
+    heap_node *kill(int i);
+};
+
+
+
+#endif

src/Prep/Terra/Makefile.orig

@@ -0,0 +1,87 @@
+#################################################################
+#
+# Configuration variables.
+# You should change these to fit your system.
+#
+
+CC = cc
+C++ = CC
+
+# For compiling on SGI's with the pre-5.3 (ie. cfront-based) compiler:
+# add '-ptr/tmp/terra_ptrepository' to OPTFLAGS
+# add '-pte.cc' to LFLAGS
+
+OPTFLAGS = -g -mips2 
+# OPTFLAGS = -O2 -mips2
+
+GUI_LIBS = -lglut -lGLU -lGL -lXmu -lX11
+LIBS = -lmalloc -lmx
+
+#
+# This defines the location of the GLUT libraries
+#
+ANIM = /afs/cs/project/anim/garland
+GLUT_FLAGS = 
+GLUT_INCDIR = $(ANIM)/include
+GLUT_LIBDIR = $(ANIM)/lib
+
+#
+# Include any other search directories you need here
+#
+INCDIR = -I$(GLUT_INCDIR)
+LIBDIR = -L$(GLUT_LIBDIR)
+
+#
+# These are the flags for compilation (CFLAGS) and linking (LFLAGS) 
+#
+CFLAGS = $(INCDIR) $(OPTFLAGS) -DSAFETY
+LFLAGS = $(LIBDIR) $(OPTFLAGS)
+
+
+#################################################################
+#
+# Rules for building the Terra programs.
+# You should not need to change anything here.
+#
+
+.SUFFIXES: .cc
+.cc.o:
+	$(C++) $(CFLAGS) -c $<
+
+.C.o:
+	$(C++) $(CFLAGS) -c $<
+
+
+BASE_SRCS = Quadedge.cc Subdivision.cc Map.cc Mask.cc cmdline.cc \
+            GreedyInsert.cc Heap.cc greedy.cc output.cc
+GUI_SRCS  = glHacks.cc gui.cc xterra.cc
+
+TERRA_SRCS = terra.cc $(BASE_SRCS)
+XTERRA_SRCS = $(GUI_SRCS) $(BASE_SRCS)
+
+TERRA_OBJS = $(TERRA_SRCS:.cc=.o)
+XTERRA_OBJS = $(XTERRA_SRCS:.cc=.o)
+
+
+TARGETS = terra xterra
+
+all: $(TARGETS)
+
+terra: $(TERRA_OBJS)
+	$(C++) $(LFLAGS) -o terra $(TERRA_OBJS) $(LIBS)
+
+xterra: $(XTERRA_OBJS)
+	$(C++) $(LFLAGS) -o xterra $(XTERRA_OBJS) $(GUI_LIBS) $(LIBS)
+
+clean :
+	/bin/rm -f $(XTERRA_OBJS) $(TERRA_OBJS) $(TARGETS)
+	/bin/rm -f -r ii_files ptrepository
+	find . -name '*~'    -print -exec rm -f {} \;
+	find . -name 'core' -print -exec rm -f {} \;
+
+depend :
+	touch Makefile.depend
+	makedepend -fMakefile.depend $(INCDIR) -I/usr/include/CC $(BASE_SRCS) $(GUI_SRCS)
+	/bin/rm -f Makefile.depend.bak
+
+sinclude Makefile.depend

src/Prep/Terra/Map.cc

@@ -0,0 +1,77 @@
+#include <math.h>
+
+#include "Geom.h"
+#include "Map.h"
+
+void Map::findLimits()
+{
+    min = HUGE;
+    max = -HUGE;
+
+    for(int i=0;i<width;i++)
+	for(int j=0;j<height;j++)
+	{
+	    real val = eval(i,j);
+
+	    if( val<min ) min = val;
+	    if( val>max ) max = val;
+
+	}
+}
+
+Map *readPGM(istream& in)
+{
+    char magicP, magicNum;
+
+    int width, height, maxval;
+
+    in >> magicP;
+    in >> magicNum;
+    in >> width >> height >> maxval;
+
+    if( magicP != 'P' )
+    {
+	cerr << "readPGM: This is not PGM data." << endl;
+	return NULL;
+    }
+
+    Map *map;
+
+    if( maxval < 256 )
+    {
+	cerr << "readPGM: Allocating a ByteMap to hold data." << endl;
+	map = new ByteMap(width, height);
+    }
+    else if( maxval < 65536 )
+    {
+	cerr << "readPGM: Allocating a ShortMap to hold data." << endl;
+	map = new ShortMap(width, height);
+    }
+    else
+    {
+	cerr << "readPGM: Allocating a WordMap to hold data." << endl;
+	map = new WordMap(width, height);
+    }
+
+    switch( magicNum )
+    {
+    case '2':
+	cerr << "readPGM: Looks like textual PGM data" << endl;
+	map->textRead(in);
+	break;
+
+    case '5':
+	cerr << "readPGM: Looks like binary PGM data" << endl;
+	in.get(magicP); // read the EOL
+	map->rawRead(in);
+	break;
+
+    default:
+	cerr << "readPGM: This is not PGM data." << endl;
+	return NULL;
+    }
+
+    map->findLimits();
+
+    return map;
+}

src/Prep/Terra/Map.h

@@ -0,0 +1,107 @@
+#ifndef MAP_INCLUDED // -*- C++ -*-
+#define MAP_INCLUDED
+
+#include <stdlib.h>
+#include <iostream.h>
+
+#include "Geom.h"
+
+class Map
+{
+public:
+
+    int width;
+    int height;
+    int depth;  // in bits
+
+    real min, max;
+
+    real operator()(int i, int j) { return eval(i,j); }
+    real operator()(real i, real j) { return eval((int)i,(int)j); }
+    real eval(real i, real j) { return eval((int)i,(int)j); }
+    virtual real eval(int i, int j) = 0;
+
+    virtual void rawRead(istream&) = 0;
+    virtual void textRead(istream&) = 0;
+    virtual void *getBlock() { return NULL; }
+    virtual void findLimits();
+};
+
+extern Map *readPGM(istream&);
+
+
+
+template<class T>
+class DirectMap : public Map
+{
+    T *data;
+
+protected:
+    inline T& ref(int i,int j)
+    {
+#ifdef SAFETY
+	assert(i>=0); assert(j>=0); assert(i<width); assert(j<height);
+#endif
+
+	return data[j*width + i];
+    }
+
+public:
+
+    DirectMap(int width, int height);
+
+    real eval(int i, int j) { return (real)ref(i,j); }
+    void *getBlock() { return data; }
+
+    void rawRead(istream&);
+    void textRead(istream&);
+};
+
+typedef DirectMap<unsigned char>  ByteMap;
+typedef DirectMap<unsigned short> ShortMap;
+typedef DirectMap<unsigned int>   WordMap;
+typedef DirectMap<real>           RealMap;
+
+
+
+template<class T>
+DirectMap<T>::DirectMap(int w, int h)
+{
+    width = w;
+    height = h;
+    depth = sizeof(T) << 3;
+
+    data = (T *)calloc(w*h, sizeof(T));
+}
+
+template<class T>
+void DirectMap<T>::rawRead(istream& in)
+{
+    char *loc = (char *)data;
+    int target = width*height*sizeof(T);
+
+    while( target>0 )
+    {
+	in.read(loc, target);
+	target -= in.gcount();
+	loc += in.gcount();
+    }
+}
+
+template<class T>
+void DirectMap<T>::textRead(istream& in)
+{
+    for(int j=0;j<height;j++)
+	for(int i=0;i<width;i++)
+	{
+	    real val;
+	    in >> val;
+	    ref(i,j) = (T)val;
+	}
+}
+
+
+
+
+
+#endif

src/Prep/Terra/Mask.cc

@@ -0,0 +1,60 @@
+#include <math.h>
+#include <stdlib.h>
+#include <iostream.h>
+
+#include "Geom.h"
+#include "Mask.h"
+
+
+RealMask *readMask(istream& in)
+{
+    char magicP, magicNum;
+
+    int width, height, maxval;
+
+    in >> magicP >> magicNum;
+    in >> width >> height >> maxval;
+
+    if( magicP != 'P' )
+    {
+	cerr << "readMask: This is not PGM data." << endl;
+	return NULL;
+    }
+
+    RealMask *mask = new RealMask(width, height);
+
+    if( magicNum == '2' )
+    {
+	for(int j=0; j<height; j++)
+	    for(int i=0; i<width; i++)
+	    {
+		real val;
+		in >> val;
+		mask->ref(i, j) = val;
+	    }
+    }
+    else if( magicNum == '5' )
+    {
+	for(int j=0; j<height; j++)
+	    for(int i=0; i<width; i++)
+	    {
+		unsigned char val;
+		in >> val;
+		mask->ref(i, j) = (real)val;
+	    }
+    }
+    else
+    {
+	cerr << "readMask: This is not PGM data." << endl;
+	return NULL;
+    }
+
+
+    real max = (real)maxval;
+
+    for(int i=0; i<width; i++)
+	for(int j=0; j<height; j++)
+	    mask->ref(i,j) /= max;
+
+    return mask;
+}

src/Prep/Terra/Mask.h

@@ -0,0 +1,47 @@
+#ifndef MASK_INCLUDED // -*- C++ -*-
+#define MASK_INCLUDED
+
+class ImportMask
+{
+
+public:
+    int width, height;
+
+    ImportMask() { width=0; height=0; }
+
+
+    virtual real apply(int /*x*/, int /*y*/, real val) { return val; }
+};
+
+
+class RealMask : public ImportMask
+{
+    real *data;
+
+public:
+    RealMask(int width, int height);
+
+    inline real& ref(int x, int y);
+    real apply(int x, int y, real val) { return ref(x,y) * val; }
+};
+
+
+inline RealMask::RealMask(int w, int h)
+{
+    width = w;
+    height = h;
+    data = (real *)calloc(w*h, sizeof(real));
+}
+
+inline real& RealMask::ref(int i, int j)
+{
+#ifdef SAFETY
+	assert(i>=0); assert(j>=0); assert(i<width); assert(j<height);
+#endif
+
+	return data[j*width + i];
+}
+
+extern RealMask *readMask(istream&);
+
+#endif

src/Prep/Terra/Quadedge.cc

@@ -0,0 +1,82 @@
+#include <stdlib.h>
+#include <iostream.h>
+
+#include "Quadedge.h"
+
+Edge::Edge(const Edge&)
+{
+    cerr << "Edge: Edge assignments are forbidden." << endl;
+    exit(1);
+}
+
+
+Edge::Edge(Edge *prev)
+{
+    qprev = prev;
+    prev->qnext = this;
+
+    lface = NULL;
+    token = 0;
+}
+
+Edge::Edge()
+{
+    Edge *e0 = this;
+    Edge *e1 = new Edge(e0);
+    Edge *e2 = new Edge(e1);
+    Edge *e3 = new Edge(e2);
+
+    qprev = e3;
+    e3->qnext = e0;
+
+    e0->next = e0;
+    e1->next = e3;
+    e2->next = e2;
+    e3->next = e1;
+
+    lface = NULL;
+    token = 0;
+}
+
+Edge::~Edge()
+{
+    if( qnext )
+    {
+	Edge *e1 = qnext;
+	Edge *e2 = qnext->qnext;
+	Edge *e3 = qprev;
+
+#ifdef SAFETY
+	qnext = NULL;
+
+	token = -69;
+	e1->token = -69;
+	e2->token = -69;
+	e3->token = -69;
+#endif
+	e1->qnext = NULL;
+	e2->qnext = NULL;
+	e3->qnext = NULL;
+
+	delete e1;
+	delete e2;
+	delete e3;
+    }
+}
+
+
+void splice(Edge *a, Edge *b)
+{
+    Edge *alpha = a->Onext()->Rot();
+    Edge *beta  = b->Onext()->Rot();
+
+    Edge *t1 = b->Onext();
+    Edge *t2 = a->Onext();
+    Edge *t3 = beta->Onext();
+    Edge *t4 = alpha->Onext();
+
+    a->next = t1;
+    b->next = t2;
+    alpha->next = t3;
+    beta->next = t4;
+}

src/Prep/Terra/Quadedge.h

@@ -0,0 +1,82 @@
+#ifndef QUADEDGE_INCLUDED // -*- C++ -*-
+#define QUADEDGE_INCLUDED
+
+#include "Geom.h"
+
+class Triangle;
+
+class Edge : public Labelled {
+private:
+    Edge *qnext, *qprev;
+
+    Edge(Edge *prev);
+
+protected:
+    Vec2 *data;
+    Edge *next;
+
+    Triangle *lface;
+
+public:
+    Edge();
+    Edge(const Edge&);
+    ~Edge();
+
+    //
+    // Primitive methods
+    //
+    Edge *Onext() const { return next; }
+    Edge *Sym() const { return qnext->qnext; }
+    Edge *Rot() const { return qnext; }
+    Edge *invRot() const { return qprev; }
+
+    //
+    // Synthesized methods
+    //
+    Edge *Oprev() const { return Rot()->Onext()->Rot(); }
+    Edge *Dnext() const { return Sym()->Onext()->Sym(); }
+    Edge *Dprev() const { return invRot()->Onext()->invRot(); }
+
+    Edge *Lnext() const { return invRot()->Onext()->Rot(); }
+    Edge *Lprev() const { return Onext()->Sym(); }
+    Edge *Rnext() const { return Rot()->Onext()->invRot(); }
+    Edge *Rprev() const { return Sym()->Onext(); }
+
+
+    Vec2& Org() const { return *data; }
+    Vec2& Dest() const { return *Sym()->data; }
+
+    Triangle *Lface() const { return lface; }
+    void set_Lface(Triangle *t) { lface = t; }
+
+    void EndPoints(Vec2& org, Vec2& dest)
+    {
+	data = &org;
+	Sym()->data = &dest;
+    }
+
+    //
+    // The fundamental topological operator
+    friend void splice(Edge *a, Edge *b);
+};
+
+
+inline boolean rightOf(const Vec2& x, const Edge *e)
+{
+    return rightOf(x, e->Org(), e->Dest());
+}
+
+inline boolean leftOf(const Vec2& x, const Edge *e)
+{
+    return leftOf(x, e->Org(), e->Dest());
+}
+
+
+#ifdef IOSTREAMH
+inline ostream& operator<<(ostream& out, const Edge *e)
+{
+    return out << "{ " << e->Org() << " ---> " << e->Dest() << " }";
+}
+#endif
+
+#endif

src/Prep/Terra/README.html

@@ -0,0 +1,341 @@
+<HTML>
+<HEAD>
+<TITLE>Terra Height Field Simplification Software</TITLE>
+<!-- Created by: Michael Garland, 17-Jan-1996 -->
+<!-- Changed by: Michael Garland, 19-Jan-1996 -->
+</HEAD>
+
+<BODY>
+<H1>Terra Height Field Simplification Software</H1>
+
+<P>This software is in the public domain and is provided <STRONG>AS
+IS</STRONG>. Use it at <STRONG>YOUR OWN RISK</STRONG>.
+
+<P><HR>
+
+<P>This is Terra, the successor to Scape.  Not only is this software
+<STRONG>UNSUPPORTED</STRONG>, but it is in the early stages of
+development.  This software is clearly incomplete and may still have
+bugs.  However, the basic mechanisms do function properly; I have run
+several large terrains through Terra, and the results were correct.
+
+<P>For up-to-date information on Terra, Scape, and related topics,
+watch:
+<PRE>
+	<A HREF="http://www.cs.cmu.edu/~garland/scape/">http://www.cs.cmu.edu/~garland/scape/</A>
+</PRE>
+
+<P><HR>
+
+<H2>Introduction</H2>
+
+<P>This is <STRONG>Terra</STRONG>, a program for generating polygonal
+approximations of terrains and other height fields.  It is based on
+algorithms described in:
+<BLOCKQUOTE>
+<CITE>
+<A HREF="http://www.cs.cmu.edu/~garland/scape/scape.ps.gz">Fast
+Polygonal Approximation of Terrains
+and Height Fields</A>,<BR>by Michael Garland and Paul Heckbert
+(Technical Report CMU-CS-95-181).
+</CITE><BR>
+A <A HREF="http://www.cs.cmu.edu/~garland/scape/color.ps.gz">color plate</A>
+is included separately.
+</BLOCKQUOTE>
+
+<P>The <B>Scape</B> package is the companion software for this paper.
+It was written with speed and memory efficiency as the primary
+concerns.  Although it was designed strictly for the testing of the
+algorithms described in our paper, we made it available so that people
+interested in our results could examine it.  We also hoped that it
+might be of interest to people who were attempting to build terrain
+approximations.  However, Scape is not particularly easy to use and
+the code is definitely less than aesthetically pleasing.
+
+<P>I wrote Terra because I was dissatisfied with Scape.  I wanted code
+which was better structured and programs which were easier to use.
+Terra will also provide more features than Scape.
+
+<P><B>Disclaimer:</B> Please remember that both Terra and Scape are
+unsupported programs that I tinker with in my spare time.  As such,
+development is typically sporadic.
+
+
+<H3>Brief feature summary</H3>
+<DL>
+<DT><EM>PGM input files</EM>
+<DD>Terra uses the PGM file format for data input.  At first, this
+might seem odd; however, it has several advantages.  PGM is a standard
+format which provides for both textual (i.e., human editable) and
+binary data files.  Since PGM is an image file format, height field
+data is directly viewable with most available image viewers.  Plus,
+there are many programs available to perform various manipulations on
+PGM files.
+
+<P>In particular, many of the tools provided by the <B>NetPBM</B> package
+can be used to manipulate PGM terrain data.  The NetPBM package can be
+found online at:
+<PRE>
+	<A HREF="http://wuarchive.wustl.edu/graphics/graphics/packages/NetPBM/">http://wuarchive.wustl.edu/graphics/graphics/packages/NetPBM/</A>
+</PRE>
+or by anonymous ftp to:
+<PRE>
+	<A HREF="ftp://wuarchive.wustl.edu/graphics/graphics/packages/NetPBM/">wuarchive.wustl.edu:/graphics/graphics/packages/NetPBM</A>
+</PRE>
+
+<P><DT><EM>Flexible output</EM>
+<DD>Terra can output its approximations in several different formats.
+The supported formats are described below.
+
+<P><DT><EM>Extended approximation features</EM>
+<DD>Terra supports preinsertion scripts and importance masks.  See
+below for details on what these are and how they work.
+
+<P><DT><EM>Portable graphics</EM>
+<DD>The interactive program, <TT>xterra</TT>, uses the GLUT library
+for windowing and OpenGL for rendering.  In theory, this should make
+it portable to machines other than SGI's.
+</DL>
+
+<H3>Currently absent features</H3>
+<P>All these features are currently missing.  Ideally, I would like
+to include these and other features.  In reality, what gets done and
+how fast it gets done might be highly variable.
+<DL>
+<DT><EM>Multiple simultaneous height fields</EM>
+<DD>Scape provided for approximating height fields with RGB textures
+applied to them.  Ultimately, Terra will support the approximation of
+arbitrarily many simultaneous height fields.  Currently, however,
+Terra will only accept the input of a <EM>single</EM> height field.
+
+<P><DT><EM>Data-dependent triangulation</EM>
+<DD>One of the significant features of Scape was the option to use
+data-dependent triangulation.  This triangulation scheme has not yet
+been ported to Terra.
+
+<P><DT><EM>Data import facilities</EM>
+<DD>I would like to be able to import USGS DEM data into Terra.
+Although correcting for the spherical mapping of USGS data is beyond
+the scope of this project, Terra is in need of some simple tools to
+facilitate conversion of USGS data.  I definitely don't have the time
+to write these tools, so if someone wants to suggest some reasonable
+ones, please let me know.
+
+<P><DT><EM>Robust PGM input routines</EM>
+<DD>The PGM input routines are rather fragile at present.  You should
+make sure that there is no extraneous white space and <EM>no
+comments</EM> in the file.  For example, the <TT>xv</TT> program will
+insert a comment in the PGM file citing itself as the creator of the
+file.  You will need to remove this comment from the file.
+</DL>
+
+<H2>Installation</H2>
+
+<OL>
+<P><LI>In order to compile the interactive version of Terra (<TT>xterra</TT>),
+       you must obtain the GLUT library.  It can be found at:
+
+<PRE>
+	<A HREF="http://www.sgi.com/Technology/openGL/glut.html">http://www.sgi.com/Technology/openGL/glut.html</A>
+</PRE>
+	       
+       <P>or by anonymous ftp at:
+
+<PRE>
+	<A HREF="ftp://sgigate.sgi.com/pub/opengl/xjournal/GLUT/">sgigate.sgi.com:/pub/opengl/xjournal/GLUT</A>
+</PRE>
+
+      <P>NOTE: For proper viewing, <TT>xterra</TT> needs to set the aspect ratio
+       of its windows.  This is currently not supported via GLUT.
+       Therefore, I've had to hack things a bit by accessing GLUT
+       internals.  The file gui.cc includes the glutint.h header.
+       This is not installed by GLUT.  You must install this manually.
+       Again, this is only of importance if you want to build <TT>xterra</TT>.
+
+
+<P><LI>Edit the Makefile for local customization.  Basically, you
+       should set the compilation flags, the libraries you need, and
+       the location of GLUT if you are compiling <TT>xterra</TT>.
+
+<P><LI>Typing '<TT>make</TT>' will build both <TT>terra</TT> and
+       <TT>xterra</TT>.  However, you
+       can selectively compile either of them (e.g., '<TT>make terra</TT>').
+</OL>
+
+
+
+<H2>Using Terra</H2>
+
+<P>The Terra software provides two programs: <TT>terra</TT>, a simple
+batch program, and <TT>xterra</TT>, an interactive application which
+uses OpenGL to display the surface approximation being
+constructed. Both programs utilize the same command line interface;
+<TT>xterra</TT> simply layers an interactive interface on top of
+<TT>terra</TT>.  This section contains an outline of the operation of
+<TT>terra</TT>.  All this information applies to <TT>xterra</TT> as well.
+
+<P>The operation of Terra goes through a simple series of phases:
+<OL>
+<LI><EM>Data input</EM>.<BR>
+This will load the terrain data and the importance mask (if specified).
+
+<LI><EM>Pre-insertion</EM>.<BR>
+If you have specified a pre-insertion script, it will be executed at
+this point.
+
+<LI><EM>Greedy insertion</EM>.<BR>
+The iterative greedy insertion procedure will begin.  Terra will
+continue selecting points until it meets the goals that you have
+specified.
+
+<LI><EM>Output</EM>.<BR>
+Terra will output the approximation it has constructed.  You have a
+choice of a handful of different formats.
+</OL>
+
+<P>Currently, all the information Terra needs to run is communicated
+on the command line.  The correct usage of Terra is:
+<PRE>
+	usage: terra &lt;options&gt; filename
+	       where &lt;options&gt; is some combination of: 
+	-e &lt;thresh&gt;      Sets the tolerable error threshold
+	-p &lt;count&gt;       Sets the maximum number of allowable points
+	-o &lt;file&gt; &lt;type&gt; When finished, output the approximation to &lt;file&gt;.
+	                 Valid types:  tin, eps, dem, obj
+	-m &lt;file&gt;        Load the importance mask from &lt;file&gt;
+	-s &lt;file&gt;        Execute preinsertion script from &lt;file&gt;
+</PRE>
+
+<P>The error threshold and point limit set the termination criteria.
+Terra will continue adding points as long as it it below the point
+limit and above the error threshold.  The default error threshold is
+0; the default point limit is the total number of points in the input
+grid.
+
+<P>The type of output desired is also specified on the command line.
+The <TT>eps</TT> output format simply generates an Encapsulated
+PostScript rendering of the approximation mesh.  The <TT>dem</TT>
+output format creates an approximate DEM from the approximate mesh.
+This can be useful for comparison with the original.  Both the
+<TT>tin</TT> and <TT>obj</TT> output formats generate 3D surfaces.
+The <TT>obj</TT> format is just the Wavefront <TT>.OBJ</TT> file
+format.  The <TT>tin</TT> format is a very simple model description;
+it is a series of lines where each line is of the form:
+<PRE>
+	t x1 y1 z1 x2 y2 z2 x3 y3 z3
+</PRE>
+Each such line describes a triangle with the three corners (x1,y1,z1)
+(x2,y2,z2) (x3,y3,z3).
+
+<P>The remaining options, involving importance masks and preinsertion
+scripts, are described in detail below.
+
+<H2>Xterra User Interface</H2>
+
+<P><TT>xterra</TT> accepts the same options as <TT>terra</TT> and
+operates in much the same way.  It introduces one extra command line option:
+<PRE>
+	-h &lt;factor&gt;      Sets the height scaling factor.  For example,
+	                 if grid points are 25m apart, use a factor of 0.04
+</PRE>
+
+<P>This is used to properly scale the display of the height field in
+3D.  The input to Terra is specified in pixel coordinates; data
+samples are taken at integral pixel coordinates.  However, the height
+values are probably not given in pixel coordinates.  So, for display
+purposes, the height values are scaled by a constant factor to account
+for this loss of units in Terra.
+
+<P>When <TT>xterra</TT> starts up, it performs any preinsertion
+actions that you request, and then it displays two windows: a mesh
+view and a surface view.  It does not perform greedy insertion until
+told to do so.  The mesh view provides a 2D view of the triangulation
+being generated for the approximation of the height field.  The
+surface view displays the approximate surface in 3D.  The interaction
+with these windows is currently quite simple and will probably change
+in the future, but here is an outline of how they work at the moment:
+<PRE>
+	Surface view:
+	   Left mouse drag   :   spin the surface
+	   Middle mouse drag :   translate the camera side to side
+	   Right mouse drag  :   move the camera in and out
+
+	Mesh view:
+	   Left mouse click   :  insert a point at the mouse location
+	   Middle mouse click :  run greedy insertion until goal is met
+	   Right mouse click  :  popup menu -- allows output
+</PRE>
+
+<H2>Providing Input for Terra</H2>
+
+<P>As stated above, Terra uses PGM files to read and write height
+field data.  Unfortunately, Terra does not as yet provide any direct
+means of acquiring PGM data.  For now, you will have to use the
+conversion software provided with Scape.  The
+<TT>STM-tools/stm2pgm</TT> utility distributed with Scape can convert
+Scape's STM file format into PGM's appropriate for use with Terra.
+Given an STM file,
+<PRE>
+	<TT>stm2pgm sample.stm exact &gt; sample.pgm</TT>
+</PRE>
+will generate a PGM file.  Note that the keyword <TT>exact</TT> is
+very important.  Don't forget it!  The resulting file will be textual,
+so you can even edit by hand if you need to.
+
+<H2>Importance Masks</H2>
+
+<P>One of the new features in Terra not found in Scape is the use of
+importance masks.  In order to determine the next point for insertion,
+Terra ranks the available points by an importance measure.  Importance
+masks allow you to bias this ranking.  For each data point, the mask
+assigns a weight in the range [0..1] which multiplies the computed
+importance value.
+
+<P>Importance masks are PGM files, just like the height field input.
+However, their interpretation is slightly different.  The input values
+are all integers.  They are scaled such that their maximum value will
+be mapped to 1.  One significant point is that this maximum value is
+taken from the PGM header, not determined from the data.  Therefore,
+by controlling the stated "maximum", you have much greater flexibility
+over the mapping of PGM pixel values to importance mask weights.
+
+<P>Currently, no means for constructing importance masks is provided.
+
+<H2>Preinsertion Scripts</H2>
+
+<P>The other new feature of Terra is its support for preinsertion
+scripts.  An important feature of the greedy insertion algorithm is
+that essentially provides for progressive refinement of the
+approximation.  Thus, the initial approximation can be arbitrary.  The
+preinsertion scripts allow you set up an approximation before greedy
+insertion begins.
+
+<P>A preinsertion script is a series of lines, each of the form:
+<PRE>
+	&lt;op&gt; X Y
+</PRE>
+The values for <TT>X</TT> and <TT>Y</TT> are the coordinates of a
+particular data point.  The currently supported opcodes are:
+<PRE>
+	<B>s</B> -- Select this point for use in the approximation and insert
+	     it into the current mesh.
+
+	<B>i</B> -- Mark this point as one to be ignored.  Terra will never
+	     process this point for insertion or evaluation.
+
+	<B>u</B> -- Mark the height value at this point as unknown.
+</PRE>
+Currently, Terra makes no distinction between points to be ignored and
+points whose height value is unknown; it ignores them equally.
+
+<P><HR>
+
+<P>January 19, 1996
+<ADDRESS>
+<A HREF="http://www.cs.cmu.edu/~garland/home.html">Michael Garland</A>
+<BR>
+<A HREF="mailto:garland@cs.cmu.edu">garland@cs.cmu.edu</A>
+</ADDRESS>
+
+</BODY>
+</HTML>

src/Prep/Terra/Subdivision.cc

@@ -0,0 +1,429 @@
+#include <stdlib.h>
+#include <iostream.h>
+#include <assert.h>
+
+#include "Subdivision.h"
+
+
+
+Edge *Subdivision::makeEdge(Vec2& org, Vec2& dest)
+{
+    Edge *e = new Edge();
+    e->EndPoints(org, dest);
+    return e;
+}
+
+Edge *Subdivision::makeEdge()
+{
+    return new Edge();
+}
+
+void Subdivision::initMesh(const Vec2& A,const Vec2& B,
+			   const Vec2& C,const Vec2& D)
+{
+    Vec2& a = A.clone();
+    Vec2& b = B.clone();
+    Vec2& c = C.clone();
+    Vec2& d = D.clone();
+
+    Edge *ea = makeEdge();
+    ea->EndPoints(a, b);
+
+    Edge *eb = makeEdge();
+    splice(ea->Sym(), eb);
+    eb->EndPoints(b, c);
+
+    Edge *ec = makeEdge();
+    splice(eb->Sym(), ec);
+    ec->EndPoints(c, d);
+
+    Edge *ed = makeEdge();
+    splice(ec->Sym(), ed);
+    ed->EndPoints(d, a);
+    splice(ed->Sym(), ea);
+
+    Edge *diag = makeEdge();
+    splice(ed->Sym(),diag);
+    splice(eb->Sym(),diag->Sym());
+    diag->EndPoints(a,c);
+
+    startingEdge = ea;
+
+    first_face = NULL;
+
+    makeFace(ea->Sym()).update(*this);
+    makeFace(ec->Sym()).update(*this);
+
+}
+
+
+
+void Subdivision::deleteEdge(Edge *e)
+{
+    splice(e, e->Oprev());
+    splice(e->Sym(), e->Sym()->Oprev());
+
+    delete e;
+}
+
+Edge *Subdivision::connect(Edge *a, Edge *b)
+{
+    Edge *e = makeEdge();
+    splice(e, a->Lnext());
+    splice(e->Sym(), b);
+    e->EndPoints(a->Dest(), b->Org());
+
+    return e;
+}
+
+void Subdivision::swap(Edge *e)
+{
+    Triangle *f1 = e->Lface();
+    Triangle *f2 = e->Sym()->Lface();
+
+    Edge *a = e->Oprev();
+    Edge *b = e->Sym()->Oprev();
+
+    splice(e, a);
+    splice(e->Sym(), b);
+    splice(e, a->Lnext());
+    splice(e->Sym(), b->Lnext());
+    e->EndPoints(a->Dest(), b->Dest());
+
+
+    f1->reshape(e);
+    f2->reshape(e->Sym());
+}
+
+
+
+//
+// Subdivision iterators
+//
+
+static unsigned int timestamp = 0;
+
+static void overEdge(Edge *e, edge_callback fn, void *closure)
+{
+    if( e->token != timestamp )
+    {
+	e->token = timestamp;
+	e->Sym()->token = timestamp;
+
+	(*fn)(e, closure);
+
+	overEdge(e->Onext(), fn, closure);
+	overEdge(e->Oprev(), fn, closure);
+	overEdge(e->Dnext(), fn, closure);
+	overEdge(e->Dprev(), fn, closure);
+    }
+}
+
+void Subdivision::overEdges(edge_callback fn, void *closure)
+{
+    if( ++timestamp == 0 )
+	timestamp = 1;
+
+    overEdge(startingEdge, fn, closure);
+}
+
+void Subdivision::overFaces(face_callback fn, void *closure)
+{
+    Triangle *t = first_face;
+
+    while( t )
+    {
+	(*fn)(*t, closure);
+	t = t->getLink();
+    }
+}
+
+
+
+
+//
+// Random predicates
+//
+
+boolean Subdivision::ccwBoundary(const Edge *e)
+{
+    return !rightOf(e->Oprev()->Dest(), e);
+}
+
+
+boolean Subdivision::onEdge(const Vec2& x, Edge *e)
+{
+    real t1, t2, t3;
+
+    t1 = (x - e->Org()).length();
+    t2 = (x - e->Dest()).length();
+
+    if (t1 < EPS || t2 < EPS)
+	return True;
+
+    t3 = (e->Org() - e->Dest()).length();
+
+    if (t1 > t3 || t2 > t3)
+	return False;
+
+    Line line(e->Org(), e->Dest());
+    return (fabs(line.eval(x)) < EPS);
+}
+
+
+boolean Subdivision::isInterior(Edge *e)
+//
+// Tests whether e is an interior edge.
+// 
+// WARNING: This topological test will not work if the boundary is
+//          a triangle.  This is not a problem here; the boundary is
+//          always a rectangle.  But if you try to adapt this code, please
+//          keep this in mind.
+{
+    return (e->Lnext()->Lnext()->Lnext() == e &&
+            e->Rnext()->Rnext()->Rnext() == e );
+}
+
+boolean Subdivision::shouldSwap(const Vec2& x, Edge *e)
+{
+        Edge *t = e->Oprev();
+        return inCircle(e->Org(), t->Dest(), e->Dest(), x);
+}
+
+
+
+
+
+Edge *Subdivision::locate(const Vec2& x, Edge *start)
+{
+    Edge *e = start;
+    real t = triArea(x, e->Dest(), e->Org());
+
+    if (t>0) {                  // x is to the right of edge e
+        t = -t;
+        e = e->Sym();
+    }
+
+
+    while (True)
+    {
+        Edge *eo = e->Onext();
+        Edge *ed = e->Dprev();
+
+        real to = triArea(x, eo->Dest(), eo->Org());
+        real td = triArea(x, ed->Dest(), ed->Org());
+
+        if (td>0)                       // x is below ed
+            if (to>0 || to==0 && t==0) {// x is interior, or origin endpoint
+                startingEdge = e;
+                return e;
+            }
+            else {                      // x is below ed, below eo
+                t = to;
+                e = eo;
+            }
+        else                            // x is on or above ed
+            if (to>0)                   // x is above eo
+                if (td==0 && t==0) {    // x is destination endpoint
+                    startingEdge = e;
+                    return e;
+                }
+                else {                  // x is on or above ed and above eo
+                    t = td;
+                    e = ed;
+                }
+            else                        // x is on or below eo
+                if (t==0 && !leftOf(eo->Dest(), e))
+		    // x on e but subdiv. is to right
+                    e = e->Sym();
+                else if (random()&1) {  // x is on or above ed and
+                    t = to;             // on or below eo; step randomly
+                    e = eo;
+                }
+                else {
+                    t = td;
+                    e = ed;
+                }
+    }
+}
+
+
+Edge *Subdivision::spoke(Vec2& x, Edge *e)
+{
+    Triangle *new_faces[4];
+    int facedex = 0;
+
+    //
+    // NOTE: e is the edge returned by locate(x)
+    //
+
+    if ( (x == e->Org()) || (x == e->Dest()) ) {
+        // point is already in the mesh
+	//
+        cerr << "WARNING: Tried to reinsert point: " << x << endl;
+	cerr << "         org: " << e->Org() << endl;
+	cerr << "        dest: " << e->Dest() << endl;
+        return NULL;
+    }
+
+    Edge *boundary_edge = NULL;
+
+    Triangle *lface = e->Lface();
+    lface->dontAnchor(e);
+    new_faces[facedex++] = lface;
+
+    if( onEdge(x,e) )
+    {
+	if( ccwBoundary(e) ) {
+	    //
+	    // e lies on the boundary
+	    // Defer deletion until after new edges are added.
+	    boundary_edge = e;
+	}
+	else {
+	    Triangle *sym_lface = e->Sym()->Lface();
+	    new_faces[facedex++] = sym_lface;
+	    sym_lface->dontAnchor(e->Sym());
+
+
+	    e = e->Oprev();
+	    deleteEdge(e->Onext());
+	}
+    }
+    else
+    {
+	// x lies within the Lface of e
+    }
+
+    Edge *base = makeEdge(e->Org(), x.clone());
+
+    splice(base, e);
+
+    startingEdge = base;
+    do {
+	base = connect(e, base->Sym());
+	e = base->Oprev();
+    } while( e->Lnext() != startingEdge );
+
+    if( boundary_edge )
+	deleteEdge(boundary_edge);
+
+    // Update all the faces in our new spoked polygon.
+    // If point x on perimeter, then don't add an exterior face
+
+
+    base = boundary_edge ? startingEdge->Rprev() : startingEdge->Sym();
+    do {
+	if( facedex )
+	    new_faces[--facedex]->reshape(base);
+	else
+	    makeFace(base);
+
+	base = base->Onext();
+
+    } while( base != startingEdge->Sym() );
+
+    return startingEdge;
+}
+
+
+//
+// s is a spoke pointing OUT from x
+//
+void Subdivision::optimize(Vec2& x, Edge *s)
+{
+    Edge *start_spoke = s;
+    Edge *spoke = s;
+
+    do {
+
+	Edge *e = spoke->Lnext();
+	Edge *t = e->Oprev();
+
+	if( isInterior(e) && shouldSwap(x, e) )
+	    swap(e);
+	else
+	{
+	    spoke = spoke->Onext();
+	    if( spoke == start_spoke )
+		break;
+	}
+
+    } while( True );
+
+    //
+    // Now, update all the triangles
+
+    spoke = start_spoke;
+
+    do {
+	Edge *e = spoke->Lnext();
+	Triangle *t = e->Lface();
+
+	if( t ) t->update(*this);
+
+	spoke = spoke->Onext();
+    } while( spoke != start_spoke );
+}
+
+Edge *Subdivision::insert(Vec2& x, Triangle *tri)
+{
+    Edge *e = tri?locate(x, tri->getAnchor()):locate(x);
+
+    Edge *start_spoke = spoke(x, e);
+
+    if( start_spoke )
+	optimize(x, start_spoke->Sym());
+
+    return start_spoke;
+}
+
+
+
+Triangle *Subdivision::allocFace(Edge *e)
+{
+    return new Triangle(e);
+}
+
+Triangle& Subdivision::makeFace(Edge *e)
+{
+    Triangle *t = allocFace(e);
+
+    first_face = t->linkTo(first_face);
+
+    return *t;
+}
+
+
+
+
+
+
+
+void Triangle::dontAnchor(Edge *e)
+{
+    if( anchor == e )
+    {
+	anchor = e->Lnext();
+    }
+}
+
+void Triangle::reshape(Edge *e)
+{
+    anchor = e;
+
+    e->set_Lface(this);
+    e->Lnext()->set_Lface(this);
+    e->Lprev()->set_Lface(this);
+
+}
+
+
+void Triangle::update(Subdivision&)
+// called by reshape to update stuff
+//
+// the default method will do nothing
+{
+
+
+}

src/Prep/Terra/Subdivision.h

@@ -0,0 +1,95 @@
+#ifndef SUBDIVISION_INCLUDED // -*- C++ -*-
+#define SUBDIVISION_INCLUDED
+
+#include "Quadedge.h"
+
+class Subdivision;
+
+class Triangle : public Labelled {
+    Edge *anchor;
+    Triangle *next_face;
+
+public:
+
+    Triangle(Edge *e, int t=0)
+    {
+	token = t;
+	reshape(e);
+    }
+
+    Triangle *linkTo(Triangle *t) { next_face = t; return this; }
+    Triangle *getLink() { return next_face; }
+    Edge *getAnchor() { return anchor; }
+    void dontAnchor(Edge *e);
+
+    void reshape(Edge *e);
+
+    virtual void update(Subdivision&); // called  to update stuff
+
+    const Vec2& point1() const { return anchor->Org(); }
+    const Vec2& point2() const { return anchor->Dest(); }
+    const Vec2& point3() const { return anchor->Lprev()->Org(); }
+};
+
+typedef void (*edge_callback)(Edge *, void *);
+typedef void (*face_callback)(Triangle&, void *);
+
+
+class Subdivision {
+private:
+    Edge *startingEdge;
+    Triangle *first_face;
+
+protected:
+    void initMesh(const Vec2&, const Vec2&, const Vec2&, const Vec2&);
+    Subdivision() { }
+
+    Edge *makeEdge();
+    Edge *makeEdge(Vec2& org, Vec2& dest);
+
+    virtual Triangle *allocFace(Edge *e);
+    Triangle& makeFace(Edge *e);
+
+
+    void deleteEdge(Edge *);
+    Edge *connect(Edge *a, Edge *b);
+    void swap(Edge *e);
+
+    //
+    // Some random functions
+    boolean ccwBoundary(const Edge *e);
+    boolean onEdge(const Vec2&, Edge *);
+
+public:
+
+    Subdivision(Vec2& a, Vec2& b, Vec2& c, Vec2& d) { initMesh(a,b,c,d); }
+
+    //
+    // virtual functions for customization
+    virtual boolean shouldSwap(const Vec2&, Edge *);
+
+
+    boolean isInterior(Edge *);
+
+    Edge *spoke(Vec2&, Edge *e);
+    void optimize(Vec2&, Edge *);
+
+    Edge *locate(const Vec2& x)  { return locate(x, startingEdge); }
+    Edge *locate(const Vec2&, Edge *hint);
+    Edge *insert(Vec2&, Triangle *t=NULL);
+
+    void overEdges(edge_callback, void *closure=NULL);
+    void overFaces(face_callback, void *closure=NULL);
+};
+
+
+#ifdef IOSTREAMH
+inline ostream& operator<<(ostream& out, Triangle& t)
+{
+    return out << "Triangle("<< t.point1() << " " << t.point2() << " "
+	       << t.point3() << ")";
+}
+#endif
+
+
+#endif

src/Prep/Terra/Vec2.h

@@ -0,0 +1,171 @@
+#ifndef VEC2_INCLUDED // -*- C++ -*-
+#define VEC2_INCLUDED
+
+class Vec2 {
+protected:
+    real elt[2];
+
+    inline void copy(const Vec2& v);
+
+public:
+    // Standard constructors
+    Vec2(real x=0, real y=0) { elt[0]=x; elt[1]=y; }
+    Vec2(const Vec2& v) { copy(v); }
+    Vec2(const real *v) { elt[0]=v[0]; elt[1]=v[1]; }
+    Vec2& clone() const { return *(new Vec2(elt[0], elt[1])); }
+
+    // Access methods
+    real& operator()(int i)             { return elt[i]; }
+    const real& operator()(int i) const { return elt[i]; }
+    real& operator[](int i)             { return elt[i]; }
+    const real& operator[](int i) const { return elt[i]; }
+
+    // Assignment methods
+    inline Vec2& operator=(const Vec2& v);
+    inline Vec2& operator+=(const Vec2& v);
+    inline Vec2& operator-=(const Vec2& v);
+    inline Vec2& operator*=(real s);
+    inline Vec2& operator/=(real s);
+
+    // Arithmetic methods
+    inline Vec2 operator+(const Vec2& v) const;
+    inline Vec2 operator-(const Vec2& v) const;
+    inline Vec2 operator-() const;
+
+    inline Vec2 operator*(real s) const;
+    inline Vec2 operator/(real s) const;
+    inline real operator*(const Vec2& v) const;
+
+#ifdef IOSTREAMH
+    // Input/Output methods
+    friend ostream& operator<<(ostream&, const Vec2&);
+    friend istream& operator>>(istream&, Vec2&);
+#endif
+
+    // Additional vector methods
+    inline real length();
+    inline real norm();
+    inline real norm2();
+
+    inline real unitize();
+
+    inline int operator==(const Vec2& v) const
+    {
+	return (*this - v).norm2() < EPS2;
+    }
+};
+
+
+
+inline void Vec2::copy(const Vec2& v)
+{
+    elt[0]=v.elt[0]; elt[1]=v.elt[1];
+}
+
+inline Vec2& Vec2::operator=(const Vec2& v)
+{
+    copy(v);
+    return *this;
+}
+
+inline Vec2& Vec2::operator+=(const Vec2& v)
+{
+    elt[0] += v[0];
+    elt[1] += v[1];
+    return *this;
+}
+
+inline Vec2& Vec2::operator-=(const Vec2& v)
+{
+    elt[0] -= v[0];
+    elt[1] -= v[1];
+    return *this;
+}
+
+inline Vec2& Vec2::operator*=(real s)
+{
+    elt[0] *= s;
+    elt[1] *= s;
+    return *this;
+}
+
+inline Vec2& Vec2::operator/=(real s)
+{
+    elt[0] /= s;
+    elt[1] /= s;
+    return *this;
+}
+
+///////////////////////
+
+inline Vec2 Vec2::operator+(const Vec2& v) const
+{
+    Vec2 w(elt[0]+v[0], elt[1]+v[1]);
+    return w;
+}
+
+inline Vec2 Vec2::operator-(const Vec2& v) const
+{
+    Vec2 w(elt[0]-v[0], elt[1]-v[1]);
+    return w;
+}
+
+inline Vec2 Vec2::operator-() const
+{
+    return Vec2(-elt[0], -elt[1]);
+}
+
+inline Vec2 Vec2::operator*(real s) const
+{
+    Vec2 w(elt[0]*s, elt[1]*s);
+    return w;
+}
+
+inline Vec2 Vec2::operator/(real s) const
+{
+    Vec2 w(elt[0]/s, elt[1]/s);
+    return w;
+}
+
+inline real Vec2::operator*(const Vec2& v) const
+{
+    return elt[0]*v[0] + elt[1]*v[1];
+}
+
+inline real Vec2::length()
+{
+    return norm();
+}
+
+inline real Vec2::norm()
+{
+    return sqrt(elt[0]*elt[0] + elt[1]*elt[1]);
+}
+
+inline real Vec2::norm2()
+{
+    return elt[0]*elt[0] + elt[1]*elt[1];
+}
+
+inline real Vec2::unitize()
+{
+    real l=norm();
+    if( l!=1.0 )
+	(*this)/=l;
+    return l;
+}
+
+#ifdef IOSTREAMH
+inline ostream& operator<<(ostream& out, const Vec2& v)
+{
+    return out << "[" << v[0] << " " << v[1] << "]";
+}
+
+inline istream& operator>>(istream& in, Vec2& v)
+{
+    return in >> "[" >> v[0] >> v[1] >> "]";
+}
+#endif
+
+
+#endif

src/Prep/Terra/Vec3.h

@@ -0,0 +1,181 @@
+#ifndef VEC3_INCLUDED // -*- C++ -*-
+#define VEC3_INCLUDED
+
+class Vec3 {
+protected:
+    real elt[3];
+
+    inline void copy(const Vec3& v);
+
+public:
+    // Standard constructors
+    Vec3(real x=0, real y=0, real z=0) { elt[0]=x; elt[1]=y; elt[2]=z; }
+    Vec3(const Vec2& v, real z) { elt[0]=v[0]; elt[1]=v[1]; elt[2]=z; }
+    Vec3(const Vec3& v) { copy(v); }
+    Vec3(const real *v) { elt[0]=v[0]; elt[1]=v[1]; elt[2]=v[2]; }
+
+    // Access methods
+    real& operator()(int i)             { return elt[i]; }
+    const real& operator()(int i) const { return elt[i]; }
+    real& operator[](int i)             { return elt[i]; }
+    const real& operator[](int i) const { return elt[i]; }
+
+    // Assignment methods
+    inline Vec3& operator=(const Vec3& v);
+    inline Vec3& operator+=(const Vec3& v);
+    inline Vec3& operator-=(const Vec3& v);
+    inline Vec3& operator*=(real s);
+    inline Vec3& operator/=(real s);
+
+    // Arithmetic methods
+    inline Vec3 operator+(const Vec3& v) const;
+    inline Vec3 operator-(const Vec3& v) const;
+    inline Vec3 operator-() const;
+
+    inline Vec3 operator*(real s) const;
+    inline Vec3 operator/(real s) const;
+    inline real operator*(const Vec3& v) const;
+    inline Vec3 operator^(const Vec3& v) const;
+
+
+#ifdef IOSTREAMH
+    // Input/Output methods
+    friend ostream& operator<<(ostream&, const Vec3&);
+    friend istream& operator>>(istream&, Vec3&);
+#endif
+
+    // Additional vector methods
+    inline real length();
+    inline real norm();
+    inline real norm2();
+
+    inline real unitize();
+};
+
+
+
+
+inline void Vec3::copy(const Vec3& v)
+{
+    elt[0]=v.elt[0]; elt[1]=v.elt[1]; elt[2]=v.elt[2];
+}
+
+inline Vec3& Vec3::operator=(const Vec3& v)
+{
+    copy(v);
+    return *this;
+}
+
+inline Vec3& Vec3::operator+=(const Vec3& v)
+{
+    elt[0] += v[0];
+    elt[1] += v[1];
+    elt[2] += v[2];
+    return *this;
+}
+
+inline Vec3& Vec3::operator-=(const Vec3& v)
+{
+    elt[0] -= v[0];
+    elt[1] -= v[1];
+    elt[2] -= v[2];
+    return *this;
+}
+
+inline Vec3& Vec3::operator*=(real s)
+{
+    elt[0] *= s;
+    elt[1] *= s;
+    elt[2] *= s;
+    return *this;
+}
+
+inline Vec3& Vec3::operator/=(real s)
+{
+    elt[0] /= s;
+    elt[1] /= s;
+    elt[2] /= s;
+    return *this;
+}
+
+///////////////////////
+
+inline Vec3 Vec3::operator+(const Vec3& v) const
+{
+    Vec3 w(elt[0]+v[0], elt[1]+v[1], elt[2]+v[2]);
+    return w;
+}
+
+inline Vec3 Vec3::operator-(const Vec3& v) const
+{
+    Vec3 w(elt[0]-v[0], elt[1]-v[1], elt[2]-v[2]);
+    return w;
+}
+
+inline Vec3 Vec3::operator-() const
+{
+    return Vec3(-elt[0], -elt[1], -elt[2]);
+}
+
+inline Vec3 Vec3::operator*(real s) const
+{
+    Vec3 w(elt[0]*s, elt[1]*s, elt[2]*s);
+    return w;
+}
+
+inline Vec3 Vec3::operator/(real s) const
+{
+    Vec3 w(elt[0]/s, elt[1]/s, elt[2]/s);
+    return w;
+}
+
+inline real Vec3::operator*(const Vec3& v) const
+{
+    return elt[0]*v[0] + elt[1]*v[1] + elt[2]*v[2];
+}
+
+inline Vec3 Vec3::operator^(const Vec3& v) const
+{
+    Vec3 w( elt[1]*v[2] - v[1]*elt[2],
+	   -elt[0]*v[2] + v[0]*elt[2],
+	    elt[0]*v[1] - v[0]*elt[1] );
+    return w;
+}
+
+inline real Vec3::length()
+{
+    return norm();
+}
+
+inline real Vec3::norm()
+{
+    return sqrt(elt[0]*elt[0] + elt[1]*elt[1] + elt[2]*elt[2]);
+}
+
+inline real Vec3::norm2()
+{
+    return elt[0]*elt[0] + elt[1]*elt[1] + elt[2]*elt[2];
+}
+
+inline real Vec3::unitize()
+{
+    real l=norm();
+    if( l!=1.0 )
+	(*this)/=l;
+    return l;
+}
+
+#ifdef IOSTREAMH
+inline ostream& operator<<(ostream& out, const Vec3& v)
+{
+    return out << "[" << v[0] << " " << v[1] << " " << v[2] << "]";
+}
+
+inline istream& operator>>(istream& in, Vec3& v)
+{
+    return in >> "[" >> v[0] >> v[1] >> v[2] >> "]";
+}
+#endif
+
+
+#endif

src/Prep/Terra/cmdline.cc

@@ -0,0 +1,173 @@
+#include <stdlib.h>
+#include <fstream.h>
+#include <string.h>
+#include "terra.h"
+
+GreedySubdivision *mesh;
+Map *DEM;
+
+static ImportMask default_mask;
+ImportMask *MASK;
+
+
+real error_threshold = 0.0;
+int  point_limit     = -1;
+real height_scale    = 1.0;
+FileFormat output_format;
+char *output_filename = NULL;
+char *mask_filename   = NULL;
+char *script_filename = NULL;
+
+static char *options = "e:p:h:o:m:s:";
+
+static char *usage_string =
+"-e <thresh>      Sets the tolerable error threshold\n"
+"-p <count>       Sets the maximum number of allowable points\n"
+"-h <factor>      Sets the height scaling factor.  For example,\n"
+"                 if grid points are 25m apart, use a factor of 0.04\n"
+"-o <file> <type> When finished, output the approximation to <file>.\n"
+"                 Valid types:  tin, eps, dem, obj\n"
+"-m <file>        Load the importance mask from <file>\n"
+"-s <file>        Execute preinsertion script from <file>\n"
+"\n";
+
+static void usage_error(char *msg = NULL)
+{
+    if( msg )
+	cerr << msg << endl;
+
+    cerr << "usage: terra <options> filename" << endl;
+    cerr << "       where <options> is some combination of: " << endl;
+    cerr << usage_string << endl;
+
+    exit(1);
+}
+
+void process_cmdline(int argc, char **argv)
+{
+    int c;
+
+    while( (c = getopt(argc, argv, options)) != EOF )
+    {
+	int errflg = False;
+	char *name;
+
+	switch( c )
+	{
+	case 'e':
+	    error_threshold = atof(optarg);
+	    cerr << "    Setting error threshold to " <<error_threshold <<endl;
+	    break;
+
+	case 'p':
+	    point_limit = atoi(optarg);
+	    cerr << "    Setting point limit to " << point_limit << endl;
+	    break;
+
+	case 'h':
+	    height_scale = atof(optarg);
+	    cerr << "    Setting height scaling factor to " << height_scale
+		 << endl;
+	    break;
+
+	case 'm':
+	    mask_filename = optarg;
+	    cerr << "    Input mask data from " << mask_filename << endl;
+	    break;
+
+	case 's':
+	    script_filename = optarg;
+	    cerr << "    Will execute script from " << script_filename << endl;
+	    break;
+
+	case 'o':
+	    output_filename = optarg;
+	    name = argv[optind++];
+
+	    if( !strcmp(name,"tin") )
+		output_format = TINfile;
+	    else if( !strcmp(name,"eps") )
+		output_format = EPSfile;
+	    else if( !strcmp(name,"dem") )
+		output_format = DEMfile;
+	    else if( !strcmp(name,"obj") )
+		output_format = OBJfile;
+	    else if( !strcmp(name,"rms") )
+		output_format = RMSfile;
+	    else
+		usage_error("Unknown output file type");
+
+	    cerr << "    Output file set to " << output_filename << endl;
+	    cerr << "    Output file type is " << name << endl;
+
+	    break;
+
+	default:
+	    errflg++;
+	    break;
+
+	}
+
+	if( errflg )
+	    usage_error();
+    }
+
+    ////////////////////////////////////////////////////////////////
+
+    if( optind==argc )
+	usage_error();
+
+    if( argv[optind][0]=='-' )
+	DEM = readPGM(cin);
+    else
+    {
+	ifstream in(argv[optind]);
+	DEM = readPGM(in);
+	in.close();
+    }
+
+    if( !DEM )
+    {
+	cerr << "Unable to load height field." << endl;
+	exit(1);
+    }
+
+    MASK = &default_mask;
+    if( mask_filename )
+    {
+	ifstream in(mask_filename);
+	MASK = readMask(in);
+	in.close();
+
+	if( !MASK )
+	{
+	    cerr << "Unable to load mask data ... reverting to identity mask.";
+	    cerr << endl;
+	    MASK = &default_mask;
+	}
+	else if( MASK->width != DEM->width  || MASK->height != DEM->height )
+	{
+	    cerr << "Mask resolution does not match DEM resolution." << endl;
+	    cerr << "     ... reverting to identity mask." << endl;
+	    MASK = &default_mask;
+	}
+    }
+    
+
+    mesh = new GreedySubdivision(DEM);
+
+
+    ////////////////////////////////////////////////////////////////
+
+    if( point_limit < 0 )
+	point_limit = DEM->width * DEM->height;
+
+
+    if( script_filename )
+    {
+	cerr << "Executing preinsertion script:" << endl;
+	ifstream in(script_filename);
+	scripted_preinsertion(in);
+	in.close();
+    }
+}

src/Prep/Terra/glHacks.cc

@@ -0,0 +1,28 @@
+#include "glHacks.h"
+
+void glGetViewport(int *x, int *y, int *w, int *h)
+{
+    GLint viewport[4];
+    glGetIntegerv(GL_VIEWPORT, viewport);
+
+    if( x ) *x = viewport[0];
+    if( y ) *y = viewport[1];
+    if( w ) *w = viewport[2];
+    if( h ) *h = viewport[3];
+}
+
+void glUnproject(int win_x, int win_y, int win_z,
+			double *x, double *y, double *z)
+{
+    GLdouble modelMatrix[16];
+    GLdouble projMatrix[16];
+    GLint viewport[4];
+
+    glGetDoublev(GL_PROJECTION_MATRIX, projMatrix);
+    glGetDoublev(GL_MODELVIEW_MATRIX, modelMatrix);
+    glGetIntegerv(GL_VIEWPORT, viewport);
+
+    gluUnProject(win_x,win_y,win_z,
+		 modelMatrix, projMatrix, viewport,
+		 x, y, z);
+}

src/Prep/Terra/glHacks.h

@@ -0,0 +1,189 @@
+#ifndef GLHACKS_INCLUDED // -*- C++ -*-
+#define GLHACKS_INCLUDED
+
+#include <GL/glx.h>
+#include <GL/gl.h>
+#include <GL/glu.h>
+
+
+/*************************************************************************
+ *
+ * Yes, here it is.  A bunch of nice overloaded routines to make it
+ * easier on the fingers to write OpenGL drawing code.
+ *
+ * The general stuff provided here is:
+ *
+ *      - glV
+ *      - glN
+ *      - glC
+ *      - glLoadMatrix/glGetMatrix
+ *
+ *************************************************************************/
+
+
+
+inline void glV(GLdouble x, GLdouble y) { glVertex2d(x,y); }
+inline void glV(GLfloat x, GLfloat y) { glVertex2f(x,y); }
+inline void glV(GLint x, GLint y) {  glVertex2i(x,y); }
+inline void glV(GLshort x, GLshort y) { glVertex2s(x,y); }
+
+inline void glV(GLdouble x, GLdouble y, GLdouble z) { glVertex3d(x,y,z); }
+inline void glV(GLfloat x, GLfloat y, GLfloat z) { glVertex3f(x,y,z); }
+inline void glV(GLint x, GLint y, GLint z) { glVertex3i(x,y,z); }
+inline void glV(GLshort x, GLshort y, GLshort z) { glVertex3s(x,y,z); }
+
+
+inline void glV(GLdouble x, GLdouble y, GLdouble z, GLdouble w)
+            { glVertex4d(x,y,z,w); }
+inline void glV(GLfloat x, GLfloat y, GLfloat z, GLfloat w)
+            { glVertex4f(x,y,z,w); }
+inline void glV(GLint x, GLint y, GLint z, GLint w) { glVertex4i(x,y,z,w); }
+inline void glV(GLshort x, GLshort y, GLshort z, GLshort w )
+            { glVertex4s(x,y,z,w); }
+
+inline void glV2(const GLdouble *v) { glVertex2dv(v); }
+inline void glV2(const GLfloat *v) { glVertex2fv(v); }
+inline void glV2(const GLint *v) { glVertex2iv(v); }
+inline void glV2(const GLshort *v) { glVertex2sv(v); }
+
+
+inline void glV3(const GLdouble *v) { glVertex3dv(v); }
+inline void glV3(const GLfloat *v) { glVertex3fv(v); }
+inline void glV3(const GLint *v) { glVertex3iv(v); }
+inline void glV3(const GLshort *v) { glVertex3sv(v); }
+
+
+inline void glV4(const GLdouble *v) { glVertex4dv(v); }
+inline void glV4(const GLfloat *v) { glVertex4fv(v); }
+inline void glV4(const GLint *v) { glVertex4iv(v); }
+inline void glV4(const GLshort *v) { glVertex4sv(v); }
+
+
+
+inline void glN(GLdouble x, GLdouble y, GLdouble z) { glNormal3d(x,y,z); }
+inline void glN(GLfloat x, GLfloat y, GLfloat z) { glNormal3f(x,y,z); }
+inline void glN(GLint x, GLint y, GLint z) { glNormal3i(x,y,z); }
+inline void glN(GLshort x, GLshort y, GLshort z) { glNormal3s(x,y,z); }
+
+inline void glN3(const GLdouble *v) { glNormal3dv(v); }
+inline void glN3(const GLfloat *v) { glNormal3fv(v); }
+inline void glN3(const GLint *v) { glNormal3iv(v); }
+inline void glN3(const GLshort *v) { glNormal3sv(v); }
+
+
+
+inline void glC(GLdouble x, GLdouble y, GLdouble z) { glColor3d(x,y,z); }
+inline void glC(GLfloat x, GLfloat y, GLfloat z) { glColor3f(x,y,z); }
+inline void glC(GLint x, GLint y, GLint z) { glColor3i(x,y,z); }
+inline void glC(GLshort x, GLshort y, GLshort z) { glColor3s(x,y,z); }
+
+inline void glC(GLdouble x, GLdouble y, GLdouble z, GLdouble w)
+            { glColor4d(x,y,z,w); }
+inline void glC(GLfloat x, GLfloat y, GLfloat z, GLfloat w)
+            { glColor4f(x,y,z,w); }
+inline void glC(GLint x, GLint y, GLint z, GLint w) { glColor4i(x,y,z,w); }
+inline void glC(GLshort x, GLshort y, GLshort z, GLshort w )
+            { glColor4s(x,y,z,w); }
+
+inline void glC3(const GLdouble *v) { glColor3dv(v); }
+inline void glC3(const GLfloat *v) { glColor3fv(v); }
+inline void glC3(const GLint *v) { glColor3iv(v); }
+inline void glC3(const GLshort *v) { glColor3sv(v); }
+
+inline void glC4(const GLdouble *v) { glColor4dv(v); }
+inline void glC4(const GLfloat *v) { glColor4fv(v); }
+inline void glC4(const GLint *v) { glColor4iv(v); }
+inline void glC4(const GLshort *v) { glColor4sv(v); }
+
+
+
+inline void glLoadMatrix(const GLdouble *m) { glLoadMatrixd(m); }
+inline void glLoadMatrix(const GLfloat *m) { glLoadMatrixf(m); }
+
+inline void glGetMatrix(GLdouble *m,GLenum src=GL_PROJECTION_MATRIX)
+{ glGetDoublev(src,m); }
+inline void glGetMatrix(GLfloat *m,GLenum src=GL_PROJECTION_MATRIX)
+{ glGetFloatv(src,m); }
+
+
+
+/*************************************************************************
+ *
+ * Here are some more generally useful things.
+ *
+ *************************************************************************/
+
+extern void glGetViewport(int *x, int *y, int *w, int *h);
+extern void glUnproject(int win_x, int win_y, int win_z,
+			double *x, double *y, double *z);
+
+typedef union {
+    unsigned int pixel;
+    struct { unsigned char r,g,b,a; } channel;
+} gfxPixel;
+
+
+inline gfxPixel *glSnapshot(int x, int y, int w, int h)
+{
+    gfxPixel *data = new gfxPixel[w*h];
+
+    glReadBuffer(GL_FRONT);
+    glReadPixels(x,y,w,h,
+		 GL_RGBA,
+		 GL_UNSIGNED_BYTE,
+		 data);
+
+    return data;
+}
+
+#ifdef __STDIO_H__
+inline void glToPPM(FILE *out)
+{
+    int x,y,w,h;
+    
+    glGetViewport(&x, &y, &w, &h);
+    gfxPixel *data = glSnapshot(x,y,w,h);
+
+    fprintf(out,"P6 %u %u 255\n",w,h);
+    int i,j;
+    for(j=h-1;j>=0;j--)
+	for(i=0;i<w;i++)
+	{
+	    int index = j*w + i;
+	    fputc(data[index].channel.r, out);
+	    fputc(data[index].channel.g, out);
+	    fputc(data[index].channel.b, out);
+	}
+
+    delete[] data;
+}
+#endif
+
+#ifdef IOSTREAMH
+inline ostream& operator<<(ostream& out,const gfxPixel& p)
+{
+    return out << p.channel.r << p.channel.g << p.channel.b;
+}
+
+inline void glToPPM(ostream& out)
+{
+    int x,y,w,h;
+    
+    glGetViewport(&x, &y, &w, &h);
+    gfxPixel *data = glSnapshot(x,y,w,h);
+
+    out << "P6 " << w <<" "<< h << " 255" << endl;
+    int i,j;
+    for(j=h-1;j>=0;j--)
+	for(i=0;i<w;i++)
+	{
+	    int index = j*w + i;
+	    out << data[index];
+	}
+
+    delete[] data;
+}
+#endif
+
+
+#endif

src/Prep/Terra/greedy.cc

@@ -0,0 +1,105 @@
+#include "terra.h"
+
+void scripted_preinsertion(istream& script)
+{
+    char op[4];
+    int x, y;
+
+    while( script.peek() != EOF )
+    {
+	script >> op >> x >> y;
+
+	switch( op[0] )
+	{
+	case 's':
+	    if( !mesh->is_used(x, y) )
+	    {
+		mesh->select(x, y);
+		mesh->is_used(x, y) = DATA_POINT_USED;
+	    }
+	    break;
+
+	case 'i':
+	    if( !mesh->is_used(x,y) )
+		mesh->is_used(x, y) = DATA_POINT_IGNORED;
+	    break;
+
+	case 'u':
+	    if( !mesh->is_used(x,y) )
+		mesh->is_used(x, y) = DATA_VALUE_UNKNOWN;
+	    break;
+
+	default:
+	    break;
+	}
+    }
+}
+
+void subsample_insertion(int target_width)
+{
+    int width = DEM->width;
+    int height = DEM->height;
+
+    // 'i' is the target width and 'j' is the target height
+
+    real i = (real)target_width;
+    real j = rint((i*height) / width);
+
+    real dx = (width-1)/(i-1);
+    real dy = (height-1)/(j-1);
+
+    real u, v;
+    int x, y;
+    for(u=0; u<i; u++)
+	for(v=0; v<j; v++)
+	{
+	    x =  (int)rint(u*dx);
+	    y =  (int)rint(v*dy);
+
+	    if( !mesh->is_used(x,y) )
+		mesh->select(x, y);
+	}
+}
+
+inline int goal_not_met()
+{
+    return mesh->maxError() > error_threshold &&
+	   mesh->pointCount() < point_limit;
+}
+
+static void announce_goal()
+{
+    cerr << "Goal conditions met:" << endl;
+    cerr << "     error=" << mesh->maxError()
+	 << " [thresh="<< error_threshold << "]" << endl;
+    cerr << "     points=" << mesh->pointCount()
+	 << " [limit=" << point_limit << "]" << endl;
+}
+
+void greedy_insertion()
+{
+
+    while( goal_not_met() )
+    {
+	if( !mesh->greedyInsert() )
+	    break;
+    }
+
+    announce_goal();
+}
+
+void display_greedy_insertion(void (*callback)())
+{
+
+    while( goal_not_met() )
+    {
+	if( !mesh->greedyInsert() )
+	{
+	    (*callback)();
+	    break;
+	}
+	(*callback)();
+    }
+
+    announce_goal();
+}

src/Prep/Terra/gui.cc

@@ -0,0 +1,456 @@
+#include <iostream.h>
+#include <fstream.h>
+
+#include <GL/glut.h>
+#include "glHacks.h"
+
+#include "terra.h"
+#include "gui.h"
+
+
+int mesh_view;
+int surf_view;
+
+int will_draw_dem = False;
+
+
+// Prototype for our hack below.
+//
+void xglutKeepAspect(float width, float height);
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// Mesh functions
+//
+////////////////////////////////////////////////////////////////////////
+#define MESH_MENU_EPS 1000
+#define MESH_MENU_TIN 1001
+#define MESH_MENU_OBJ 1002
+#define MESH_MENU_DEM 1003
+
+static void mesh_main_menu(int value)
+{
+    switch( value )
+    {
+    case MESH_MENU_EPS:
+	generate_output("out.eps", EPSfile);
+	break;
+
+    case MESH_MENU_TIN:
+	generate_output("out.tin", TINfile);
+	break;
+
+    case MESH_MENU_OBJ:
+	generate_output("out.obj", OBJfile);
+	break;
+
+    case MESH_MENU_DEM:
+	generate_output("out.pgm", DEMfile);
+
+    default:
+	break;
+    }
+
+}
+
+static void mesh_toggle_menu(int value)
+{
+    int *val = (int *)value;
+
+    *val = !(*val);
+
+    glutPostRedisplay();
+}
+
+
+
+static void create_mesh_menus()
+{
+    int toggle = glutCreateMenu(mesh_toggle_menu);
+
+    glutAddMenuEntry("Draw DEM data", (int)&will_draw_dem);
+
+
+    int main = glutCreateMenu(mesh_main_menu);
+    glutAddSubMenu("Toggle", toggle);
+
+    glutAddMenuEntry("Output Mesh EPS", MESH_MENU_EPS);
+    glutAddMenuEntry("Output Surface TIN", MESH_MENU_TIN);
+    glutAddMenuEntry("Output Surface OBJ", MESH_MENU_OBJ);
+    glutAddMenuEntry("Output Approximate DEM", MESH_MENU_DEM);
+
+    glutAttachMenu(GLUT_RIGHT_BUTTON);
+}
+
+
+
+static void draw_edge(Edge *e, void *)
+{
+    Vec2& org = e->Org();
+    Vec2& dest = e->Dest();
+    
+    glV(org[X], org[Y]);
+    glV(dest[X], dest[Y]);
+}
+
+static void draw_dem()
+{
+    Map *bits = DEM;
+
+    glPixelZoom((float)glutGet(GLUT_WINDOW_WIDTH) / (float)DEM->width ,
+		-(float)glutGet(GLUT_WINDOW_HEIGHT) / (float)DEM->height);
+    glRasterPos3f(0, 0, 0);
+
+
+    static GLenum type_map[] = {GL_UNSIGNED_BYTE,
+				GL_UNSIGNED_SHORT,
+				GL_UNSIGNED_INT };
+    glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
+
+    float scale = (float)((1<<bits->depth) - 1) / (DEM->max - DEM->min);
+    glPixelTransferf(GL_RED_SCALE, scale);
+    glPixelTransferf(GL_GREEN_SCALE, scale);
+    glPixelTransferf(GL_BLUE_SCALE, scale);
+
+    float bias  = -DEM->min/(DEM->max - DEM->min);
+    glPixelTransferf(GL_RED_BIAS, bias);
+    glPixelTransferf(GL_GREEN_BIAS, bias);
+    glPixelTransferf(GL_BLUE_BIAS, bias);
+
+    //!!
+    // WARNING:
+    // We are assuming that the map will never be a RealMap
+    //
+    GLenum type = type_map[((bits->depth)>>3)-1];
+    glDrawPixels(bits->width, bits->height,
+		 GL_LUMINANCE,
+		 type,
+		 bits->getBlock());
+}
+
+void mesh_display()
+{
+    glClear( GL_COLOR_BUFFER_BIT );
+
+    if( will_draw_dem )
+	draw_dem();
+
+
+    glBegin(GL_LINES);
+      glC(1.0, 0.15, 0.15);
+      mesh->overEdges(draw_edge);
+    glEnd();
+}
+
+static inline void redisplay_all(int other)
+{
+    glutPostRedisplay();
+    glutSetWindow(other);
+    glutPostRedisplay();
+}
+
+void mesh_mouse(int button, int state, int x, int y)
+{
+    if( state == GLUT_DOWN )
+    {
+	if( button == GLUT_LEFT_BUTTON )
+	{
+	    GLdouble u, v, z;
+	    int win_height = glutGet(GLUT_WINDOW_HEIGHT);
+
+	    glUnproject(x, win_height-y, 0, &u, &v, &z);
+	    if( u<0 ) u = 0;
+	    if( v<0 ) v = 0;
+	    if( u>(DEM->width-1) ) u = DEM->width-1;
+	    if( v>(DEM->height-1) ) v = DEM->height-1;
+
+	    Vec2 p(u,v);
+	    cout << "Inserting point: " << p << endl;
+
+	    mesh->insert(p);
+	    redisplay_all(surf_view);
+	}
+	else if( button == GLUT_MIDDLE_BUTTON )
+	{
+	    display_greedy_insertion(mesh_display);
+	    redisplay_all(surf_view);
+	}
+    }
+}
+
+
+
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// Surface functions
+//
+////////////////////////////////////////////////////////////////////////
+static void synthesize_normal(const Vec3& p1, const Vec3& p2,const Vec3& p3)
+{
+    // We explicitly declare these (rather than putting them in a
+    // Vector) so that they can be allocated into registers.
+    real v1x = p2[X]-p1[X],
+        v1y = p2[Y]-p1[Y],
+        v1z = p2[Z]-p1[Z];
+    real v2x = p3[X]-p2[X],
+         v2y = p3[Y]-p2[Y],
+         v2z = p3[Z]-p2[Z];
+
+    //
+    // NOTE: We do not unitize the normal vectors here because
+    //       they will be normalized in the graphics pipeline.
+    //       We have to enable GL_NORMALIZE because of the height
+    //       scaling we're doing.  So why bother normalizing now
+    //       when they're just going to get munged?
+    //
+
+    glN(v1y*v2z - v1z*v2y,
+	v1z*v2x - v1x*v2z,
+	v1x*v2y - v1y*v2x);
+}
+
+void render_face(Triangle& T, void *)
+{
+    const Vec2& p1 = T.point1();
+    const Vec2& p2 = T.point2();
+    const Vec2& p3 = T.point3();
+
+    Vec3 v1(p1, DEM->eval(p1[X],p1[Y]));
+    Vec3 v2(p2, DEM->eval(p2[X],p2[Y]));
+    Vec3 v3(p3, DEM->eval(p3[X],p3[Y]));
+
+    synthesize_normal(v1, v2, v3);
+
+    glV(v1[X], v1[Y], v1[Z]);
+    glV(v2[X], v2[Y], v2[Z]);
+    glV(v3[X], v3[Y], v3[Z]);
+}
+
+
+
+static long statex, statey;
+static enum {NoDrag, Spin, Dolly, Zoom } dragging = NoDrag;
+static float rotx = 0, roty = 0, rotz = 0;
+
+void surf_display()
+{
+    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+
+    glPushMatrix();
+      float w  = (float)DEM->width;
+      float h  = (float)DEM->height;
+      float cx = w / 2.0f;
+      float cy = h / 2.0f;
+      float cz = (float)DEM->min + (float)(DEM->max - DEM->min) / 2.0f;
+
+      glRotatef(rotx, 1.0, 0.0, 0.0);
+      glRotatef(roty, 0.0, 1.0, 0.0);
+      glRotatef(rotz, 0.0, 0.0, 1.0);
+      glScalef(1.0f, 1.0f, (float)height_scale);
+      glTranslatef(-cx, -cy, -cz);
+
+
+      glScalef(1.0f, -1.0f, 1.0f);
+      glTranslatef(0.0f, 1-h, 0.0f);
+
+      glBegin(GL_TRIANGLES);
+        mesh->overFaces(render_face);
+      glEnd();
+    glPopMatrix();
+
+    glutSwapBuffers();
+}
+
+
+void surf_mouse(int button, int state, int x, int y)
+{
+    if( state == GLUT_UP )
+    {
+	dragging = NoDrag;
+	return;
+    }
+
+    statex = x;
+    statey = y;
+
+    switch( button )
+    {
+    case GLUT_LEFT_BUTTON:
+	dragging = Spin;
+	break;
+
+    case GLUT_MIDDLE_BUTTON:
+	dragging = Dolly;
+	break;
+
+    case GLUT_RIGHT_BUTTON:
+	dragging = Zoom;
+	break;
+    }
+}
+
+
+void surf_motion(int x, int y)
+{
+    switch( dragging )
+    {
+    case Spin:
+        rotx += (float)(y - statey);
+        roty += (float)(x - statex);
+
+        statex = x;
+        statey = y;
+
+        glutPostRedisplay();
+	break;
+
+    case Dolly:
+	glMatrixMode(GL_PROJECTION);
+	glTranslatef((float)(x - statex), (float)(statey - y), 0.0f);
+        statex = x;
+        statey = y;
+	glMatrixMode(GL_MODELVIEW);
+
+        glutPostRedisplay();
+	break;
+
+    case Zoom:
+	glMatrixMode(GL_PROJECTION);
+	glTranslatef(0, 0, (float)(y - statey));
+        statey = y;
+	glMatrixMode(GL_MODELVIEW);
+
+        glutPostRedisplay();
+	break;
+
+    default:
+	break;
+    }
+}
+
+
+void init_surf_view(GreedySubdivision& mesh)
+{
+    glClearColor(0.3125, 0.3125, 1.0, 0.0);
+
+    glEnable(GL_DEPTH_TEST);
+    glEnable(GL_NORMALIZE);
+
+    //////////////////////
+    //
+    // Define viewing parameters
+    //
+    Map& map = mesh.getData();
+    real width = (real)map.width;
+    real height = (real)map.height;
+    real depth = map.max - map.min;
+
+    xglutKeepAspect(width, height);
+
+    glMatrixMode(GL_PROJECTION);
+    gluPerspective(45.0,
+		   width/height,
+		   0.1,
+		   10*depth);
+
+    gluLookAt(0, 0, depth,
+	      0, 0, 0,
+	      0, 1, 0);
+
+    glMatrixMode(GL_MODELVIEW);
+
+
+    //////////////////////
+    //
+    // Define lighting parameters
+    //
+    glShadeModel(GL_SMOOTH);
+
+    GLfloat mat_ambient[] = { 0.2, 0.2, 0.2, 1.0 };
+    GLfloat mat_diffuse[] = { 0.6, 0.6, 0.6, 1.0 };
+    // GLfloat mat_specular[] = { 0.5, 0.5, 0.5, 1.0 };
+
+    glMaterialfv(GL_FRONT, GL_AMBIENT, mat_ambient);
+    glMaterialfv(GL_FRONT, GL_DIFFUSE, mat_diffuse);
+    // glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
+
+    glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, 0);
+
+    GLfloat light_pos[] = { 1.0, 0.0, 1.0, 0.0 };
+    glLightfv(GL_LIGHT0, GL_POSITION, light_pos);
+
+    glEnable(GL_LIGHTING);
+    glEnable(GL_LIGHT0);
+}
+
+
+
+void gui_init()
+{
+    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
+    surf_view = glutCreateWindow("TERRA: Surface");
+    init_surf_view(*mesh);
+
+    glutDisplayFunc(surf_display);
+    glutMouseFunc(surf_mouse);
+    glutMotionFunc(surf_motion);
+
+    // ---------------------------------------------------------------------
+
+    glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
+    mesh_view = glutCreateWindow("TERRA: Mesh");
+
+    xglutKeepAspect(DEM->width, DEM->height);
+
+    create_mesh_menus();
+
+    glutDisplayFunc(mesh_display);
+    glutMouseFunc(mesh_mouse);
+
+    glMatrixMode(GL_PROJECTION);
+    gluOrtho2D(-1, DEM->width,DEM->height, -1);
+    // gluOrtho2D(0, DEM->width-1, DEM->height-1, 0);
+    glMatrixMode(GL_MODELVIEW);
+    glClearColor(0.0, 0.0, 0.0, 0.0);
+}
+
+void gui_interact()
+{
+    glutMainLoop();
+}
+
+
+
+////////////////////////////////////////////////////////////////////////
+//
+// WARNING: Entering hack zone!
+//
+// GLUT does not provide a function for setting the aspect
+// ratio of a window.  We want to do this.  Hence the following.
+//
+////////////////////////////////////////////////////////////////////////
+
+extern "C" {
+#include <GL/glutint.h>
+}
+
+void xglutKeepAspect(float width, float height)
+{
+    Window win;
+    XSizeHints hints;
+
+    if( __glutCurrentWindow )
+    {
+        win = __glutCurrentWindow->win;
+
+
+        hints.flags = PAspect;
+        hints.min_aspect.x = hints.max_aspect.x = (int)(1000*width);
+        hints.min_aspect.y = hints.max_aspect.y = (int)(1000*height);
+
+        XSetWMNormalHints(__glutDisplay, win, &hints);
+    }
+}

src/Prep/Terra/gui.h

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

src/Prep/Terra/output.cc

@@ -0,0 +1,192 @@
+#include "terra.h"
+#include <fstream.h>
+
+
+void generate_output(char *filename, FileFormat format)
+{
+    if( !filename )
+	filename = output_filename;
+    if( format==NULLfile )
+	format = output_format;
+
+    //
+    // If the user doesn't want output, don't do it.
+    if( !filename )
+	return;
+
+    ofstream out(filename);
+    switch( format )
+    {
+    case TINfile:
+	output_tin(out);
+	break;
+
+    case EPSfile:
+	output_eps(out);
+	break;
+
+    case DEMfile:
+	output_dem(out);
+	break;
+
+    case OBJfile:
+	output_obj(out);
+	break;
+
+    case RMSfile:
+	cout << mesh->pointCount() << "\t\t" << mesh->rmsError() << endl;
+	break;
+
+    default:
+	// Do nothing for unknown formats
+	break;
+    }
+    out.close();
+}
+
+
+
+
+static ostream *output_stream;
+
+
+static void obj_face(Triangle& T, void *closure)
+{
+    ostream& out = *output_stream;
+    array2<int>& vert_id = *(array2<int> *)closure;
+
+    const Vec2& p1 = T.point1();
+    const Vec2& p2 = T.point2();
+    const Vec2& p3 = T.point3();
+
+    out << "f ";
+    out << vert_id((int)p1[X], (int)p1[Y]) << " ";
+    out << vert_id((int)p2[X], (int)p2[Y]) << " ";
+    out << vert_id((int)p3[X], (int)p3[Y]) << endl;
+}
+
+static void obj_vertex(ostream& out, int x,int y)
+{
+    out << "v " << x <<" "<< y <<" "<< DEM->eval(x,y) << endl;
+}
+
+void output_obj(ostream& out)
+{
+    output_stream = &out;
+    int width = DEM->width;
+    int height = DEM->height;
+
+    array2<int> vert_id(width, height);
+    int index = 1;
+
+    for(int x=0;x<width;x++)
+	for(int y=0;y<height;y++)
+	{
+	    if( mesh->is_used(x,y) == DATA_POINT_USED )
+	    {
+		vert_id(x,y) = index++;
+		obj_vertex(out, x,y);
+	    }
+	    else
+		vert_id(x,y) = 0;
+	}
+
+    mesh->overFaces(obj_face, &vert_id);
+}
+
+
+
+
+static void tin_face(Triangle& T, void *)
+{
+    ostream& out = *output_stream;
+
+    const Vec2& p1 = T.point1();
+    const Vec2& p2 = T.point2();
+    const Vec2& p3 = T.point3();
+
+    out << "t ";
+
+    out << p1[X] <<" "<< p1[Y] <<" "<< DEM->eval(p1[X],p1[Y]) <<" ";
+    out << p2[X] <<" "<< p2[Y] <<" "<< DEM->eval(p2[X],p2[Y]) <<" ";
+    out << p3[X] <<" "<< p3[Y] <<" "<< DEM->eval(p3[X],p3[Y]) << endl;
+}
+
+void output_tin(ostream& out)
+{
+    output_stream = &out;
+
+    mesh->overFaces(tin_face);
+}
+
+
+static void eps_prologue(ostream& out)
+{
+    out << "%!PS-Adobe-2.0 EPSF-2.0" << endl;
+    out << "%%Creator: TERRA" << endl;
+    out << "%%BoundingBox: 0 0 " << DEM->width-1 <<" "<< DEM->height-1 << endl;
+    out << "%%EndComments" << endl;
+
+    out << endl;
+    out << "% -- The following transformation is a hack to change" << endl;
+    out << "% -- from image coordinates to page coordinates" << endl;
+    out << "% -- (in other words, flip things upside down)" << endl;
+    out << "1 -1 scale" << endl;
+    out << "0 -" << DEM->height-1 << " translate" << endl;
+    out << endl;
+
+    out << "/L {moveto lineto stroke} bind def" << endl;
+    out << "/C {newpath 2.5 0 360 arc closepath fill} bind def" << endl;
+}
+
+static void ps_edge(Edge *e, void *)
+{
+    ostream& out = *output_stream;
+
+    const Vec2& a = e->Org();
+    const Vec2& b = e->Dest();
+
+    out << a[X] << " " << a[Y] << " "
+	<< b[X] << " " << b[Y]
+	<< " L" << endl;
+}
+
+void output_eps(ostream& out)
+{
+    eps_prologue(out);
+
+    output_stream = &out;
+
+    out << ".2 setlinewidth" << endl;
+    mesh->overEdges(ps_edge);
+    out << "showpage" << endl;
+}
+
+void output_dem(ostream& out)
+{
+    cerr << ">> Writing approximate DEM output ... ";
+    cerr << "this may take a couple minutes." << endl;
+
+    int width = DEM->width;
+    int height = DEM->height;
+
+    out << "P2 " << width << " " << height << " " << DEM->max << endl;
+
+    streampos start = out.tellp();
+
+    for(int j=0;j<height;j++)
+	for(int i=0;i<width;i++)
+	{
+	    out << rint(mesh->eval(i,j));
+
+	    if( out.tellp() - start < 60 )
+		out << " ";
+	    else
+	    {
+		out << endl;
+		start = out.tellp();
+	    }
+	}
+    out << endl;
+    cerr << ">> Done." << endl;
+}

src/Prep/Terra/terra.cc

@@ -0,0 +1,10 @@
+#include "terra.h"
+
+main(int argc, char **argv)
+{
+    process_cmdline(argc, argv);
+
+    greedy_insertion();
+
+    generate_output();
+}

src/Prep/Terra/terra.h

@@ -0,0 +1,37 @@
+#ifndef TERRA_INCLUDED // -*- C++ -*-
+#define TERRA_INCLUDED
+
+#include "GreedyInsert.h"
+#include "Map.h"
+#include "Mask.h"
+
+
+extern GreedySubdivision *mesh;
+extern Map *DEM;
+extern ImportMask *MASK;
+
+extern real error_threshold;
+extern int point_limit;
+extern real height_scale;
+enum FileFormat {NULLfile, TINfile, EPSfile, DEMfile, OBJfile, RMSfile};
+extern FileFormat output_format;
+extern char *output_filename;
+extern char *script_filename;
+
+extern int goal_not_met();
+extern void greedy_insertion();
+extern void display_greedy_insertion(void (*callback)());
+extern void scripted_preinsertion(istream&);
+extern void subsample_insertion(int target_width);
+
+extern void generate_output(char *filename=NULL,
+			    FileFormat format=NULLfile);
+extern void output_tin(ostream&);
+extern void output_eps(ostream&);
+extern void output_obj(ostream&);
+extern void output_dem(ostream&);
+
+extern void process_cmdline(int argc, char **argv);
+
+
+#endif

src/Prep/Terra/version.h

@@ -0,0 +1,8 @@
+#ifndef VERSION_INCLUDED
+#define VERSION_INCLUDED
+
+#define terra_version 0.7
+#define terra_version_string "0.7"
+#define terra_tag_string "RELEASE_0_7"
+
+#endif

src/Prep/Terra/xterra.cc

@@ -0,0 +1,18 @@
+#include <GL/glut.h>
+
+#include "terra.h"
+#include "gui.h"
+
+
+
+main(int argc, char **argv)
+{
+    glutInit(&argc, argv);
+    process_cmdline(argc, argv);
+
+
+    gui_init();
+
+
+    gui_interact();
+}