This is a first pass at removing the nurbs++ dependency from TerraGear (in

favor of newmat11 which is much simpler, and seems to compile well on modern
OS's.)  I need to do some further testing of genapts and until then, don't
assume the new mechanism is working perfectly.

README.newmat

@@ -0,0 +1,25 @@
+You *must* have the newmat package installed on your system to build
+the airport generator portion of TerraGear.
+
+You can get the latest version of newmat from:
+
+    http://www.robertnz.net/nm_intro.htm
+
+This package is distributed with several make files.  Read the docs
+(online) and follow those instructions to build the package.  For
+instance, to build on a typical "gnu" system:
+
+- extract the package into a new subdirectory
+- run "make -f nm_gnu.mak"
+- copy the newly created libnewmat.a to /usr/local/lib
+- make a new directory called /usr/local/include/newmat
+- copy all the *.h files to /usr/local/include/newmat
+
+Licensing terms:
+
+The newmat11 documentation states: "There are no restrictions on the
+use of newmat except that I take no liability for any problems that
+may arise from this use."
+
+
+

README.nurbs++

@@ -1,46 +0,0 @@
-You *must* have the libnurbs++ package installed on your system to
-build portions of TerraGear".
-
-You can get the latest version of libnurbs++ from:
-
-    http://libnurbs.sourceforge.net/
-
-IMPORTANT!  If you are using gcc-3.4.x (or newer) you will need to apply
-a patch to your nurbs++ tree to successfully compile.  Look for the file
-called PATCH-nurbs++-gcc-3.4.x in the current directory.  Now cd to the
-top level of your fresh nurbs++-3.0.11 source tree and run:
-
-    patch -p0 < /path/to/patch/file
-
-Good luck!
-
-
-More information:
-
-  Non-Uniform Rational B-Splines (NURBS) curves and surface are
-  parametric functions which can represent any type of curves or
-  surfaces. This C++ library hides the basic mathematics of
-  NURBS. This allows the user to focus on the more challenging parts
-  of their projects. The library also offers a lot of features to help
-  generate NURBS from data points.
-
-  The NURBS++ package includes a matrix library, an image manipulation
-  library, a numerical library and a NURBS library. They can all be
-  used on their own but they are all developped to support my NURBS
-  needs.
-
-  NURBS++ is being revived under sourceforge as the nurbs++ project
-  (it is under the libnurbs directory for obvious reasons). Please use
-  the sourceforge tools to send bug reports, feature requests, etc.
-
-  This library is now used in the Mind's Eyes project and also in the
-  Innovation3D project. Both projects aim at creating a free 3D
-  modeller for UNIX machines. If you know other projects which use the
-  library please let me know. It helps justify the effort I put into
-  it.
-
-  The library is under the GNU Library Public Licence term. This means
-  it's free. You can copy it, modify it and even enjoy it. You do have
-  to add my name as the author of the work, after all this is still
-  copyrighted material. If the above license is bad for you, talk to
-  me. I'm sure we can work something out.

configure.ac

@@ -247,9 +247,9 @@ if test "x$ac_cv_header_plib_pu_h" != "xyes"; then
     exit
 fi
 
-dnl Check for "libnurbs++" without which we cannot go on
+dnl Check for "libnewmat" without which we cannot build airport surfaces
-AC_CHECK_HEADER(nurbs++/nurbsS.h)
+AC_CHECK_HEADER(newmat/newmat.h)
-AM_CONDITIONAL(HAVE_NURBS, test "x$ac_cv_header_nurbspp_nurbsS_h" = "xyes" )
+AM_CONDITIONAL(HAVE_NEWMAT, test "x$ac_cv_header_newmat_newmat_h" = "xyes" )
 AC_LANG_POP
 
 AC_CHECK_HEADER(gts.h)
@@ -442,11 +442,11 @@ else
     echo "Debug messages: yes"
 fi
 
-if test "x$ac_cv_header_nurbspp_nurbsS_h" != "xyes"; then
+if test "x$ac_cv_header_newmat_newmat_h" != "xyes"; then
     echo
-    echo "You must have the nurbs++ library installed on your system to build"
+    echo "You must have the newmat library installed on your system to build"
-    echo "the GenAirport utility. This program will not be build now."
+    echo "the GenAirport utility. This program will not be built."
     echo
-    echo "Please see README.nurbs++ for more details."
+    echo "Please see README.newmat for more details."
     echo
 fi

src/Airports/GenAirports/Makefile.am

@@ -23,7 +23,7 @@
 #---------------------------------------------------------------------------
 
 
-bin_PROGRAMS = genapts testnurbs
+bin_PROGRAMS = genapts
 
 genapts_SOURCES = \
 	apt_surface.hxx apt_surface.cxx \
@@ -53,13 +53,10 @@ genapts_LDADD = \
 	$(top_builddir)/src/Lib/TriangleJRS/libTriangleJRS.a \
 	-lsgbucket -lsgdebug -lsgio -lsgmath -lsgmisc -lsgstructure -lsgxml \
 	-lgenpolyclip \
-	-lnurbsd -lmatrixI -lmatrixN -lmatrix \
+	-lnewmat \
 	-lz \
 	$(base_LIBS)
 
-testnurbs_SOURCES = testnurbs.cxx
-testnurbs_LDADD = -lnurbsf -lmatrixI -lmatrixN -lmatrix
-
 INCLUDES = \
 	-I$(top_srcdir) \
 	-I$(top_srcdir)/src \

src/Airports/GenAirports/apt_surface.cxx

@@ -22,9 +22,14 @@
 // $Id$
 //
 
+#include <simgear/compiler.h>
 
-#include <nurbs++/nurbsS.h>
+// libnewmat includes and defines
-#include <nurbs++/nurbsSub.h>
+#define WANT_STREAM		// include.h will get stream fns
+#define WANT_MATH		// include.h will get math fns
+				// newmatap.h will get include.h
+#include <newmat/newmatap.h>	// need matrix applications
+#include <newmat/newmatio.h>	// need matrix output routines
 
 #include <simgear/constants.h>
 #include <simgear/math/sg_geodesy.hxx>
@@ -37,17 +42,15 @@
 #include "global.hxx"
 #include "apt_surface.hxx"
 
-SG_USING_NAMESPACE( PLib );
 
-
+static bool limit_slope( SimpleMatrix *Pts, int i1, int j1, int i2, int j2,
-static bool limit_slope( Matrix_Point3Dd &Pts, int i1, int j1, int i2, int j2,
                          double average_elev_m )
 {
     bool slope_error = false;
 
-    Point3Dd p1, p2;
+    Point3D p1, p2;
-    p1 = Pts(i1,j1);
+    p1 = Pts->element(i1,j1);
-    p2 = Pts(i2,j2);
+    p2 = Pts->element(i2,j2);
 
     double az1, az2, dist;
     double slope;
@@ -73,19 +76,19 @@ static bool limit_slope( Matrix_Point3Dd &Pts, int i1, int j1, int i2, int j2,
         if ( e1 > e2 ) {
             // cout << "  p1 error larger" << endl;
             if ( slope > 0 ) {
-                p1.z() = p2.z() - (dist * slope_max);
+	      p1.setz( p2.z() - (dist * slope_max) );
             } else {
-                p1.z() = p2.z() + (dist * slope_max);
+	      p1.setz( p2.z() + (dist * slope_max) );
             }
-            Pts(i1,j1) = p1;
+            Pts->set(i1, j1, p1);
         } else {
             // cout << "  p2 error larger" << endl;
             if ( slope > 0 ) {
-                p2.z() = p1.z() + (dist * slope_max);
+	      p2.setz( p1.z() + (dist * slope_max) );
             } else {
-                p2.z() = p1.z() - (dist * slope_max);
+	      p2.setz( p1.z() - (dist * slope_max) );
             }
-            Pts(i2,j2) = p2;
+            Pts->set(i2, j2, p2);
         }
         // cout << "   z1 = " << p1.z() << "  z2 = " << p2.z() << endl;
     }
@@ -126,18 +129,12 @@ TGAptSurface::TGAptSurface( const string& path,
     int xdivs = (int)(x_m / coarse_grid) + 1;
     int ydivs = (int)(y_m / coarse_grid) + 1;
 
-#if defined( _NURBS_GLOBAL_INTERP )
+    // set an arbitrary minumum number of divisions to keep things
-    if ( xdivs < 3 ) { xdivs = 3; }
+    // interesting
-    if ( ydivs < 3 ) { ydivs = 3; }
-#elif defined( _NURBS_LEAST_SQUARES )
-    // Minimum divs appears to need to be at least 5 before the
-    // leastsquares nurbs surface approximation stops crashing.
     if ( xdivs < 6 ) { xdivs = 6; }
     if ( ydivs < 6 ) { ydivs = 6; }
-#else
-# error "Need to define _NURBS_GLOBAL_INTER or _NURBS_LEAST_SQUARES"
-#endif
     SG_LOG(SG_GENERAL, SG_INFO, "  M(" << ydivs << "," << xdivs << ")");
+
     double dlon = x_deg / xdivs;
     double dlat = y_deg / ydivs;
 
@@ -147,14 +144,15 @@ TGAptSurface::TGAptSurface( const string& path,
     // Build the extra res input grid (shifted SW by half (dlon,dlat)
     // with an added major row column on the NE sides.)
     int mult = 10;
-    Matrix_Point3Dd dPts( (ydivs + 1) * mult + 1, (xdivs + 1) * mult + 1 );
+    SimpleMatrix dPts( (ydivs + 1) * mult + 1, (xdivs + 1) * mult + 1 );
     for ( int j = 0; j < dPts.cols(); ++j ) {
         for ( int i = 0; i < dPts.rows(); ++i ) {
-            dPts(i,j) = Point3Dd( min_deg.lon() - dlon_h
+	  dPts.set(i, j, Point3D( min_deg.lon() - dlon_h
-                                    + j * (dlon / (double)mult),
+				    + j * (dlon / (double)mult),
-                                  min_deg.lat() - dlat_h
+				   min_deg.lat() - dlat_h
-                                    + i * (dlat / (double)mult),
+				    + i * (dlat / (double)mult),
-                                  -9999 );
+				   -9999 )
+		   );
         }
     }
 
@@ -183,21 +181,21 @@ TGAptSurface::TGAptSurface( const string& path,
 #endif
 
     // Build the normal res input grid from the double res version
-    Matrix_Point3Dd Pts(ydivs + 1, xdivs + 1);
+    Pts = new SimpleMatrix(ydivs + 1, xdivs + 1);
     double ave_divider = (mult+1) * (mult+1);
-    for ( int j = 0; j < Pts.cols(); ++j ) {
+    for ( int j = 0; j < Pts->cols(); ++j ) {
-        for ( int i = 0; i < Pts.rows(); ++i ) {
+        for ( int i = 0; i < Pts->rows(); ++i ) {
             SG_LOG(SG_GENERAL, SG_DEBUG, i << "," << j);
             double accum = 0.0;
             double lon_accum = 0.0;
             double lat_accum = 0.0;
             for ( int jj = 0; jj <= mult; ++jj ) {
                 for ( int ii = 0; ii <= mult; ++ii ) {
-                    double value = dPts(mult*i + ii, mult*j + jj).z();
+                    double value = dPts.element(mult*i + ii, mult*j + jj).z();
                     SG_LOG( SG_GENERAL, SG_DEBUG, "value = " << value );
                     accum += value;
-                    lon_accum += dPts(mult*i + ii, mult*j + jj).x();
+                    lon_accum += dPts.element(mult*i + ii, mult*j + jj).x();
-                    lat_accum += dPts(mult*i + ii, mult*j + jj).y();
+                    lat_accum += dPts.element(mult*i + ii, mult*j + jj).y();
                 }
             }
             double val_ave = accum / ave_divider;
@@ -205,9 +203,10 @@ TGAptSurface::TGAptSurface( const string& path,
             double lat_ave = lat_accum / ave_divider;
 
             SG_LOG( SG_GENERAL, SG_DEBUG, "  val_ave = " << val_ave );
-            Pts(i,j) = Point3Dd( min_deg.lon() + j * dlon,
+            Pts->set(i, j, Point3D( min_deg.lon() + j * dlon,
-                                 min_deg.lat() + i * dlat,
+				   min_deg.lat() + i * dlat,
-                                 val_ave );
+				   val_ave )
+		    );
             SG_LOG( SG_GENERAL, SG_DEBUG, "lon_ave = " << lon_ave
                     << "  lat_ave = " << lat_ave );
             SG_LOG( SG_GENERAL, SG_DEBUG, "lon = " << min_deg.lon() + j * dlon
@@ -216,8 +215,8 @@ TGAptSurface::TGAptSurface( const string& path,
     }
 
 #ifdef DEBUG
-    for ( int j = 0; j < Pts.cols(); ++j ) {
+    for ( int j = 0; j < Pts->cols(); ++j ) {
-        for ( int i = 0; i < Pts.rows(); ++i ) {
+        for ( int i = 0; i < Pts->rows(); ++i ) {
             printf("%.5f %.5f %.1f\n", Pts(i,j).x(), Pts(i,j).y(),
                    Pts(i,j).z() );
         }
@@ -229,8 +228,8 @@ TGAptSurface::TGAptSurface( const string& path,
         SG_LOG( SG_GENERAL, SG_DEBUG, "start of slope processing pass" );
         slope_error = false;
         // Add some "slope" sanity to the resulting surface grid points
-        for ( int j = 0; j < Pts.cols() - 1; ++j ) {
+        for ( int j = 0; j < Pts->cols() - 1; ++j ) {
-            for ( int i = 0; i < Pts.rows() - 1; ++i ) {
+            for ( int i = 0; i < Pts->rows() - 1; ++i ) {
                 if ( limit_slope( Pts, i, j, i+1, j, average_elev_m ) ) {
                     slope_error = true;
                 }
@@ -245,52 +244,36 @@ TGAptSurface::TGAptSurface( const string& path,
     }
 
 #ifdef DEBUG
-    for ( int j = 0; j < Pts.cols(); ++j ) {
+    for ( int j = 0; j < Pts->cols(); ++j ) {
-        for ( int i = 0; i < Pts.rows(); ++i ) {
+        for ( int i = 0; i < Pts->rows(); ++i ) {
             printf("%.5f %.5f %.1f\n", Pts(i,j).x(), Pts(i,j).y(),
                    Pts(i,j).z() );
         }
     }
 #endif
 
-    // Create the nurbs surface
+    // Create the fitted surface
-
+    SG_LOG(SG_GENERAL, SG_DEBUG, "ready to create fitted surface");
-    SG_LOG(SG_GENERAL, SG_DEBUG, "ready to create nurbs surface");
+    fit();
-    apt_surf = new PlNurbsSurfaced;
-#if defined( _NURBS_GLOBAL_INTERP )
-    apt_surf->globalInterp( Pts, 3, 3);
-#elif defined( _NURBS_LEAST_SQUARES )
-    cout << "Col = " << Pts.cols() << " Rows = " << Pts.rows() << endl;
-    int nU = Pts.rows() / 2; if ( nU < 4 ) { nU = 4; }
-    int nV = Pts.cols() / 2; if ( nV < 4 ) { nV = 4; }
-    cout << "nU = " << nU << " nV = " << nV << endl;
-    apt_surf->leastSquares( Pts, 3, 3, nU, nV );
 
     // sanity check: I'm finding that leastSquares() can produce nan
     // surfaces.  We test for this and fall back to globalInterp() if
     // the least squares fails.
-    double result =  query_solver( (min_deg.lon() + max_deg.lon()) / 2.0,
+    double result =  query( (min_deg.lon() + max_deg.lon()) / 2.0,
-                                   (min_deg.lat() + max_deg.lat()) / 2.0 );
+			    (min_deg.lat() + max_deg.lat()) / 2.0 );
-    Point3Dd p = apt_surf->pointAt( 0.5, 0.5 );
+    if ( result > -9000.0 ) {
-
-    if ( (result > -9000.0) && (p.z() <= 0.0 || p.z() >= 0.0) ) {
         // ok, a valid number
     } else {
-        // no, sorry, a nan is not <= 0.0 or >= 0.0
         SG_LOG(SG_GENERAL, SG_WARN,
-               "leastSquares() nurbs interpolation failed!!!");
+               "leastSquares() fit seemed to fail!!!");
         char command[256];
         sprintf( command,
                  "echo least squares nurbs interpolation failed, using globalInterp() >> last_apt" );
         system( command );
 
-        // we could fall back to globalInterp() rather than aborting
+	// abort and force developer to debug for now
-        // if we wanted to ...
+	exit(-1);
-        apt_surf->globalInterp( Pts, 3, 3);
     }
-#else
-# error "Need to define _NURBS_GLOBAL_INTER or _NURBS_LEAST_SQUARES"
-#endif
     SG_LOG(SG_GENERAL, SG_DEBUG, "  successful.");
 
 #ifdef DEBUG
@@ -315,7 +298,110 @@ TGAptSurface::TGAptSurface( const string& path,
 
 
 TGAptSurface::~TGAptSurface() {
-    delete apt_surf;
+    delete Pts;
+}
+
+
+static ColumnVector qr_method( Real* y,
+			       Real* t1, Real* t2, Real* t3, Real* t4,
+			       Real* t5, Real* t6, Real* t7, Real* t8,
+			       int nobs, int npred )
+{
+  cout << "QR triangularisation" << endl;;
+
+  // QR triangularisation method
+ 
+  // load data - 1s into col 1 of matrix
+  int npred1 = npred+1;
+  Matrix X(nobs,npred1); ColumnVector Y(nobs);
+  X.column(1) = 1.0;
+  X.column(2) << t1;
+  X.column(3) << t2;
+  X.column(4) << t3;
+  X.column(5) << t4;
+  X.column(6) << t5;
+  X.column(7) << t6;
+  X.column(8) << t7;
+  X.column(9) << t8;
+  Y << y;
+
+  // do Householder triangularisation
+  // no need to deal with constant term separately
+  Matrix X1 = X;                 // Want copy of matrix
+  ColumnVector Y1 = Y;
+  UpperTriangularMatrix U; ColumnVector M;
+  QRZ(X1, U); QRZ(X1, Y1, M);    // Y1 now contains resids
+  ColumnVector A = U.i() * M;
+  ColumnVector Fitted = X * A;
+  Real ResVar = sum_square(Y1) / (nobs-npred1);
+
+  // get variances of estimates
+  U = U.i(); DiagonalMatrix D; D << U * U.t();
+
+  // Get diagonals of Hat matrix
+  DiagonalMatrix Hat;  Hat << X1 * X1.t();
+
+  // print out answers
+  cout << "\nEstimates and their standard errors\n\n";
+  ColumnVector SE(npred1);
+  for (int i=1; i<=npred1; i++) SE(i) = sqrt(D(i)*ResVar);
+  cout << setw(11) << setprecision(5) << (A | SE) << endl;
+  cout << "\nObservations, fitted value, residual value, hat value\n";
+  cout << setw(9) << setprecision(3) << 
+    (X.columns(2,3) | Y | Fitted | Y1 | Hat.as_column());
+  cout << "\n\n";
+
+  return A;
+}
+
+
+// Use a linear least squares method to fit a 3d polynomial to the
+// sampled surface data
+void TGAptSurface::fit() {
+
+  // the fit function is: f(x,y) = a*x + b*x^2 + c*y + d*y^2 + e*x*y +
+  // f*x*y^2 + g*x^2*y + h*x^2*y^2
+  int nobs = Pts->cols() * Pts->rows();	// number of observations
+  int npred = 8;		        // number of predictor values
+
+  Real z[nobs];
+  Real t1[nobs];
+  Real t2[nobs];
+  Real t3[nobs];
+  Real t4[nobs];
+  Real t5[nobs];
+  Real t6[nobs];
+  Real t7[nobs];
+  Real t8[nobs];
+
+  // generate the required fit data
+  for ( int j = 0; j < Pts->cols(); j++ ) {
+    for ( int i = 0; i <= Pts->rows(); i++ ) {
+      Point3D p = Pts->element( i, j );
+      int index = ( j * Pts->rows() ) + i;
+      Real xi = p.x();
+      Real yi = p.y();
+      z[index] = p.z();
+      t1[index] = xi;
+      t2[index] = xi*xi;
+      t3[index] = yi;
+      t4[index] = yi*yi;
+      t5[index] = xi*yi;
+      t6[index] = xi*yi*yi;
+      t7[index] = xi*xi*yi;
+      t8[index] = xi*xi*yi*yi;
+    }
+  }
+
+  // we want to find the values of a,b,c to give the best
+  // fit
+
+  Try {
+    surface_coefficients
+      = qr_method(z, t1, t2, t3, t4, t5, t6, t7, t8, nobs, npred);
+  }
+  CatchAll { cout << BaseException::what(); }
+
 }
 
 
@@ -325,198 +411,21 @@ double TGAptSurface::query( double lon_deg, double lat_deg ) {
     if ( lon_deg < min_deg.lon() || lon_deg > max_deg.lon() ||
          lat_deg < min_deg.lat() || lat_deg > max_deg.lat() )
     {
-        SG_LOG(SG_GENERAL, SG_WARN, "Warning: query out of bounds for NURBS surface!");
+        SG_LOG(SG_GENERAL, SG_WARN,
+	       "Warning: query out of bounds for fitted surface!");
         return -9999.0;
     }
 
-    double lat_range = max_deg.lat() - min_deg.lat();
+    // compute the function with solved coefficients
-    double lon_range = max_deg.lon() - min_deg.lon();
 
-    // convert lon,lat to nurbs space (NOTE: that this is a dumb
+    // the fit function is: f(x,y) = a*x + b*x^2 + c*y + d*y^2 + e*x*y +
-    // approximation that assumes that nurbs surface space exactly
+    //                               f*x*y^2 + g*x^2*y + h*x^2*y^2
-    // corresponds to real world x,y space which it doesn't, but maybe
-    // the error is small enough that it doesn't matter?)
-    double u = (lat_deg - min_deg.lat()) / lat_range;
-    double v = (lon_deg - min_deg.lon()) / lon_range;
 
-    Point3Dd p = apt_surf->pointAt( u, v );
+    double x = lon_deg;
+    double y = lat_deg;
+    ColumnVector A = surface_coefficients;
+    double result = A(0)*x + A(1)*x*x + A(2)*y + A(3)*y*y + A(4)*x*y
+      + A(5)*x*y*y + A(6)*x*x*y + A(7)*x*x*y*y;
 
-    // double az1, az2, dist;
+    return result;
-    // geo_inverse_wgs_84( 0, lat_deg, lon_deg, p.y(), p.x(),
-    //                     &az1, &az2, &dist );
-    // cout << "query distance error = " << dist << endl;
-    
-    return p.z();
-}
-
-
-// Query the elevation of a point, return -9999 if out of range
-double TGAptSurface::query_solver( double lon_deg, double lat_deg ) {
-    // sanity check
-    if ( lon_deg < min_deg.lon() || lon_deg > max_deg.lon() ||
-         lat_deg < min_deg.lat() || lat_deg > max_deg.lat() )
-    {
-        SG_LOG(SG_GENERAL, SG_WARN, "Warning: query out of bounds for NURBS surface!");
-        return -9999.0;
-    }
-
-    double lat_range = max_deg.lat() - min_deg.lat();
-    double lon_range = max_deg.lon() - min_deg.lon();
-
-    // convert lon,lat to nurbs space
-    double u = (lat_deg - min_deg.lat()) / lat_range;
-    double v = (lon_deg - min_deg.lon()) / lon_range;
-
-    // cout << "running quick solver ..." << endl;
-
-    int count;
-    double dx = 1000;
-    double dy = 1000;
-    Point3Dd p;
-
-    // cout << " solving for u,v simultaneously" << endl;
-    count = 0;
-    while ( count < 0 /* disable fast solver for now */ ) {
-        if ( u < 0.0 || u > 1.0 || v < 0.0 || v > 1.0 ) {
-            cout << "oops, something goofed in the solver" << endl;
-            exit(-1);
-        }
-
-        p = apt_surf->pointAt( u, v );
-        // cout << "  querying for " << u << ", " << v << " = " << p.z()
-        //      << endl;
-        // cout << "    found " << p.x() << ", " << p.y() << " = " << p.z()
-        //      << endl;
-        dx = lon_deg - p.x();
-        dy = lat_deg - p.y();
-        // cout << "      dx = " << dx << " dy = " << dy << endl;
-
-        if ( (fabs(dx) < nurbs_eps) && (fabs(dy) < nurbs_eps) ) {
-            return p.z();
-        } else {
-            u = u + (dy/lat_range);
-            v = v + (dx/lon_range);
-            if ( u < 0.0 ) { u = 0.0; }
-            if ( u > 1.0 ) { u = 1.0; }
-            if ( v < 0.0 ) { v = 0.0; }
-            if ( v > 1.0 ) { v = 1.0; }
-        }
-
-        ++count;
-    }
-
-    // cout << "quick query solver failed..." << endl;
-
-    // failed to converge with fast approach, try a binary partition
-    // scheme, however we have to solve each axis independently
-    // because when we move in one axis we may change our position in
-    // the other axis.
-
-    double min_u = 0.0; 
-    double max_u = 1.0;
-    double min_v = 0.0; 
-    double max_v = 1.0;
-
-    u = (max_u + min_u) / 2.0;
-    v = (max_v + min_v) / 2.0;
-
-    dx = 1000.0;
-    dy = 1000.0;
-
-    int gcount = 0;
-
-    while ( true ) {
-        // cout << "solving for u" << endl;
-        
-        min_u = 0.0; 
-        max_u = 1.0;
-        count = 0;
-        while ( true ) {
-            p = apt_surf->pointAt( u, v );
-            // cout << "  binary querying for " << u << ", " << v << " = "
-            //      << p.z() << endl;
-            // cout << "    found " << p.x() << ", " << p.y() << " = " << p.z()
-            //      << endl;
-            dx = lon_deg - p.x();
-            dy = lat_deg - p.y();
-            // cout << "      dx = " << dx << " dy = " << dy << " z = " << p.z() << endl;
-
-            // double az1, az2, dist;
-            // geo_inverse_wgs_84( 0, lat_deg, lon_deg, p.y(), p.x(),
-            //                     &az1, &az2, &dist );
-            // cout << "      query distance error = " << dist << endl;
-
-            if ( fabs(dy) < nurbs_eps ) {
-                break;
-            } else {
-                if ( dy >= nurbs_eps ) {
-                    min_u = u;
-                } else if ( dy <= nurbs_eps ) {
-                    max_u = u;
-                }
-            }
-
-            u = (max_u + min_u) / 2.0;
-
-            if ( count > 100 ) {
-                // solver failed
-                cout << "binary solver failed..." << endl;
-                return -9999.0;
-            }
-
-            ++count;
-        }
-
-        // cout << "solving for v" << endl;
-
-        min_v = 0.0; 
-        max_v = 1.0;
-        count = 0;
-        while ( true ) {
-            p = apt_surf->pointAt( u, v );
-            // cout << "  binary querying for " << u << ", " << v << " = "
-            //      << p.z() << endl;
-            // cout << "    found " << p.x() << ", " << p.y() << " = " << p.z()
-            //      << endl;
-            dx = lon_deg - p.x();
-            dy = lat_deg - p.y();
-            // cout << "      dx = " << dx << " dy = " << dy << " z = " << p.z() << endl;
-
-            // double az1, az2, dist;
-            // geo_inverse_wgs_84( 0, lat_deg, lon_deg, p.y(), p.x(),
-            //                     &az1, &az2, &dist );
-            // cout << "      query distance error = " << dist << endl;
-            if ( fabs(dx) < nurbs_eps ) {
-                break;
-            } else {
-                if ( dx >= nurbs_eps ) {
-                    min_v = v;
-                } else if ( dx <= nurbs_eps ) {
-                    max_v = v;
-                }
-            }
-            v = (max_v + min_v) / 2.0;
-
-            if ( count > 100 ) {
-                // solver failed
-                cout << "binary solver failed..." << endl;
-                return -9999.0;
-            }
-
-            ++count;
-        }
-
-        // cout << "query count = " << gcount << "  dist = " << dx << ", "
-        //      << dy << endl;
-
-        if ( (fabs(dx) < nurbs_eps) && (fabs(dy) < nurbs_eps) ) {
-            // double az1, az2, dist;
-            // geo_inverse_wgs_84( 0, lat_deg, lon_deg, p.y(), p.x(),
-            //                     &az1, &az2, &dist );
-            // cout << "        final query distance error = " << dist << endl;
-            return p.z();
-        }
-    }
-
-    return p.z();
 }

src/Airports/GenAirports/apt_surface.hxx

@@ -29,23 +29,64 @@
 
 #include <string>
 
-#include <nurbs++/nurbsS.h>
+// libnewmat includes and defines
+#define WANT_STREAM		// include.h will get stream fns
+#define WANT_MATH		// include.h will get math fns
+				// newmatap.h will get include.h
+#include <newmat/newmatap.h>	// need matrix applications
+#include <newmat/newmatio.h>	// need matrix output routines
 
 #include <simgear/math/point3d.hxx>
 
 
+/***
+ * A dirt simple matrix class for our convenience based on top of Point3D
+ */
+
+#define SIMPLE_MATRIX_MAX_SIZE 200
+
+class SimpleMatrix {
+
+private:
+
+  int _rows;
+  int _cols;
+  point_list m;
+
+public:
+
+  inline SimpleMatrix( int columns, int rows ) {
+    _cols = columns;
+    _rows = rows;
+    m.resize( _cols * _rows );
+  }
+
+  inline Point3D element( int col, int row ) {
+    int index = ( col * _rows ) + row;
+    return m[index];
+  }
+
+  inline void set( int col, int row, Point3D p ) {
+    int index = ( col * _rows ) + row;
+    m[index] = p;
+  }
+
+  inline int cols() const { return _cols; }
+  inline int rows() const { return _rows; }
+};
+
+
 /***
  * Note of explanation.  When a TGAptSurface instance is created, you
  * must specify a min and max lon/lat containing the entire airport
  * area.  The class will divide up that area into a reasonably sized
  * regular grid.  It will then look up the elevation of each point on
- * the grid from the DEM/Array data.  Finally it will build a nurbs
+ * the grid from the DEM/Array data.  Finally it will fit do a linear
- * surface from this grid.  Each vertex of the actual airport model is
+ * least squares polygonal surface approximation from this grid.  Each
- * drapped over this nurbs surface rather than over the underlying
+ * vertex of the actual airport model is drapped over this fitted
- * terrain data.  This provides a) smoothing of noisy terrain data, b)
+ * surface rather than over the underlying terrain data.  This
- * natural rises and dips in the airport surface, c) chance to impress
+ * provides a) smoothing of noisy terrain data and b) natural rises
- * colleages by using something or other nurbs surfaces -- or whatever
+ * and dips in the airport surface.
- * they are called. :-)
  */
 
 class TGAptSurface {
@@ -53,7 +94,8 @@ class TGAptSurface {
 private:
 
     // The actual nurbs surface approximation for the airport
-    PlNurbsSurfaced *apt_surf;
+    SimpleMatrix *Pts;
+    ColumnVector surface_coefficients;
 
     Point3D min_deg, max_deg;
 
@@ -71,18 +113,15 @@ public:
     // Destructor
     ~TGAptSurface();
 
+    // Use a linear least squares method to fit a 3d polynomial to the
+    // sampled surface data
+    void fit();
+
     // Query the elevation of a point, return -9999 if out of range.
     // This routine makes a simplistic assumption that X,Y space is
     // proportional to u,v space on the nurbs surface which it isn't.
     double query( double lon_deg, double lat_deg );
 
-    // Query the elevation of a point, return -9999 if out of range.
-    // This routine incorporates a complex solver that attempts to
-    // find the exact u,v coordinates correspondinging to the
-    // requested lat, lon.  But there seems to be problems with my
-    // solver routine that I haven't tracked down yet so I'm testing
-    // the _dumb version above.
-    double query_solver( double lon_deg, double lat_deg );
 };
 
 

src/Airports/GenAirports/build.cxx

@@ -191,7 +191,7 @@ static point_list calc_elevations( TGAptSurface &surf,
 {
     point_list result = geod_nodes;
     for ( unsigned int i = 0; i < result.size(); ++i ) {
-        double elev = surf.query_solver( result[i].lon(), result[i].lat() );
+        double elev = surf.query( result[i].lon(), result[i].lat() );
         result[i].setelev( elev + offset );
     }
 

src/Airports/GenAirports/elevations.cxx

@@ -23,8 +23,12 @@
 //
 
 
-#include <nurbs++/nurbsS.h>
+// libnewmat includes and defines
-#include <nurbs++/nurbsSub.h>
+#define WANT_STREAM		// include.h will get stream fns
+#define WANT_MATH		// include.h will get math fns
+				// newmatap.h will get include.h
+#include <newmat/newmatap.h>	// need matrix applications
+#include <newmat/newmatio.h>	// need matrix output routines
 
 #include <simgear/constants.h>
 #include <simgear/math/sg_geodesy.hxx>
@@ -36,8 +40,6 @@
 #include "global.hxx"
 #include "apt_surface.hxx"
 
-SG_USING_NAMESPACE( PLib );
-
 
 // lookup node elevations for each point in the point_list.  Returns
 // average of all points.  Doesn't modify the original list.
@@ -64,7 +66,7 @@ double tgAverageElevation( const string &root, const string_list elev_src,
 
     while ( !done ) {
 	// find first node with -9999 elevation
-        Point3D first;
+        Point3D first(0.0);
         bool found_one = false;
         for ( i = 0; i < points.size(); ++i ) {
             if ( points[i].z() < -9000.0 && !found_one ) {
@@ -142,7 +144,8 @@ double tgAverageElevation( const string &root, const string_list elev_src,
 // matrix.  Returns average of all points.
 
 void tgCalcElevations( const string &root, const string_list elev_src,
-                       Matrix_Point3Dd &Pts, const double average ) {
+                       SimpleMatrix &Pts, const double average )
+{
     bool done = false;
     int i, j;
     TGArray array;
@@ -154,19 +157,20 @@ void tgCalcElevations( const string &root, const string_list elev_src,
 
     // set all elevations to -9999
     for ( j = 0; j < Pts.cols(); ++j ) {
-        for ( i = 0; i < Pts.rows(); ++i ) {
+      for ( i = 0; i < Pts.rows(); ++i ) {
-            Point3Dd p = Pts(i,j);
+	Point3D p = Pts.element(i, j);
-            p.z() = -9999.0;
+	p.setz( -9999.0 );
-        }
+	Pts.set(i, j, p);
+      }
     }
 
     while ( !done ) {
 	// find first node with -9999 elevation
-        Point3Dd first;
+        Point3D first(0.0);
         bool found_one = false;
         for ( j = 0; j < Pts.cols(); ++j ) {
             for ( i = 0; i < Pts.rows(); ++i ) {
-                Point3Dd p = Pts(i,j);
+                Point3D p = Pts.element(i,j);
                 if ( p.z() < -9000.0 && !found_one ) {
                     first = p;
                     found_one = true;
@@ -206,17 +210,17 @@ void tgCalcElevations( const string &root, const string_list elev_src,
 	    done = true;
             for ( j = 0; j < Pts.cols(); ++j ) {
                 for ( i = 0; i < Pts.rows(); ++i ) {
-                    Point3Dd p = Pts(i,j);
+                    Point3D p = Pts.element(i,j);
                     if ( p.z() < -9000.0 ) {
                         done = false;
                         elev = array.altitude_from_grid( p.x() * 3600.0,
                                                          p.y() * 3600.0 );
                         if ( elev > -9000 ) {
-                            p.z() = elev;
+			  p.setz( elev );
-                            Pts(i,j) = p;
+			  Pts.set(i, j, p);
-                            // cout << "interpolating for " << p << endl;
+			  // cout << "interpolating for " << p << endl;
-                            // cout << p.x() << " " << p.y() << " " << p.z()
+			  // cout << p.x() << " " << p.y() << " " << p.z()
-                            //      << endl;
+			  //      << endl;
                         }
                     }
                 }
@@ -235,7 +239,7 @@ void tgCalcElevations( const string &root, const string_list elev_src,
     int count = 0;
     for ( j = 0; j < Pts.cols(); ++j ) {
         for ( i = 0; i < Pts.rows(); ++i ) {
-            Point3Dd p = Pts(i,j);
+            Point3D p = Pts.element(i,j);
             total += p.z();
             count++;
         }
@@ -248,7 +252,7 @@ void tgCalcElevations( const string &root, const string_list elev_src,
 
 // clamp all elevations to the specified range
 
-void tgClampElevations( Matrix_Point3Dd &Pts,
+void tgClampElevations( SimpleMatrix &Pts,
                         double center_m, double max_clamp_m )
 {
     int i, j;
@@ -257,18 +261,18 @@ void tgClampElevations( Matrix_Point3Dd &Pts,
     // +/-max_m of the center_m elevation.
     for ( j = 0; j < Pts.cols(); ++j ) {
         for ( i = 0; i < Pts.rows(); ++i ) {
-            Point3Dd p = Pts(i,j);
+            Point3D p = Pts.element(i,j);
             if ( p.z() < center_m - max_clamp_m ) {
                 SG_LOG(SG_GENERAL, SG_DEBUG, "   clamping " << p.z()
                        << " to " << center_m - max_clamp_m );
-                p.z() = center_m - max_clamp_m;
+                p.setz( center_m - max_clamp_m );
-                Pts(i,j) = p;
+                Pts.set(i, j, p);
             }
             if ( p.z() > center_m + max_clamp_m ) {
                 SG_LOG(SG_GENERAL, SG_DEBUG, "   clamping " << p.z()
                        << " to " << center_m + max_clamp_m );
-                p.z() = center_m + max_clamp_m;
+                p.setz( center_m + max_clamp_m );
-                Pts(i,j) = p;
+                Pts.set(i, j, p);
             }
         }
     }

src/Airports/GenAirports/elevations.hxx

@@ -23,8 +23,12 @@
 //
 
 
-#include <nurbs++/nurbsS.h>
+// libnewmat includes and defines
-#include <nurbs++/nurbsSub.h>
+#define WANT_STREAM		// include.h will get stream fns
+#define WANT_MATH		// include.h will get math fns
+				// newmatap.h will get include.h
+#include <newmat/newmatap.h>	// need matrix applications
+#include <newmat/newmatio.h>	// need matrix output routines
 
 #include <simgear/constants.h>
 #include <simgear/math/sg_geodesy.hxx>
@@ -36,8 +40,6 @@
 #include "global.hxx"
 #include "apt_surface.hxx"
 
-SG_USING_NAMESPACE( PLib );
-
 
 // lookup node elevations for each point in the point_list.  Returns
 // average of all points.  Doesn't modify the original list.
@@ -47,8 +49,8 @@ double tgAverageElevation( const string &root, const string_list elev_src,
 // lookup node elevations for each point in the specified nurbs++
 // matrix.
 void tgCalcElevations( const string &root, const string_list elev_src,
-                       Matrix_Point3Dd &Pts, double average );
+                       SimpleMatrix &Pts, double average );
 
 // clamp all elevations to the specified range
-void tgClampElevations( Matrix_Point3Dd &Pts,
+void tgClampElevations( SimpleMatrix &Pts,
-                        double center_m, double max_clamp_m );
+			double center_m, double max_clamp_m );

src/Airports/GenAirports/global.hxx

@@ -42,8 +42,4 @@ const double slope_eps = 0.00001;
 // nurbs query/search epsilon
 const double nurbs_eps = 0.0000001;
 
-// Define only one of the following
-// #define _NURBS_GLOBAL_APPROX 1
-#define _NURBS_LEAST_SQUARES 1
-
 #endif // _GEN_AIRPORT_GLOBAL_HXX

src/Airports/GenAirports/testnurbs.cxx

@@ -1,98 +0,0 @@
-#include <math.h>
-
-#include <nurbs++/nurbsS.h>
-#include <nurbs++/nurbsSub.h>
-
-
-int main(){
-    Matrix_Point3Df Pts(4,5) ;
-    int i,j ;
-
-    using namespace PLib ; 
-
-
-    for(i=0;i<Pts.rows();++i){
-        for(j=0;j<Pts.cols();++j){
-            Pts(i,j) = Point3Df(i,j,j) ;
-        }
-    }
-
-    PlNurbsSurfacef surf ;
-
-    surf.globalInterp(Pts,3,3) ;
-
-    cerr << "Point at (0.5,0.5) = " << surf(0.5,0.5) << endl ;
-    cerr << "\t=  " << surf.pointAt(0.5,0.5) << endl ;
-    cerr << "\t normal = " << surf.normal(0.5,0.5) << endl ;
-
-    for(i=0;i<Pts.rows();++i)
-        for(j=0;j<Pts.cols();++j)
-            Pts(i,j) = Point3Df(i*20,j*20,(1+cos(i+j))*20) ;
-
-    surf.globalInterp(Pts,3,3) ;
-
-    cerr << "Point at (0.5,0.5) = " << surf(0.5,0.5) << endl ;
-    cerr << "\t=  " << surf.pointAt(0.5,0.5) << endl ;
-    cerr << "\t normal = " << surf.normal(0.5,0.5) << endl ;
-
-    surf.write("tnurbsS.ns") ;
-    surf.read("tnurbsS.ns");
-    surf.writeVRML("tnurbsS.wrl") ;
-    surf.writePS("tnurbsS.ps",5,5,Point3Df(10,10,10),Point3Df(0,0,0)) ;
-  
-
-    PLib::NurbsSubSurface<float> sub(surf) ;
-    sub.drawSubdivisionVRML("tnurbsSb.wrl",0.5);
-
-    surf.writeVRML("tnurbsS2.wrl",Color(255,255,0),float(0.1)) ;
-
-
-    NurbsCurvef curve ;
-    curve.makeCircle(Point3Df(0.5,0.5,0),0.3);
-
-
-    Vector_Point3Df pnt(100) ;
-
-    for(i=0;i<100;++i){
-        Point3Df param = curve.pointAt(float(i)/99.0) ;
-        pnt[i] = surf.pointAt(param.x(),param.y());
-    }
-
-    NurbsCurvef curve2 ;
-
-    curve2.globalInterp(pnt,3) ; 
-
-    ofstream fout ;
-
-    fout.open("tnurbsST.wrl");
-    surf.writeVRML(fout,Color(255,255,255)) ;
-
-    curve2.writeVRML(fout,0.1,30,Color(255,0,0),6,60) ;
-
-    curve.writeVRML(fout,0.1,30,Color(0,255,0),6,60) ;
-
-    curve.degreeElevate(7) ;
-
-    for(i=0;i<curve.ctrlPnts().n();++i){
-        HPoint3Df p = surf(curve.ctrlPnts(i).x(),curve.ctrlPnts(i).y()) ;
-        curve.modCP(i,p) ;
-    }
-
-    curve.writeVRML(fout,0.1,30,Color(0,0,255),6,60);
-
-    fout.close();
-
-
-    //surf.writePOVRAY(0.1f,"tnurbsS.pov",Color(255,255,0),Point3Df(0,1,0),Point3Df(0,0,1)) ;
-
-    cerr << "Testing the normal computation at the corner points\n" ;
-    cerr << "normal(0,0) = " << surf.normal(0,0).unitLength() << endl ;
-    cerr << "normal(0,1) = " << surf.normal(0,1).unitLength() << endl ;
-    cerr << "normal(1,0) = " << surf.normal(1,0).unitLength() << endl ;
-    cerr << "normal(1,1) = " << surf.normal(1,1).unitLength() << endl ;
-
-    surf.makeTorus(Point3Df(0,0,0),4,1) ;
-    surf.writeVRML("torus.wrl");
-
-    return 0 ;
-}

src/Airports/Makefile.am

@@ -1,4 +1,4 @@
-if HAVE_NURBS
+if HAVE_NEWMAT
 GEN_AIRPORTS_DIR = GenAirports
 else
 GEN_AIRPORTS_DIR =