From bf7d86732d478104b5e22a26739d29b0e512d23a Mon Sep 17 00:00:00 2001 From: Christian Schmitt Date: Wed, 31 Oct 2012 13:13:28 +0100 Subject: [PATCH] Remove the newmat dependency and bring in TNT templates as a substitute --- CMakeLists.txt | 1 - CMakeModules/FindNewMat11.cmake | 25 - README.newmat | 25 - projects/VC100/newmat11/newmat11.sln | 26 - projects/VC100/newmat11/newmat11.vcproj | 412 ---------------- projects/VC100/newmat11/newmat11.vcxproj | 188 -------- .../VC100/newmat11/newmat11.vcxproj.filters | 105 ---- src/Airports/CMakeLists.txt | 4 +- src/Airports/GenAirports850/CMakeLists.txt | 3 - src/Airports/GenAirports850/airport.cxx | 2 +- src/Airports/GenAirports850/apt_surface.cxx | 203 ++------ src/Airports/GenAirports850/apt_surface.hxx | 29 +- src/Airports/GenAirports850/elevations.cxx | 21 +- src/Airports/GenAirports850/elevations.hxx | 14 +- src/Lib/Geometry/TNT/jama_qr.h | 326 +++++++++++++ src/Lib/Geometry/TNT/tnt_array1d.h | 301 ++++++++++++ src/Lib/Geometry/TNT/tnt_array2d.h | 451 ++++++++++++++++++ src/Lib/Geometry/TNT/tnt_i_refvec.h | 243 ++++++++++ src/Lib/Geometry/TNT/tnt_math_utils.h | 37 ++ 19 files changed, 1430 insertions(+), 986 deletions(-) delete mode 100644 CMakeModules/FindNewMat11.cmake delete mode 100644 README.newmat delete mode 100644 projects/VC100/newmat11/newmat11.sln delete mode 100644 projects/VC100/newmat11/newmat11.vcproj delete mode 100644 projects/VC100/newmat11/newmat11.vcxproj delete mode 100644 projects/VC100/newmat11/newmat11.vcxproj.filters create mode 100644 src/Lib/Geometry/TNT/jama_qr.h create mode 100644 src/Lib/Geometry/TNT/tnt_array1d.h create mode 100644 src/Lib/Geometry/TNT/tnt_array2d.h create mode 100644 src/Lib/Geometry/TNT/tnt_i_refvec.h create mode 100644 src/Lib/Geometry/TNT/tnt_math_utils.h diff --git a/CMakeLists.txt b/CMakeLists.txt index dde89670..fe346ecf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -71,7 +71,6 @@ find_package(Threads REQUIRED) find_package(SimGear 2.9.0 REQUIRED) find_package(GDAL REQUIRED) find_package(TIFF REQUIRED) # needed for SRTM -find_package(NewMat11 REQUIRED) find_library(POCO_FOUNDATION PocoFoundation REQUIRED) find_library(POCO_NET PocoNet REQUIRED) diff --git a/CMakeModules/FindNewMat11.cmake b/CMakeModules/FindNewMat11.cmake deleted file mode 100644 index ed7788a4..00000000 --- a/CMakeModules/FindNewMat11.cmake +++ /dev/null @@ -1,25 +0,0 @@ - -FIND_PATH(NEWMAT_INCLUDE_DIR newmat/newmat.h - HINTS $ENV{NEWMAT_DIR} - PATH_SUFFIXES include - PATHS - /usr/local - /usr - /opt -) - - -FIND_LIBRARY(NEWMAT_LIBRARY - NAMES newmat - HINTS - $ENV{NEWMAT_DIR} - PATH_SUFFIXES lib64 lib libs64 libs libs/Win32 libs/Win64 - PATHS - /usr/local - /usr - /opt -) - -include(FindPackageHandleStandardArgs) -FIND_PACKAGE_HANDLE_STANDARD_ARGS(newmat DEFAULT_MSG - NEWMAT_LIBRARY NEWMAT_INCLUDE_DIR) diff --git a/README.newmat b/README.newmat deleted file mode 100644 index 95570b6b..00000000 --- a/README.newmat +++ /dev/null @@ -1,25 +0,0 @@ -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." - - - diff --git a/projects/VC100/newmat11/newmat11.sln b/projects/VC100/newmat11/newmat11.sln deleted file mode 100644 index 09355c32..00000000 --- a/projects/VC100/newmat11/newmat11.sln +++ /dev/null @@ -1,26 +0,0 @@ - -Microsoft Visual Studio Solution File, Format Version 11.00 -# Visual Studio 2010 -Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "newmat11", "newmat11.vcxproj", "{BAAE2014-7AB0-4AE7-BA70-AFD26E7EB15A}" -EndProject -Global - GlobalSection(SolutionConfigurationPlatforms) = preSolution - Debug|Win32 = Debug|Win32 - Debug|x64 = Debug|x64 - Release|Win32 = Release|Win32 - Release|x64 = Release|x64 - EndGlobalSection - GlobalSection(ProjectConfigurationPlatforms) = postSolution - {BAAE2014-7AB0-4AE7-BA70-AFD26E7EB15A}.Debug|Win32.ActiveCfg = Debug|Win32 - {BAAE2014-7AB0-4AE7-BA70-AFD26E7EB15A}.Debug|Win32.Build.0 = Debug|Win32 - {BAAE2014-7AB0-4AE7-BA70-AFD26E7EB15A}.Debug|x64.ActiveCfg = Debug|x64 - {BAAE2014-7AB0-4AE7-BA70-AFD26E7EB15A}.Debug|x64.Build.0 = Debug|x64 - {BAAE2014-7AB0-4AE7-BA70-AFD26E7EB15A}.Release|Win32.ActiveCfg = Release|Win32 - {BAAE2014-7AB0-4AE7-BA70-AFD26E7EB15A}.Release|Win32.Build.0 = Release|Win32 - {BAAE2014-7AB0-4AE7-BA70-AFD26E7EB15A}.Release|x64.ActiveCfg = Release|x64 - {BAAE2014-7AB0-4AE7-BA70-AFD26E7EB15A}.Release|x64.Build.0 = Release|x64 - EndGlobalSection - GlobalSection(SolutionProperties) = preSolution - HideSolutionNode = FALSE - EndGlobalSection -EndGlobal diff --git a/projects/VC100/newmat11/newmat11.vcproj b/projects/VC100/newmat11/newmat11.vcproj deleted file mode 100644 index 8ab824b3..00000000 --- a/projects/VC100/newmat11/newmat11.vcproj +++ /dev/null @@ -1,412 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/projects/VC100/newmat11/newmat11.vcxproj b/projects/VC100/newmat11/newmat11.vcxproj deleted file mode 100644 index a6e30d14..00000000 --- a/projects/VC100/newmat11/newmat11.vcxproj +++ /dev/null @@ -1,188 +0,0 @@ - - - - - Debug - Win32 - - - Debug - x64 - - - Release - Win32 - - - Release - x64 - - - - {BAAE2014-7AB0-4AE7-BA70-AFD26E7EB15A} - Win32Proj - - - - StaticLibrary - MultiByte - - - StaticLibrary - MultiByte - - - StaticLibrary - MultiByte - - - StaticLibrary - MultiByte - - - - - - - - - - - - - - - - - - - - - - - <_ProjectFileVersion>10.0.40219.1 - $(Platform)\$(Configuration)\ - $(Platform)\$(Configuration)\ - $(Platform)\$(Configuration)\ - $(Platform)\$(Configuration)\ - $(Platform)\$(Configuration)\ - $(Platform)\$(Configuration)\ - $(Platform)\$(Configuration)\ - $(Platform)\$(Configuration)\ - AllRules.ruleset - - - AllRules.ruleset - - - AllRules.ruleset - - - AllRules.ruleset - - - $(ProjectName)d - $(ProjectName)d - - - - Disabled - WIN32;_DEBUG;_LIB;_CRT_SECURE_NO_WARNINGS;%(PreprocessorDefinitions) - true - EnableFastChecks - MultiThreadedDebugDLL - - - Level3 - EditAndContinue - - - $(OutDir)newmat11d.lib - - - - - X64 - - - Disabled - WIN32;_DEBUG;_LIB;_CRT_SECURE_NO_WARNINGS;%(PreprocessorDefinitions) - true - EnableFastChecks - MultiThreadedDebugDLL - - - Level3 - ProgramDatabase - - - $(OutDir)newmat11d.lib - - - - - WIN32;NDEBUG;_LIB;_CRT_SECURE_NO_WARNINGS;%(PreprocessorDefinitions) - MultiThreadedDLL - true - - - Level3 - ProgramDatabase - - - $(OutDir)newmat11.lib - - - - - X64 - - - WIN32;NDEBUG;_LIB;_CRT_SECURE_NO_WARNINGS;%(PreprocessorDefinitions) - MultiThreadedDLL - true - - - Level3 - ProgramDatabase - - - $(OutDir)newmat11.lib - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ No newline at end of file diff --git a/projects/VC100/newmat11/newmat11.vcxproj.filters b/projects/VC100/newmat11/newmat11.vcxproj.filters deleted file mode 100644 index a7454e11..00000000 --- a/projects/VC100/newmat11/newmat11.vcxproj.filters +++ /dev/null @@ -1,105 +0,0 @@ - - - - - {4FC737F1-C7A5-4376-A066-2A32D752A2FF} - cpp;c;cxx;def;odl;idl;hpj;bat;asm;asmx - - - {93995380-89BD-4b04-88EB-625FBE52EBFB} - h;hpp;hxx;hm;inl;inc;xsd - - - {67DA6AB6-F800-4c08-8B7A-83BB121AAD01} - rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx - - - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - Source Files - - - - - Header Files - - - Header Files - - - Header Files - - - Header Files - - - Header Files - - - \ No newline at end of file diff --git a/src/Airports/CMakeLists.txt b/src/Airports/CMakeLists.txt index 47ac41e1..c7e2bd1c 100644 --- a/src/Airports/CMakeLists.txt +++ b/src/Airports/CMakeLists.txt @@ -1,6 +1,4 @@ include_directories(${PROJECT_SOURCE_DIR}/src/Lib) include_directories(${PROJECT_SOURCE_DIR}/src/BuildTiles) -if (NEWMAT_FOUND) -add_subdirectory(GenAirports850) -endif (NEWMAT_FOUND) +add_subdirectory(GenAirports850) \ No newline at end of file diff --git a/src/Airports/GenAirports850/CMakeLists.txt b/src/Airports/GenAirports850/CMakeLists.txt index 3e241ee2..576d8778 100644 --- a/src/Airports/GenAirports850/CMakeLists.txt +++ b/src/Airports/GenAirports850/CMakeLists.txt @@ -1,5 +1,3 @@ -include_directories(${NEWMAT_INCLUDE_DIR}) - add_executable(genapts850 airport.hxx airport.cxx apt_math.hxx apt_math.cxx @@ -30,7 +28,6 @@ target_link_libraries(genapts850 ${GDAL_LIBRARY} ${SIMGEAR_CORE_LIBRARIES} ${SIMGEAR_CORE_LIBRARY_DEPENDENCIES} - ${NEWMAT_LIBRARY} ${RT_LIBRARY}) install(TARGETS genapts850 RUNTIME DESTINATION bin) diff --git a/src/Airports/GenAirports850/airport.cxx b/src/Airports/GenAirports850/airport.cxx index 40baa3ad..0a00e5cf 100644 --- a/src/Airports/GenAirports850/airport.cxx +++ b/src/Airports/GenAirports850/airport.cxx @@ -1586,7 +1586,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src ) // add light points superpoly_list tmp_light_list; tmp_light_list.clear(); - typedef map < string, double, less > elev_map_type; + typedef std::map < string, double, std::less > elev_map_type; typedef elev_map_type::const_iterator const_elev_map_iterator; elev_map_type elevation_map; diff --git a/src/Airports/GenAirports850/apt_surface.cxx b/src/Airports/GenAirports850/apt_surface.cxx index 1bc1f4f2..795550a6 100644 --- a/src/Airports/GenAirports850/apt_surface.cxx +++ b/src/Airports/GenAirports850/apt_surface.cxx @@ -19,29 +19,19 @@ // along with this program; if not, write to the Free Software // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. // -// $Id: apt_surface.cxx,v 1.31 2005-12-19 16:51:25 curt Exp $ -// #ifdef HAVE_CONFIG_H # include #endif +#include + #include - -// 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 // need matrix applications -#include // need matrix output routines - #include #include #include #include -#include - #include "elevations.hxx" #include "global.hxx" #include "apt_surface.hxx" @@ -103,7 +93,7 @@ static bool limit_slope( SimpleMatrix *Pts, int i1, int j1, int i2, int j2, // Constructor, specify min and max coordinates of desired area in // lon/lat degrees -TGAptSurface::TGAptSurface( const string& path, +TGAptSurface::TGAptSurface( const std::string& path, const string_list& elev_src, Point3D _min_deg, Point3D _max_deg, double _average_elev_m ) @@ -300,145 +290,60 @@ TGAptSurface::~TGAptSurface() { } -static ColumnVector qr_method( Real* y, - Real* t1, Real* t2, Real* t3, Real* t4, - Real* t5, Real* t6, Real* t7, Real* t8, - Real* t9, Real* t10, Real* t11, Real* t12, - Real* t13, Real* t14, Real* t15, - int nobs, int npred ) -{ - SG_LOG(SG_GENERAL, SG_DEBUG, "QR triangularisation" ); - - // QR triangularisation method - - // load data - 1s into col 1 of matrix - int npred1 = npred+1; - SG_LOG(SG_GENERAL, SG_DEBUG, "nobs = " << nobs << " npred1 = " << npred1 ); - - 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; - X.column(10) << t9; - X.column(11) << t10; - X.column(12) << t11; - X.column(13) << t12; - X.column(14) << t13; - X.column(15) << t14; - X.column(16) << t15; - 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; - - // get variances of estimates - U = U.i(); DiagonalMatrix D; D << U * U.t(); - - // Get diagonals of Hat matrix - DiagonalMatrix Hat; Hat << X1 * X1.t(); - -#ifdef DEBUG - cout << "A vector = " << A << endl; - cout << "A rows = " << A.nrows() << endl; - - // 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,4) | Y | Fitted | Y1 | Hat.as_column()); - cout << "\n\n"; -#endif - - 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) = A1*x + A2*x*y + A3*y + - // A4*x*x + A5+x*x*y + A6*x*x*y*y + A7*y*y + A8*x*y*y + - // A9*x*x*x + A10*x*x*x*y + A11*x*x*x*y*y + A12*x*x*x*y*y*y + - // A13*y*y*y + A14*x*y*y*y + A15*x*x*y*y*y + // the fit function is: + // f(x,y) = A1*x + A2*x*y + A3*y + + // A4*x*x + A5+x*x*y + A6*x*x*y*y + A7*y*y + A8*x*y*y + + // A9*x*x*x + A10*x*x*x*y + A11*x*x*x*y*y + A12*x*x*x*y*y*y + + // A13*y*y*y + A14*x*y*y*y + A15*x*x*y*y*y - int nobs = Pts->cols() * Pts->rows(); // number of observations - int npred = 15; // number of predictor values A[n] + int nobs = Pts->cols() * Pts->rows(); // number of observations - vector tz(nobs); - vector t1(nobs); - vector t2(nobs); - vector t3(nobs); - vector t4(nobs); - vector t5(nobs); - vector t6(nobs); - vector t7(nobs); - vector t8(nobs); - vector t9(nobs); - vector t10(nobs); - vector t11(nobs); - vector t12(nobs); - vector t13(nobs); - vector t14(nobs); - vector t15(nobs); + SG_LOG(SG_GENERAL, SG_DEBUG, "QR triangularisation" ); - // generate the required fit data - for ( int j = 0; j < Pts->rows(); j++ ) { - for ( int i = 0; i < Pts->cols(); i++ ) { - Point3D p = Pts->element( i, j ); - int index = ( j * Pts->cols() ) + i; - Real x = p.x() - offset.x(); - Real y = p.y() - offset.y(); - Real z = p.z() - offset.z(); - //cout << "pt = " << x << "," << y << "," << z << endl; - tz[index] = z; - t1[index] = x; - t2[index] = x*y; - t3[index] = y; - t4[index] = x*x; - t5[index] = x*x*y; - t6[index] = x*x*y*y; - t7[index] = y*y; - t8[index] = x*y*y; - t9[index] = x*x*x; - t10[index] = x*x*x*y; - t11[index] = x*x*x*y*y; - t12[index] = x*x*x*y*y*y; - t13[index] = y*y*y; - t14[index] = x*y*y*y; - t15[index] = x*x*y*y*y; + // Create an array (matrix) with 16 columns (predictor values) A[n] + TNT::Array2D mat(nobs,16); + + // put all elevation values into a second array + TNT::Array1D zmat(nobs); + + // generate the required fit data + for ( int j = 0; j < Pts->rows(); j++ ) { + for ( int i = 0; i < Pts->cols(); i++ ) { + Point3D p = Pts->element( i, j ); + int index = ( j * Pts->cols() ) + i; + double x = p.x() - offset.x(); + double y = p.y() - offset.y(); + double z = p.z() - offset.z(); + + zmat[index] = z; + + mat[index][0] = 1.0; + mat[index][1] = x; + mat[index][2] = x*y; + mat[index][3] = y; + mat[index][4] = x*x; + mat[index][5] = x*x*y; + mat[index][6] = x*x*y*y; + mat[index][7] = y*y; + mat[index][8] = x*y*y; + mat[index][9] = x*x*x; + mat[index][10] = x*x*x*y; + mat[index][11] = x*x*x*y*y; + mat[index][12] = x*x*x*y*y*y; + mat[index][13] = y*y*y; + mat[index][14] = x*y*y*y; + mat[index][15] = x*x*y*y*y; + } } - } - // we want to find the values of a,b,c to give the best - // fit - - Try { - surface_coefficients - = qr_method( &tz[0], - &t1[0], &t2[0], &t3[0], &t4[0], &t5[0], &t6[0], &t7[0], &t8[0], - &t9[0], &t10[0], &t11[0], &t12[0], &t13[0], &t14[0], &t15[0], - nobs, npred - ); - SG_LOG(SG_GENERAL, SG_DEBUG, "surface_coefficients size = " << surface_coefficients.nrows() ); - } - CatchAll { cout << BaseException::what(); } + // Do QR decompostion + JAMA::QR qr( mat ); + // find the least squares solution using the QR factors + surface_coefficients = qr.solve(zmat); } @@ -463,15 +368,13 @@ double TGAptSurface::query( double lon_deg, double lat_deg ) { double x = lon_deg - offset.x(); double y = lat_deg - offset.y(); - ColumnVector A = surface_coefficients; + TNT::Array1D A = surface_coefficients; - double result = A(1) + A(2)*x + A(3)*x*y + A(4)*y + A(5)*x*x + A(6)*x*x*y - + A(7)*x*x*y*y + A(8)*y*y + A(9)*x*y*y + A(10)*x*x*x + A(11)*x*x*x*y - + A(12)*x*x*x*y*y + A(13)*x*x*x*y*y*y + A(14)*y*y*y + A(15)*x*y*y*y - + A(16)*x*x*y*y*y; + double result = A[0] + A[1]*x + A[2]*x*y + A[3]*y + A[4]*x*x + A[5]*x*x*y + + A[6]*x*x*y*y + A[7]*y*y + A[8]*x*y*y + A[9]*x*x*x + A[10]*x*x*x*y + + A[11]*x*x*x*y*y + A[12]*x*x*x*y*y*y + A[13]*y*y*y + A[14]*x*y*y*y + + A[15]*x*x*y*y*y; result += offset.z(); - // printf("result = %.6f %.6f %.2f\n", lon_deg, lat_deg, result); - return result; } diff --git a/src/Airports/GenAirports850/apt_surface.hxx b/src/Airports/GenAirports850/apt_surface.hxx index 8c9c4224..ddd68bc3 100644 --- a/src/Airports/GenAirports850/apt_surface.hxx +++ b/src/Airports/GenAirports850/apt_surface.hxx @@ -19,8 +19,6 @@ // along with this program; if not, write to the Free Software // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. // -// $Id: apt_surface.hxx,v 1.9 2005-09-09 20:47:04 curt Exp $ -// #ifndef _APT_SURFACE_HXX @@ -28,15 +26,10 @@ #include +#include -// 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 // need matrix applications -#include // need matrix output routines - -#include +#include +#include /*** @@ -62,12 +55,10 @@ public: inline Point3D element( int col, int row ) { int index = ( row * _cols ) + col; if ( col < 0 || col >= _cols ) { - cout << "column out of bounds on read (" << col << " >= " << _cols << ")" - << endl; + SG_LOG(SG_GENERAL, SG_WARN, "column out of bounds on read (" << col << " >= " << _cols << ")"); int *p = 0; *p = 1; // force crash } else if ( row < 0 || row >= _rows ) { - cout << "row out of bounds on read (" << row << " >= " << _rows << ")" - << endl; + SG_LOG(SG_GENERAL, SG_WARN, "row out of bounds on read (" << row << " >= " << _rows << ")"); int *p = 0; *p = 1; // force crash } return m[index]; @@ -76,12 +67,10 @@ public: inline void set( int col, int row, Point3D p ) { int index = ( row * _cols ) + col; if ( col < 0 || col >= _cols ) { - cout << "column out of bounds on set (" << col << " >= " << _cols << ")" - << endl; + SG_LOG(SG_GENERAL, SG_WARN,"column out of bounds on set (" << col << " >= " << _cols << ")"); int *p = 0; *p = 1; // force crash } else if ( row < 0 || row >= _rows ) { - cout << "row out of bounds on set (" << row << " >= " << _rows << ")" - << endl; + SG_LOG(SG_GENERAL, SG_WARN,"row out of bounds on set (" << row << " >= " << _rows << ")"); int *p = 0; *p = 1; // force crash } m[index] = p; @@ -111,7 +100,7 @@ private: // The actual nurbs surface approximation for the airport SimpleMatrix *Pts; - ColumnVector surface_coefficients; + TNT::Array1D surface_coefficients; Point3D min_deg, max_deg; @@ -127,7 +116,7 @@ public: // Constructor, specify min and max coordinates of desired area in // lon/lat degrees, also please specify an "average" airport // elevations in meters. - TGAptSurface( const string &path, const string_list& elev_src, + TGAptSurface( const std::string &path, const string_list& elev_src, Point3D _min_deg, Point3D _max_deg, double _average_elev_m ); // Destructor diff --git a/src/Airports/GenAirports850/elevations.cxx b/src/Airports/GenAirports850/elevations.cxx index f507ef3f..d2595fca 100644 --- a/src/Airports/GenAirports850/elevations.cxx +++ b/src/Airports/GenAirports850/elevations.cxx @@ -19,20 +19,11 @@ // along with this program; if not, write to the Free Software // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. // -// $Id: elevations.cxx,v 1.8 2005-12-19 16:51:25 curt Exp $ -// #ifdef HAVE_CONFIG_H # include #endif -// 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 // need matrix applications -#include // need matrix output routines - #include #include #include @@ -47,7 +38,7 @@ // lookup node elevations for each point in the point_list. Returns // average of all points. Doesn't modify the original list. -double tgAverageElevation( const string &root, const string_list elev_src, +double tgAverageElevation( const std::string &root, const string_list elev_src, const point_list points_source ) { bool done = false; @@ -82,13 +73,13 @@ double tgAverageElevation( const string &root, const string_list elev_src, if ( found_one ) { SGBucket b( first.x(), first.y() ); - string base = b.gen_base_path(); + std::string base = b.gen_base_path(); // try the various elevation sources i = 0; bool found_file = false; while ( !found_file && i < elev_src.size() ) { - string array_path = root + "/" + elev_src[i] + "/" + base + "/" + b.gen_index_str(); + std::string array_path = root + "/" + elev_src[i] + "/" + base + "/" + b.gen_index_str(); if ( array.open(array_path) ) { found_file = true; @@ -145,7 +136,7 @@ double tgAverageElevation( const string &root, const string_list elev_src, // lookup node elevations for each point in the specified simple // matrix. Returns average of all points. -void tgCalcElevations( const string &root, const string_list elev_src, +void tgCalcElevations( const std::string &root, const string_list elev_src, SimpleMatrix &Pts, const double average ) { bool done = false; @@ -182,13 +173,13 @@ void tgCalcElevations( const string &root, const string_list elev_src, if ( found_one ) { SGBucket b( first.x(), first.y() ); - string base = b.gen_base_path(); + std::string base = b.gen_base_path(); // try the various elevation sources j = 0; bool found_file = false; while ( !found_file && j < (int)elev_src.size() ) { - string array_path = root + "/" + elev_src[j] + "/" + base + "/" + b.gen_index_str(); + std::string array_path = root + "/" + elev_src[j] + "/" + base + "/" + b.gen_index_str(); if ( array.open(array_path) ) { found_file = true; diff --git a/src/Airports/GenAirports850/elevations.hxx b/src/Airports/GenAirports850/elevations.hxx index deb5d416..f1f8f292 100644 --- a/src/Airports/GenAirports850/elevations.hxx +++ b/src/Airports/GenAirports850/elevations.hxx @@ -19,16 +19,6 @@ // along with this program; if not, write to the Free Software // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. // -// $Id: elevations.hxx,v 1.4 2005-09-09 15:05:15 curt Exp $ -// - - -// 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 // need matrix applications -#include // need matrix output routines #include #include @@ -43,12 +33,12 @@ // lookup node elevations for each point in the point_list. Returns // average of all points. Doesn't modify the original list. -double tgAverageElevation( const string &root, const string_list elev_src, +double tgAverageElevation( const std::string &root, const string_list elev_src, const point_list points_source ); // lookup node elevations for each point in the specified nurbs++ // matrix. -void tgCalcElevations( const string &root, const string_list elev_src, SimpleMatrix &Pts, double average ); +void tgCalcElevations( const std::string &root, const string_list elev_src, SimpleMatrix &Pts, double average ); // clamp all elevations to the specified range void tgClampElevations( SimpleMatrix &Pts, double center_m, double max_clamp_m ); diff --git a/src/Lib/Geometry/TNT/jama_qr.h b/src/Lib/Geometry/TNT/jama_qr.h new file mode 100644 index 00000000..98d37f53 --- /dev/null +++ b/src/Lib/Geometry/TNT/jama_qr.h @@ -0,0 +1,326 @@ +#ifndef JAMA_QR_H +#define JAMA_QR_H + +#include "tnt_array1d.h" +#include "tnt_array2d.h" +#include "tnt_math_utils.h" + +namespace JAMA +{ + +/** +

+ Classical QR Decompisition: + for an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n + orthogonal matrix Q and an n-by-n upper triangular matrix R so that + A = Q*R. +

+ The QR decompostion always exists, even if the matrix does not have + full rank, so the constructor will never fail. The primary use of the + QR decomposition is in the least squares solution of nonsquare systems + of simultaneous linear equations. This will fail if isFullRank() + returns 0 (false). + +

+ The Q and R factors can be retrived via the getQ() and getR() + methods. Furthermore, a solve() method is provided to find the + least squares solution of Ax=b using the QR factors. + +

+ (Adapted from JAMA, a Java Matrix Library, developed by jointly + by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama). +*/ + +template +class QR { + + + /** Array for internal storage of decomposition. + @serial internal array storage. + */ + + TNT::Array2D QR_; + + /** Row and column dimensions. + @serial column dimension. + @serial row dimension. + */ + int m, n; + + /** Array for internal storage of diagonal of R. + @serial diagonal of R. + */ + TNT::Array1D Rdiag; + + +public: + +/** + Create a QR factorization object for A. + + @param A rectangular (m>=n) matrix. +*/ + QR(const TNT::Array2D &A) /* constructor */ + { + QR_ = A.copy(); + m = A.dim1(); + n = A.dim2(); + Rdiag = TNT::Array1D(n); + int i=0, j=0, k=0; + + // Main loop. + for (k = 0; k < n; k++) { + // Compute 2-norm of k-th column without under/overflow. + Real nrm = 0; + for (i = k; i < m; i++) { + nrm = TNT::hypot(nrm,QR_[i][k]); + } + + if (nrm != 0.0) { + // Form k-th Householder vector. + if (QR_[k][k] < 0) { + nrm = -nrm; + } + for (i = k; i < m; i++) { + QR_[i][k] /= nrm; + } + QR_[k][k] += 1.0; + + // Apply transformation to remaining columns. + for (j = k+1; j < n; j++) { + Real s = 0.0; + for (i = k; i < m; i++) { + s += QR_[i][k]*QR_[i][j]; + } + s = -s/QR_[k][k]; + for (i = k; i < m; i++) { + QR_[i][j] += s*QR_[i][k]; + } + } + } + Rdiag[k] = -nrm; + } + } + + +/** + Flag to denote the matrix is of full rank. + + @return 1 if matrix is full rank, 0 otherwise. +*/ + int isFullRank() const + { + for (int j = 0; j < n; j++) + { + if (Rdiag[j] == 0) + return 0; + } + return 1; + } + + + + + /** + + Retreive the Householder vectors from QR factorization + @returns lower trapezoidal matrix whose columns define the reflections + */ + + TNT::Array2D getHouseholder (void) const + { + TNT::Array2D H(m,n); + + /* note: H is completely filled in by algorithm, so + initializaiton of H is not necessary. + */ + for (int i = 0; i < m; i++) + { + for (int j = 0; j < n; j++) + { + if (i >= j) { + H[i][j] = QR_[i][j]; + } else { + H[i][j] = 0.0; + } + } + } + return H; + } + + + + /** Return the upper triangular factor, R, of the QR factorization + @return R + */ + + TNT::Array2D getR() const + { + TNT::Array2D R(n,n); + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + if (i < j) { + R[i][j] = QR_[i][j]; + } else if (i == j) { + R[i][j] = Rdiag[i]; + } else { + R[i][j] = 0.0; + } + } + } + return R; + } + + + + + + /** + Generate and return the (economy-sized) orthogonal factor + @param Q the (ecnomy-sized) orthogonal factor (Q*R=A). + */ + + TNT::Array2D getQ() const + { + int i=0, j=0, k=0; + + TNT::Array2D Q(m,n); + for (k = n-1; k >= 0; k--) { + for (i = 0; i < m; i++) { + Q[i][k] = 0.0; + } + Q[k][k] = 1.0; + for (j = k; j < n; j++) { + if (QR_[k][k] != 0) { + Real s = 0.0; + for (i = k; i < m; i++) { + s += QR_[i][k]*Q[i][j]; + } + s = -s/QR_[k][k]; + for (i = k; i < m; i++) { + Q[i][j] += s*QR_[i][k]; + } + } + } + } + return Q; + } + + + /** Least squares solution of A*x = b + @param B m-length array (vector). + @return x n-length array (vector) that minimizes the two norm of Q*R*X-B. + If B is non-conformant, or if QR.isFullRank() is false, + the routine returns a null (0-length) vector. + */ + + TNT::Array1D solve(const TNT::Array1D &b) const + { + if (b.dim1() != m) /* arrays must be conformant */ + return TNT::Array1D(); + + if ( !isFullRank() ) /* matrix is rank deficient */ + { + return TNT::Array1D(); + } + + TNT::Array1D x = b.copy(); + + // Compute Y = transpose(Q)*b + for (int k = 0; k < n; k++) + { + Real s = 0.0; + for (int i = k; i < m; i++) + { + s += QR_[i][k]*x[i]; + } + s = -s/QR_[k][k]; + for (int i = k; i < m; i++) + { + x[i] += s*QR_[i][k]; + } + } + // Solve R*X = Y; + for (int k = n-1; k >= 0; k--) + { + x[k] /= Rdiag[k]; + for (int i = 0; i < k; i++) { + x[i] -= x[k]*QR_[i][k]; + } + } + + + /* return n x nx portion of X */ + TNT::Array1D x_(n); + for (int i=0; i solve(const TNT::Array2D &B) const + { + if (B.dim1() != m) /* arrays must be conformant */ + return TNT::Array2D(0,0); + + if ( !isFullRank() ) /* matrix is rank deficient */ + { + return TNT::Array2D(0,0); + } + + int nx = B.dim2(); + TNT::Array2D X = B.copy(); + int i=0, j=0, k=0; + + // Compute Y = transpose(Q)*B + for (k = 0; k < n; k++) { + for (j = 0; j < nx; j++) { + Real s = 0.0; + for (i = k; i < m; i++) { + s += QR_[i][k]*X[i][j]; + } + s = -s/QR_[k][k]; + for (i = k; i < m; i++) { + X[i][j] += s*QR_[i][k]; + } + } + } + // Solve R*X = Y; + for (k = n-1; k >= 0; k--) { + for (j = 0; j < nx; j++) { + X[k][j] /= Rdiag[k]; + } + for (i = 0; i < k; i++) { + for (j = 0; j < nx; j++) { + X[i][j] -= X[k][j]*QR_[i][k]; + } + } + } + + + /* return n x nx portion of X */ + TNT::Array2D X_(n,nx); + for (i=0; i +#include + +#ifdef TNT_BOUNDS_CHECK +#include +#endif + + +#include "tnt_i_refvec.h" + +namespace TNT +{ + +template +class Array1D +{ + + private: + + /* ... */ + i_refvec v_; + int n_; + T* data_; /* this normally points to v_.begin(), but + * could also point to a portion (subvector) + * of v_. + */ + + void copy_(T* p, const T* q, int len) const; + void set_(T* begin, T* end, const T& val); + + + public: + + typedef T value_type; + + + Array1D(); + explicit Array1D(int n); + Array1D(int n, const T &a); + Array1D(int n, T *a); + inline Array1D(const Array1D &A); + inline operator T*(); + inline operator const T*(); + inline Array1D & operator=(const T &a); + inline Array1D & operator=(const Array1D &A); + inline Array1D & ref(const Array1D &A); + Array1D copy() const; + Array1D & inject(const Array1D & A); + inline T& operator[](int i); + inline T& operator()(int i); + inline const T& operator[](int i) const; + inline const T& operator()(int i) const; + inline int dim1() const; + inline int dim() const; + ~Array1D(); + + + /* ... extended interface ... */ + + inline int ref_count() const; + inline Array1D subarray(int i0, int i1); + +}; + + + + +template +Array1D::Array1D() : v_(), n_(0), data_(0) {} + +template +Array1D::Array1D(const Array1D &A) : v_(A.v_), n_(A.n_), + data_(A.data_) +{ +#ifdef TNT_DEBUG + std::cout << "Created Array1D(const Array1D &A) \n"; +#endif + +} + +template +Array1D::Array1D(int n) : v_(n), n_(n), data_(v_.begin()) +{ +#ifdef TNT_DEBUG + std::cout << "Created Array1D(int n) \n"; +#endif +} + +template +Array1D::Array1D(int n, const T &val) : v_(n), n_(n), data_(v_.begin()) +{ +#ifdef TNT_DEBUG + std::cout << "Created Array1D(int n, const T& val) \n"; +#endif + set_(data_, data_+ n, val); + +} + +template +Array1D::Array1D(int n, T *a) : v_(a), n_(n) , data_(v_.begin()) +{ +#ifdef TNT_DEBUG + std::cout << "Created Array1D(int n, T* a) \n"; +#endif +} + +template +inline Array1D::operator T*() +{ + return &(v_[0]); +} + + +template +inline Array1D::operator const T*() +{ + return &(v_[0]); +} + + + +template +inline T& Array1D::operator[](int i) +{ +#ifdef TNT_BOUNDS_CHECK + assert(i>= 0); + assert(i < n_); +#endif + return data_[i]; +} + +template +inline const T& Array1D::operator[](int i) const +{ +#ifdef TNT_BOUNDS_CHECK + assert(i>= 0); + assert(i < n_); +#endif + return data_[i]; +} + + + +template +inline T& Array1D::operator()(int i) +{ +#ifdef TNT_BOUNDS_CHECK + assert(i>= 0); + assert(i < n_); +#endif + return data_[i-1]; +} + +template +inline const T& Array1D::operator()(int i) const +{ +#ifdef TNT_BOUNDS_CHECK + assert(i>= 0); + assert(i < n_); +#endif + return data_[i-1]; +} + + + + +template +Array1D & Array1D::operator=(const T &a) +{ + set_(data_, data_+n_, a); + return *this; +} + +template +Array1D Array1D::copy() const +{ + Array1D A( n_); + copy_(A.data_, data_, n_); + + return A; +} + + +template +Array1D & Array1D::inject(const Array1D &A) +{ + if (A.n_ == n_) + copy_(data_, A.data_, n_); + + return *this; +} + + + + + +template +Array1D & Array1D::ref(const Array1D &A) +{ + if (this != &A) + { + v_ = A.v_; /* operator= handles the reference counting. */ + n_ = A.n_; + data_ = A.data_; + + } + return *this; +} + +template +Array1D & Array1D::operator=(const Array1D &A) +{ + return ref(A); +} + +template +inline int Array1D::dim1() const { return n_; } + +template +inline int Array1D::dim() const { return n_; } + +template +Array1D::~Array1D() {} + + +/* ............................ exented interface ......................*/ + +template +inline int Array1D::ref_count() const +{ + return v_.ref_count(); +} + +template +inline Array1D Array1D::subarray(int i0, int i1) +{ + if ((i0 >= 0) && (i1 < n_) || (i0 <= i1)) + { + Array1D X(*this); /* create a new instance of this array. */ + X.n_ = i1-i0+1; + X.data_ += i0; + + return X; + } + else + { + return Array1D(); + } +} + + +/* private internal functions */ + + +template +void Array1D::set_(T* begin, T* end, const T& a) +{ + for (T* p=begin; p +void Array1D::copy_(T* p, const T* q, int len) const +{ + T *end = p + len; + while (p +#include +#ifdef TNT_BOUNDS_CHECK +#include +#endif + +#include "tnt_array1d.h" + +namespace TNT +{ + + +/** + Ttwo-dimensional numerical array which + looks like a conventional C multiarray. + Storage corresponds to C (row-major) ordering. + Elements are accessed via A[i][j] notation for 0-based indexing, + and A(i,j) for 1-based indexing.. + +

+ Array assignment is by reference (i.e. shallow assignment). + That is, B=A implies that the A and B point to the + same array, so modifications to the elements of A + will be reflected in B. If an independent copy + is required, then B = A.copy() can be used. Note + that this facilitates returning arrays from functions + without relying on compiler optimizations to eliminate + extensive data copying. + +

+ The indexing and layout of this array object makes + it compatible with C and C++ algorithms that utilize + the familiar C[i][j] notation. This includes numerous + + textbooks, such as Numercial Recipes, and various + public domain codes. + +*/ +template +class Array2D +{ + + + private: + + + + Array1D data_; + Array1D v_; + int m_; + int n_; + + public: + + + +/** + Used to determined the data type of array entries. This is most + commonly used when requiring scalar temporaries in templated algorithms + that have TNT arrays as input. For example, + 

+  template < class ArrayTwoD >
+  void foo(ArrayTwoD &A)
+  {
+    A::value_type first_entry = A[0][0];
+    ...
+  }
+
+*/ + typedef T value_type; + +/** + Create a null array. This is not the same + as Array2D(0,0), which consumes some memory overhead. +*/ + + Array2D(); + + +/** + Create a new (m x n) array, without initalizing elements. (This + encurs an O(1) operation cost, rather than a O(m*n) cost.) + + @param m the first (row) dimension of the new matrix. + @param n the second (column) dimension of the new matrix. +*/ + Array2D(int m, int n); + + +/** + Create a new (m x n) array, as a view of an existing one-dimensional + array stored in row-major order, i.e. right-most dimension varying fastest. + Note that the storage for this pre-existing array will + never be destroyed by TNT. + + @param m the first (row) dimension of the new matrix. + @param n the second (column) dimension of the new matrix. + @param a the one dimensional C array to use as data storage for + the array. +*/ + Array2D(int m, int n, T *a); + + + +/** + Create a new (m x n) array, initializing array elements to + constant specified by argument. Most often used to + create an array of zeros, as in A(m, n, 0.0). + + @param m the first (row) dimension of the new matrix. + @param n the second (column) dimension of the new matrix. + @param val the constant value to set all elements of the new array to. +*/ + + Array2D(int m, int n, const T &val); + + +/** + Copy constructor. Array data is NOT copied, but shared. + Thus, in Array2D B(A), subsequent changes to A will + be reflected in B. For an indepent copy of A, use + Array2D B(A.copy()), or B = A.copy(), instead. +*/ + inline Array2D(const Array2D &A); + + +/** + Convert 2D array into a regular multidimensional C pointer. Most often + called automatically when calling C interfaces that expect things like + double** rather than Array2D. + +*/ + inline operator T**(); + +/** + Convert a const 2D array into a const multidimensional C pointer. + Most often called automatically when calling C interfaces that expect + things like "const double**" rather than "const Array2D&". + +*/ + + inline operator const T**() const; + + +/** + Assign all elements of array the same value. + + @param val the value to assign each element. +*/ + inline Array2D & operator=(const T &val); + + +/** + Assign one Array2D to another. (This is a shallow-assignement operation, + and it is the identical semantics to ref(A). + + @param A the array to assign this one to. +*/ + inline Array2D & operator=(const Array2D &A); + + + inline Array2D & ref(const Array2D &A); + Array2D copy() const; + Array2D & inject(const Array2D & A); + inline T* operator[](int i); + inline const T* operator[](int i) const; + inline int dim1() const; + inline int dim2() const; + ~Array2D(); + + /* extended interface (not part of the standard) */ + + + inline int ref_count(); + inline int ref_count_data(); + inline int ref_count_dim1(); + Array2D subarray(int i0, int i1, int j0, int j1); + +}; + + +/** + Create a new (m x n) array, WIHOUT initializing array elements. + To create an initialized array of constants, see Array2D(m,n,value). + +

+ This version avoids the O(m*n) initialization overhead and + is used just before manual assignment. + + @param m the first (row) dimension of the new matrix. + @param n the second (column) dimension of the new matrix. +*/ +template +Array2D::Array2D() : data_(), v_(), m_(0), n_(0) {} + + +template +Array2D::Array2D(const Array2D &A) : data_(A.data_), v_(A.v_), + m_(A.m_), n_(A.n_) {} + + + + +template +Array2D::Array2D(int m, int n) : data_(m*n), v_(m), m_(m), n_(n) +{ + if (m>0 && n>0) + { + T* p = &(data_[0]); + for (int i=0; i +Array2D::Array2D(int m, int n, const T &val) : data_(m*n), v_(m), + m_(m), n_(n) +{ + if (m>0 && n>0) + { + data_ = val; + T* p = &(data_[0]); + for (int i=0; i +Array2D::Array2D(int m, int n, T *a) : data_(m*n, a), v_(m), m_(m), n_(n) +{ + if (m>0 && n>0) + { + T* p = &(data_[0]); + + for (int i=0; i +inline T* Array2D::operator[](int i) +{ +#ifdef TNT_BOUNDS_CHECK + assert(i >= 0); + assert(i < m_); +#endif + +return v_[i]; + +} + + +template +inline const T* Array2D::operator[](int i) const +{ +#ifdef TNT_BOUNDS_CHECK + assert(i >= 0); + assert(i < m_); +#endif + +return v_[i]; + +} + +template +Array2D & Array2D::operator=(const T &a) +{ + /* non-optimzied, but will work with subarrays in future verions */ + + for (int i=0; i +Array2D Array2D::copy() const +{ + Array2D A(m_, n_); + + for (int i=0; i +Array2D & Array2D::inject(const Array2D &A) +{ + if (A.m_ == m_ && A.n_ == n_) + { + for (int i=0; i +Array2D & Array2D::ref(const Array2D &A) +{ + if (this != &A) + { + v_ = A.v_; + data_ = A.data_; + m_ = A.m_; + n_ = A.n_; + + } + return *this; +} + + + +template +Array2D & Array2D::operator=(const Array2D &A) +{ + return ref(A); +} + +template +inline int Array2D::dim1() const { return m_; } + +template +inline int Array2D::dim2() const { return n_; } + + +template +Array2D::~Array2D() {} + + + + +template +inline Array2D::operator T**() +{ + return &(v_[0]); +} +template +inline Array2D::operator const T**() const +{ + return static_cast(&(v_[0])); +} + +/* ............... extended interface ............... */ +/** + Create a new view to a subarray defined by the boundaries + [i0][i0] and [i1][j1]. The size of the subarray is + (i1-i0) by (j1-j0). If either of these lengths are zero + or negative, the subarray view is null. + +*/ +template +Array2D Array2D::subarray(int i0, int i1, int j0, int j1) +{ + Array2D A; + int m = i1-i0+1; + int n = j1-j0+1; + + /* if either length is zero or negative, this is an invalide + subarray. return a null view. + */ + if (m<1 || n<1) + return A; + + A.data_ = data_; + A.m_ = m; + A.n_ = n; + A.v_ = Array1D(m); + T* p = &(data_[0]) + i0 * n_ + j0; + for (int i=0; i +inline int Array2D::ref_count() +{ + return ref_count_data(); +} + + + +template +inline int Array2D::ref_count_data() +{ + return data_.ref_count(); +} + +template +inline int Array2D::ref_count_dim1() +{ + return v_.ref_count(); +} + + + + +} /* namespace TNT */ + +#endif +/* TNT_ARRAY2D_H */ + diff --git a/src/Lib/Geometry/TNT/tnt_i_refvec.h b/src/Lib/Geometry/TNT/tnt_i_refvec.h new file mode 100644 index 00000000..5a67eb57 --- /dev/null +++ b/src/Lib/Geometry/TNT/tnt_i_refvec.h @@ -0,0 +1,243 @@ +/* +* +* Template Numerical Toolkit (TNT) +* +* Mathematical and Computational Sciences Division +* National Institute of Technology, +* Gaithersburg, MD USA +* +* +* This software was developed at the National Institute of Standards and +* Technology (NIST) by employees of the Federal Government in the course +* of their official duties. Pursuant to title 17 Section 105 of the +* United States Code, this software is not subject to copyright protection +* and is in the public domain. NIST assumes no responsibility whatsoever for +* its use by other parties, and makes no guarantees, expressed or implied, +* about its quality, reliability, or any other characteristic. +* +*/ + + + +#ifndef TNT_I_REFVEC_H +#define TNT_I_REFVEC_H + +#include +#include + +#ifdef TNT_BOUNDS_CHECK +#include +#endif + +#ifndef NULL +#define NULL 0 +#endif + +namespace TNT +{ +/* + Internal representation of ref-counted array. The TNT + arrays all use this building block. + +

+ If an array block is created by TNT, then every time + an assignment is made, the left-hand-side reference + is decreased by one, and the right-hand-side refernce + count is increased by one. If the array block was + external to TNT, the refernce count is a NULL pointer + regardless of how many references are made, since the + memory is not freed by TNT. + + + +*/ +template +class i_refvec +{ + + + private: + T* data_; + int *ref_count_; + + + public: + + i_refvec(); + explicit i_refvec(int n); + inline i_refvec(T* data); + inline i_refvec(const i_refvec &v); + inline T* begin(); + inline const T* begin() const; + inline T& operator[](int i); + inline const T& operator[](int i) const; + inline i_refvec & operator=(const i_refvec &V); + void copy_(T* p, const T* q, const T* e); + void set_(T* p, const T* b, const T* e); + inline int ref_count() const; + inline int is_null() const; + inline void destroy(); + ~i_refvec(); + +}; + +template +void i_refvec::copy_(T* p, const T* q, const T* e) +{ + for (T* t=p; q +i_refvec::i_refvec() : data_(NULL), ref_count_(NULL) {} + +/** + In case n is 0 or negative, it does NOT call new. +*/ +template +i_refvec::i_refvec(int n) : data_(NULL), ref_count_(NULL) +{ + if (n >= 1) + { +#ifdef TNT_DEBUG + std::cout << "new data storage.\n"; +#endif + data_ = new T[n]; + ref_count_ = new int; + *ref_count_ = 1; + } +} + +template +inline i_refvec::i_refvec(const i_refvec &V): data_(V.data_), + ref_count_(V.ref_count_) +{ + if (V.ref_count_ != NULL) + (*(V.ref_count_))++; +} + + +template +i_refvec::i_refvec(T* data) : data_(data), ref_count_(NULL) {} + +template +inline T* i_refvec::begin() +{ + return data_; +} + +template +inline const T& i_refvec::operator[](int i) const +{ + return data_[i]; +} + +template +inline T& i_refvec::operator[](int i) +{ + return data_[i]; +} + + +template +inline const T* i_refvec::begin() const +{ + return data_; +} + + + +template +i_refvec & i_refvec::operator=(const i_refvec &V) +{ + if (this == &V) + return *this; + + + if (ref_count_ != NULL) + { + (*ref_count_) --; + if ((*ref_count_) == 0) + destroy(); + } + + data_ = V.data_; + ref_count_ = V.ref_count_; + + if (V.ref_count_ != NULL) + (*(V.ref_count_))++; + + return *this; +} + +template +void i_refvec::destroy() +{ + if (ref_count_ != NULL) + { +#ifdef TNT_DEBUG + std::cout << "destorying data... \n"; +#endif + delete ref_count_; + +#ifdef TNT_DEBUG + std::cout << "deleted ref_count_ ...\n"; +#endif + if (data_ != NULL) + delete []data_; +#ifdef TNT_DEBUG + std::cout << "deleted data_[] ...\n"; +#endif + data_ = NULL; + } +} + +/* +* return 1 is vector is empty, 0 otherwise +* +* if is_null() is false and ref_count() is 0, then +* +*/ +template +int i_refvec::is_null() const +{ + return (data_ == NULL ? 1 : 0); +} + +/* +* returns -1 if data is external, +* returns 0 if a is NULL array, +* otherwise returns the positive number of vectors sharing +* this data space. +*/ +template +int i_refvec::ref_count() const +{ + if (data_ == NULL) + return 0; + else + return (ref_count_ != NULL ? *ref_count_ : -1) ; +} + +template +i_refvec::~i_refvec() +{ + if (ref_count_ != NULL) + { + (*ref_count_)--; + + if (*ref_count_ == 0) + destroy(); + } +} + + +} /* namespace TNT */ + + + + + +#endif +/* TNT_I_REFVEC_H */ + diff --git a/src/Lib/Geometry/TNT/tnt_math_utils.h b/src/Lib/Geometry/TNT/tnt_math_utils.h new file mode 100644 index 00000000..4d81950f --- /dev/null +++ b/src/Lib/Geometry/TNT/tnt_math_utils.h @@ -0,0 +1,37 @@ +#ifndef MATH_UTILS_H +#define MATH_UTILS_H + + +/* needed for abs(), sqrt() below */ +#include + + + +namespace TNT +{ + +using namespace std; +/** + @returns hypotenuse of real (non-complex) scalars a and b by + avoiding underflow/overflow + using (a * sqrt( 1 + (b/a) * (b/a))), rather than + sqrt(a*a + b*b). +*/ +template +Real hypot(const Real &a, const Real &b) +{ + + if (a== 0) + return abs(b); + else + { + Real c = b/a; + return abs(a) * sqrt(1 + c*c); + } +} +} /* TNT namespace */ + + + +#endif +/* MATH_UTILS_H */