library reorg step 1
- no more Geometry - delete unused files
This commit is contained in:
Peter Sadrozinski
parent
4b97381559
commit
ce1f89f832
60 changed files with 209 additions and 6338 deletions
|
|
@ -26,7 +26,7 @@ add_executable(genapts850
|
|||
)
|
||||
|
||||
target_link_libraries(genapts850
|
||||
Polygon Geometry
|
||||
Polygon
|
||||
Array
|
||||
${GDAL_LIBRARY}
|
||||
${SIMGEAR_CORE_LIBRARIES}
|
||||
|
|
|
|||
|
|
@ -288,8 +288,6 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
|
|||
// runway lights
|
||||
tglightcontour_list rwy_lights;
|
||||
|
||||
//Point3D p;
|
||||
|
||||
bool make_shapefiles = false;
|
||||
|
||||
// parse main airport information
|
||||
|
|
@ -1064,114 +1062,6 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
|
|||
TGAptSurface apt_surf( root, elev_src, aptBounds, average );
|
||||
GENAPT_LOG(SG_GENERAL, SG_DEBUG, "Airport surface created");
|
||||
|
||||
// add base skirt (to hide potential cracks)
|
||||
// this has to happen after we've calculated the node elevations
|
||||
// but before we convert to wgs84 coordinates
|
||||
|
||||
#if 0 // do we still need a skirt?
|
||||
int uindex, lindex;
|
||||
for ( unsigned int i = 0; i < divided_base.Contours(); ++i )
|
||||
{
|
||||
strip_v.clear();
|
||||
strip_n.clear();
|
||||
strip_tc.clear();
|
||||
|
||||
// prime the pump ...
|
||||
uindex = nodes.find( divided_base.GetNode( i, 0 ) );
|
||||
if ( uindex >= 0 )
|
||||
{
|
||||
Point3D lower = geod_nodes[uindex] - Point3D(0, 0, 20);
|
||||
GENAPT_LOG(SG_GENERAL, SG_DEBUG, geod_nodes[uindex] << " <-> " << lower);
|
||||
lindex = nodes.simple_add( lower );
|
||||
geod_nodes.push_back( lower );
|
||||
strip_v.push_back( uindex );
|
||||
strip_v.push_back( lindex );
|
||||
|
||||
// use 'the' normal. We are pushing on two nodes so we
|
||||
// need to push on two normals.
|
||||
index = normals.unique_add( vn );
|
||||
strip_n.push_back( index );
|
||||
strip_n.push_back( index );
|
||||
}
|
||||
else
|
||||
{
|
||||
string message = "Ooops missing node when building skirt (in init)";
|
||||
GENAPT_LOG( SG_GENERAL, SG_ALERT, message << " " << p );
|
||||
throw sg_exception( message );
|
||||
}
|
||||
|
||||
// loop through the list
|
||||
for ( unsigned int j = 1; j < divided_base.ContourSize(i); ++j )
|
||||
{
|
||||
p = Point3D::fromSGGeod( divided_base.GetNode( i, j ) );
|
||||
uindex = nodes.find( p );
|
||||
if ( uindex >= 0 )
|
||||
{
|
||||
Point3D lower = geod_nodes[uindex] - Point3D(0, 0, 20);
|
||||
GENAPT_LOG(SG_GENERAL, SG_DEBUG, geod_nodes[uindex] << " <-> " << lower);
|
||||
lindex = nodes.simple_add( lower );
|
||||
geod_nodes.push_back( lower );
|
||||
strip_v.push_back( uindex );
|
||||
strip_v.push_back( lindex );
|
||||
|
||||
index = normals.unique_add( vn );
|
||||
strip_n.push_back( index );
|
||||
strip_n.push_back( index );
|
||||
}
|
||||
else
|
||||
{
|
||||
string message = "Ooops missing node when building skirt (in loop)";
|
||||
GENAPT_LOG( SG_GENERAL, SG_ALERT, message << " " << p );
|
||||
throw sg_exception( message );
|
||||
}
|
||||
}
|
||||
|
||||
// close off the loop
|
||||
p = Point3D::fromSGGeod( divided_base.GetNode( i, 0 ) );
|
||||
uindex = nodes.find( p );
|
||||
if ( uindex >= 0 )
|
||||
{
|
||||
Point3D lower = geod_nodes[uindex] - Point3D(0, 0, 20);
|
||||
GENAPT_LOG(SG_GENERAL, SG_DEBUG, geod_nodes[uindex] << " <-> " << lower);
|
||||
lindex = nodes.simple_add( lower );
|
||||
geod_nodes.push_back( lower );
|
||||
strip_v.push_back( uindex );
|
||||
strip_v.push_back( lindex );
|
||||
|
||||
index = normals.unique_add( vn );
|
||||
strip_n.push_back( index );
|
||||
strip_n.push_back( index );
|
||||
}
|
||||
else
|
||||
{
|
||||
string message = "Ooops missing node when building skirt (at end)";
|
||||
GENAPT_LOG( SG_GENERAL, SG_ALERT, message << " " << p );
|
||||
throw sg_exception( message );
|
||||
}
|
||||
|
||||
strips_v.push_back( strip_v );
|
||||
strips_n.push_back( strip_n );
|
||||
strip_materials.push_back( "Grass" );
|
||||
|
||||
std::vector < SGGeod > geodNodes;
|
||||
for ( unsigned int j = 0; j < nodes.get_node_list().size(); j++ )
|
||||
{
|
||||
Point3D node = nodes.get_node_list()[j];
|
||||
geodNodes.push_back( SGGeod::fromDegM( node.x(), node.y(), node.z() ) );
|
||||
}
|
||||
base_txs.clear();
|
||||
base_txs = sgCalcTexCoords( b, geodNodes, strip_v );
|
||||
|
||||
base_tc.clear();
|
||||
for ( unsigned int j = 0; j < base_txs.size(); ++j )
|
||||
{
|
||||
index = texcoords.simple_add( base_txs[j] );
|
||||
base_tc.push_back( index );
|
||||
}
|
||||
strips_tc.push_back( base_tc );
|
||||
}
|
||||
#endif
|
||||
|
||||
// add light points
|
||||
// pass one, calculate raw elevations from Array
|
||||
for ( unsigned int i = 0; i < rwy_lights.size(); ++i ) {
|
||||
|
|
@ -1179,28 +1069,6 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src )
|
|||
double light_elevation = calc_elevation( apt_surf, rwy_lights[i].GetNode(j), 0.0 );
|
||||
rwy_lights[i].SetElevation(j, light_elevation);
|
||||
}
|
||||
|
||||
// TODO : It doesn't look like this does anything... got back in history and make sure...
|
||||
//string flag = rwy_lights[i].GetFlag();
|
||||
//if ( flag != (string)"" )
|
||||
//{
|
||||
// GENAPT_LOG(SG_GENERAL, SG_INFO, " flag " << flag);
|
||||
// double max = -9999;
|
||||
// const_elev_map_iterator it = elevation_map.find( flag );
|
||||
// if ( it != elevation_map.end() )
|
||||
// {
|
||||
// max = elevation_map[flag];
|
||||
// }
|
||||
// for ( unsigned int j = 0; j < geod_light_nodes.size(); ++j )
|
||||
// {
|
||||
// if ( geod_light_nodes[j].z() > max )
|
||||
// {
|
||||
// max = geod_light_nodes[j].z();
|
||||
// }
|
||||
// }
|
||||
// elevation_map[flag] = max;
|
||||
// GENAPT_LOG( SG_GENERAL, SG_INFO, " " << flag << " max = " << max );
|
||||
//}
|
||||
}
|
||||
|
||||
GENAPT_LOG(SG_GENERAL, SG_INFO, "Done with lighting calc_elevations() num light polys is " << rwy_lights.size() );
|
||||
|
|
|
|||
|
|
@ -24,7 +24,7 @@
|
|||
# include <config.h>
|
||||
#endif
|
||||
|
||||
#include <Lib/Geometry/TNT/jama_qr.h>
|
||||
#include <Lib/Polygon/TNT/jama_qr.h>
|
||||
|
||||
#include <simgear/compiler.h>
|
||||
#include <simgear/constants.h>
|
||||
|
|
|
|||
|
|
@ -26,7 +26,7 @@
|
|||
|
||||
|
||||
#include <string>
|
||||
#include <Lib/Geometry/TNT/tnt_array2d.h>
|
||||
#include <Lib/Polygon/TNT/tnt_array2d.h>
|
||||
|
||||
#include <simgear/debug/logstream.hxx>
|
||||
#include <Lib/Polygon/polygon.hxx>
|
||||
|
|
|
|||
|
|
@ -3,7 +3,6 @@
|
|||
#include <simgear/debug/logstream.hxx>
|
||||
|
||||
#include <Polygon/polygon.hxx>
|
||||
#include <Geometry/poly_support.hxx>
|
||||
|
||||
#include "global.hxx"
|
||||
#include "beznode.hxx"
|
||||
|
|
|
|||
|
|
@ -17,8 +17,6 @@
|
|||
#include <simgear/debug/logstream.hxx>
|
||||
#include <simgear/math/sg_geodesy.hxx>
|
||||
|
||||
#include <Geometry/poly_support.hxx>
|
||||
|
||||
#include "global.hxx"
|
||||
#include "apt_math.hxx"
|
||||
#include "helipad.hxx"
|
||||
|
|
|
|||
|
|
@ -1,7 +1,5 @@
|
|||
#include <stdlib.h>
|
||||
|
||||
#include <Geometry/poly_support.hxx>
|
||||
|
||||
#include "global.hxx"
|
||||
#include "beznode.hxx"
|
||||
#include "linearfeature.hxx"
|
||||
|
|
|
|||
|
|
@ -22,7 +22,6 @@
|
|||
#include <simgear/misc/sg_path.hxx>
|
||||
#include <simgear/misc/sgstream.hxx>
|
||||
|
||||
#include <Geometry/poly_support.hxx>
|
||||
#include <Include/version.h>
|
||||
|
||||
#include "scheduler.hxx"
|
||||
|
|
@ -306,7 +305,7 @@ int main(int argc, char **argv)
|
|||
}
|
||||
|
||||
// Initialize shapefile support (for debugging)
|
||||
tgShapefileInit();
|
||||
tgShapefile::Init();
|
||||
|
||||
sg_gzifstream in( input_file );
|
||||
if ( !in.is_open() )
|
||||
|
|
|
|||
|
|
@ -5,7 +5,6 @@
|
|||
#include <simgear/bucket/newbucket.hxx>
|
||||
#include <simgear/math/SGGeodesy.hxx>
|
||||
|
||||
#include <Geometry/poly_support.hxx>
|
||||
#include <Polygon/polygon.hxx>
|
||||
|
||||
#include "global.hxx"
|
||||
|
|
|
|||
|
|
@ -26,8 +26,6 @@
|
|||
|
||||
#include <boost/lexical_cast.hpp>
|
||||
|
||||
#include <Geometry/poly_support.hxx>
|
||||
|
||||
#include "global.hxx"
|
||||
#include "beznode.hxx"
|
||||
#include "runway.hxx"
|
||||
|
|
|
|||
|
|
@ -10,7 +10,7 @@
|
|||
#include <simgear/timing/timestamp.hxx>
|
||||
#include <simgear/threads/SGThread.hxx>
|
||||
#include <simgear/threads/SGQueue.hxx>
|
||||
#include <Geometry/rectangle.hxx>
|
||||
#include <Polygon/rectangle.hxx>
|
||||
#include "airport.hxx"
|
||||
|
||||
#define P_STATE_INIT (0)
|
||||
|
|
|
|||
|
|
@ -27,7 +27,7 @@ set_target_properties(tg-construct PROPERTIES
|
|||
"DEFAULT_USGS_MAPFILE=\"${PKGDATADIR}/usgsmap.txt\";DEFAULT_PRIORITIES_FILE=\"${PKGDATADIR}/default_priorities.txt\"" )
|
||||
|
||||
target_link_libraries(tg-construct
|
||||
Polygon Geometry
|
||||
Polygon
|
||||
Array landcover
|
||||
${Boost_LIBRARIES}
|
||||
${GDAL_LIBRARY}
|
||||
|
|
@ -43,7 +43,7 @@ add_executable(cliptst
|
|||
cliptst.cxx)
|
||||
|
||||
target_link_libraries(cliptst
|
||||
Polygon Geometry
|
||||
Polygon
|
||||
${GDAL_LIBRARY}
|
||||
${SIMGEAR_CORE_LIBRARIES}
|
||||
${SIMGEAR_CORE_LIBRARY_DEPENDENCIES}
|
||||
|
|
|
|||
|
|
@ -14,10 +14,6 @@
|
|||
|
||||
#include <Polygon/clipper.hpp>
|
||||
#include <Polygon/polygon.hxx>
|
||||
#include <Geometry/poly_support.hxx>
|
||||
|
||||
//#include "windows.h"
|
||||
//#include <conio.h>
|
||||
|
||||
//---------------------------------------------------------------------------
|
||||
|
||||
|
|
@ -263,7 +259,7 @@ int main(int argc, char* argv[])
|
|||
}
|
||||
|
||||
if (show_svg) {
|
||||
tgShapefileInit();
|
||||
tgShapefile::Init();
|
||||
clipper_to_shapefile( subject, "./clptst_subject" );
|
||||
clipper_to_shapefile( clip, "./clptst_clip" );
|
||||
}
|
||||
|
|
|
|||
|
|
@ -24,7 +24,6 @@
|
|||
|
||||
#include <simgear/debug/logstream.hxx>
|
||||
#include <Include/version.h>
|
||||
#include <Geometry/poly_support.hxx>
|
||||
|
||||
#include "tgconstruct.hxx"
|
||||
#include "usgs.hxx"
|
||||
|
|
@ -73,7 +72,7 @@ int main(int argc, char **argv) {
|
|||
sglog().setLogLevels( SG_ALL, SG_INFO );
|
||||
|
||||
// Initialize shapefile support (for debugging)
|
||||
tgShapefileInit();
|
||||
tgShapefile::Init();
|
||||
|
||||
//
|
||||
// Parse the command-line arguments.
|
||||
|
|
|
|||
|
|
@ -34,7 +34,7 @@
|
|||
#define TG_MAX_AREA_TYPES 128
|
||||
|
||||
#include <Array/array.hxx>
|
||||
#include <Geometry/tg_nodes.hxx>
|
||||
#include <Polygon/tg_nodes.hxx>
|
||||
#include <landcover/landcover.hxx>
|
||||
|
||||
#include "tglandclass.hxx"
|
||||
|
|
|
|||
|
|
@ -30,9 +30,6 @@
|
|||
//#include <simgear/compiler.h>
|
||||
#include <simgear/debug/logstream.hxx>
|
||||
|
||||
#include <Geometry/poly_support.hxx>
|
||||
#include <Geometry/poly_extra.hxx>
|
||||
|
||||
#include "tgconstruct.hxx"
|
||||
|
||||
using std::string;
|
||||
|
|
|
|||
|
|
@ -27,7 +27,6 @@
|
|||
#include <sstream>
|
||||
|
||||
#include <simgear/debug/logstream.hxx>
|
||||
#include <Geometry/poly_support.hxx>
|
||||
|
||||
#include "tgconstruct.hxx"
|
||||
|
||||
|
|
|
|||
|
|
@ -25,7 +25,6 @@
|
|||
#endif
|
||||
|
||||
#include <simgear/debug/logstream.hxx>
|
||||
#include <Geometry/poly_support.hxx>
|
||||
|
||||
#include "tgconstruct.hxx"
|
||||
|
||||
|
|
@ -81,7 +80,7 @@ void TGConstruct::calc_normals( std::vector<SGGeod>& geod_nodes, std::vector<SGV
|
|||
SGVec3d v2 = wgs84_nodes[ poly.GetTriIdx( tri, 1 ) ];
|
||||
SGVec3d v3 = wgs84_nodes[ poly.GetTriIdx( tri, 2 ) ];
|
||||
|
||||
area = triangle_area( g1, g2, g3 );
|
||||
area = tgTriangle::area( g1, g2, g3 );
|
||||
normal = calc_normal( area, v1, v2, v3 );
|
||||
|
||||
poly.SetTriFaceArea( tri, area );
|
||||
|
|
|
|||
|
|
@ -29,8 +29,6 @@
|
|||
#include <simgear/misc/sg_dir.hxx>
|
||||
#include <simgear/debug/logstream.hxx>
|
||||
|
||||
#include <Geometry/poly_support.hxx>
|
||||
|
||||
#include "tgconstruct.hxx"
|
||||
|
||||
using std::string;
|
||||
|
|
|
|||
|
|
@ -30,8 +30,6 @@
|
|||
#include <simgear/debug/logstream.hxx>
|
||||
#include <simgear/io/lowlevel.hxx>
|
||||
|
||||
#include <Geometry/poly_support.hxx>
|
||||
|
||||
#include "tgconstruct.hxx"
|
||||
|
||||
using std::string;
|
||||
|
|
@ -134,7 +132,7 @@ void TGConstruct::WriteNeighborFaces( gzFile& fp, const SGGeod& pt ) const
|
|||
SGVec3d const& wgs_p2 = nodes[ poly.GetTriIdx( tri, 0) ].GetWgs84();
|
||||
SGVec3d const& wgs_p3 = nodes[ poly.GetTriIdx( tri, 0) ].GetWgs84();
|
||||
|
||||
double face_area = triangle_area( p1, p2, p3 );
|
||||
double face_area = tgTriangle::area( p1, p2, p3 );
|
||||
SGVec3f face_normal = calc_normal( face_area, wgs_p1, wgs_p2, wgs_p3 );
|
||||
|
||||
sgWriteDouble( fp, face_area );
|
||||
|
|
|
|||
|
|
@ -28,7 +28,7 @@
|
|||
#include <simgear/math/sg_types.hxx>
|
||||
#include <simgear/misc/sgstream.hxx>
|
||||
|
||||
#include <Lib/Geometry/point3d.hxx>
|
||||
#include <Lib/Polygon/point3d.hxx>
|
||||
|
||||
class TGArray {
|
||||
|
||||
|
|
|
|||
|
|
@ -2,7 +2,6 @@ include_directories(${PROJECT_SOURCE_DIR}/src/Lib)
|
|||
|
||||
add_subdirectory(Array)
|
||||
add_subdirectory(DEM)
|
||||
add_subdirectory(Geometry)
|
||||
add_subdirectory(HGT)
|
||||
add_subdirectory(Polygon)
|
||||
add_subdirectory(landcover)
|
||||
|
|
|
|||
|
|
@ -1,22 +0,0 @@
|
|||
include_directories(${GDAL_INCLUDE_DIR})
|
||||
|
||||
add_library(Geometry STATIC
|
||||
contour_tree.cxx
|
||||
line.cxx
|
||||
poly_support.cxx
|
||||
poly_extra.cxx
|
||||
poly_cgal.cxx
|
||||
rectangle.cxx
|
||||
tg_nodes.cxx
|
||||
trinodes.cxx
|
||||
trisegs.cxx
|
||||
contour_tree.hxx
|
||||
line.hxx
|
||||
point3d.hxx
|
||||
poly_support.hxx
|
||||
poly_extra.hxx
|
||||
rectangle.hxx
|
||||
tg_nodes.hxx
|
||||
trinodes.hxx
|
||||
trisegs.hxx
|
||||
)
|
||||
|
|
@ -1,326 +0,0 @@
|
|||
#ifndef JAMA_QR_H
|
||||
#define JAMA_QR_H
|
||||
|
||||
#include "tnt_array1d.h"
|
||||
#include "tnt_array2d.h"
|
||||
#include "tnt_math_utils.h"
|
||||
|
||||
namespace JAMA
|
||||
{
|
||||
|
||||
/**
|
||||
<p>
|
||||
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.
|
||||
<P>
|
||||
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).
|
||||
|
||||
<p>
|
||||
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.
|
||||
|
||||
<p>
|
||||
(Adapted from JAMA, a Java Matrix Library, developed by jointly
|
||||
by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama).
|
||||
*/
|
||||
|
||||
template <class Real>
|
||||
class QR {
|
||||
|
||||
|
||||
/** Array for internal storage of decomposition.
|
||||
@serial internal array storage.
|
||||
*/
|
||||
|
||||
TNT::Array2D<Real> 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<Real> Rdiag;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
/**
|
||||
Create a QR factorization object for A.
|
||||
|
||||
@param A rectangular (m>=n) matrix.
|
||||
*/
|
||||
QR(const TNT::Array2D<Real> &A) /* constructor */
|
||||
{
|
||||
QR_ = A.copy();
|
||||
m = A.dim1();
|
||||
n = A.dim2();
|
||||
Rdiag = TNT::Array1D<Real>(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<Real> getHouseholder (void) const
|
||||
{
|
||||
TNT::Array2D<Real> 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<Real> getR() const
|
||||
{
|
||||
TNT::Array2D<Real> 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<Real> getQ() const
|
||||
{
|
||||
int i=0, j=0, k=0;
|
||||
|
||||
TNT::Array2D<Real> 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<Real> solve(const TNT::Array1D<Real> &b) const
|
||||
{
|
||||
if (b.dim1() != m) /* arrays must be conformant */
|
||||
return TNT::Array1D<Real>();
|
||||
|
||||
if ( !isFullRank() ) /* matrix is rank deficient */
|
||||
{
|
||||
return TNT::Array1D<Real>();
|
||||
}
|
||||
|
||||
TNT::Array1D<Real> 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<Real> x_(n);
|
||||
for (int i=0; i<n; i++)
|
||||
x_[i] = x[i];
|
||||
|
||||
return x_;
|
||||
}
|
||||
|
||||
/** Least squares solution of A*X = B
|
||||
@param B m x k Array (must conform).
|
||||
@return X n x k Array 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 (0x0) array.
|
||||
*/
|
||||
|
||||
TNT::Array2D<Real> solve(const TNT::Array2D<Real> &B) const
|
||||
{
|
||||
if (B.dim1() != m) /* arrays must be conformant */
|
||||
return TNT::Array2D<Real>(0,0);
|
||||
|
||||
if ( !isFullRank() ) /* matrix is rank deficient */
|
||||
{
|
||||
return TNT::Array2D<Real>(0,0);
|
||||
}
|
||||
|
||||
int nx = B.dim2();
|
||||
TNT::Array2D<Real> 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<Real> X_(n,nx);
|
||||
for (i=0; i<n; i++)
|
||||
for (j=0; j<nx; j++)
|
||||
X_[i][j] = X[i][j];
|
||||
|
||||
return X_;
|
||||
}
|
||||
|
||||
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
// namespace JAMA
|
||||
|
||||
#endif
|
||||
// JAMA_QR__H
|
||||
|
||||
|
|
@ -1,301 +0,0 @@
|
|||
/*
|
||||
*
|
||||
* 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_ARRAY1D_H
|
||||
#define TNT_ARRAY1D_H
|
||||
|
||||
//#include <cstdlib>
|
||||
#include <iostream>
|
||||
|
||||
#ifdef TNT_BOUNDS_CHECK
|
||||
#include <assert.h>
|
||||
#endif
|
||||
|
||||
|
||||
#include "tnt_i_refvec.h"
|
||||
|
||||
namespace TNT
|
||||
{
|
||||
|
||||
template <class T>
|
||||
class Array1D
|
||||
{
|
||||
|
||||
private:
|
||||
|
||||
/* ... */
|
||||
i_refvec<T> 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<T> subarray(int i0, int i1);
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
Array1D<T>::Array1D() : v_(), n_(0), data_(0) {}
|
||||
|
||||
template <class T>
|
||||
Array1D<T>::Array1D(const Array1D<T> &A) : v_(A.v_), n_(A.n_),
|
||||
data_(A.data_)
|
||||
{
|
||||
#ifdef TNT_DEBUG
|
||||
std::cout << "Created Array1D(const Array1D<T> &A) \n";
|
||||
#endif
|
||||
|
||||
}
|
||||
|
||||
template <class T>
|
||||
Array1D<T>::Array1D(int n) : v_(n), n_(n), data_(v_.begin())
|
||||
{
|
||||
#ifdef TNT_DEBUG
|
||||
std::cout << "Created Array1D(int n) \n";
|
||||
#endif
|
||||
}
|
||||
|
||||
template <class T>
|
||||
Array1D<T>::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 <class T>
|
||||
Array1D<T>::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 <class T>
|
||||
inline Array1D<T>::operator T*()
|
||||
{
|
||||
return &(v_[0]);
|
||||
}
|
||||
|
||||
|
||||
template <class T>
|
||||
inline Array1D<T>::operator const T*()
|
||||
{
|
||||
return &(v_[0]);
|
||||
}
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
inline T& Array1D<T>::operator[](int i)
|
||||
{
|
||||
#ifdef TNT_BOUNDS_CHECK
|
||||
assert(i>= 0);
|
||||
assert(i < n_);
|
||||
#endif
|
||||
return data_[i];
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline const T& Array1D<T>::operator[](int i) const
|
||||
{
|
||||
#ifdef TNT_BOUNDS_CHECK
|
||||
assert(i>= 0);
|
||||
assert(i < n_);
|
||||
#endif
|
||||
return data_[i];
|
||||
}
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
inline T& Array1D<T>::operator()(int i)
|
||||
{
|
||||
#ifdef TNT_BOUNDS_CHECK
|
||||
assert(i>= 0);
|
||||
assert(i < n_);
|
||||
#endif
|
||||
return data_[i-1];
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline const T& Array1D<T>::operator()(int i) const
|
||||
{
|
||||
#ifdef TNT_BOUNDS_CHECK
|
||||
assert(i>= 0);
|
||||
assert(i < n_);
|
||||
#endif
|
||||
return data_[i-1];
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
Array1D<T> & Array1D<T>::operator=(const T &a)
|
||||
{
|
||||
set_(data_, data_+n_, a);
|
||||
return *this;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
Array1D<T> Array1D<T>::copy() const
|
||||
{
|
||||
Array1D A( n_);
|
||||
copy_(A.data_, data_, n_);
|
||||
|
||||
return A;
|
||||
}
|
||||
|
||||
|
||||
template <class T>
|
||||
Array1D<T> & Array1D<T>::inject(const Array1D &A)
|
||||
{
|
||||
if (A.n_ == n_)
|
||||
copy_(data_, A.data_, n_);
|
||||
|
||||
return *this;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
Array1D<T> & Array1D<T>::ref(const Array1D<T> &A)
|
||||
{
|
||||
if (this != &A)
|
||||
{
|
||||
v_ = A.v_; /* operator= handles the reference counting. */
|
||||
n_ = A.n_;
|
||||
data_ = A.data_;
|
||||
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
Array1D<T> & Array1D<T>::operator=(const Array1D<T> &A)
|
||||
{
|
||||
return ref(A);
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline int Array1D<T>::dim1() const { return n_; }
|
||||
|
||||
template <class T>
|
||||
inline int Array1D<T>::dim() const { return n_; }
|
||||
|
||||
template <class T>
|
||||
Array1D<T>::~Array1D() {}
|
||||
|
||||
|
||||
/* ............................ exented interface ......................*/
|
||||
|
||||
template <class T>
|
||||
inline int Array1D<T>::ref_count() const
|
||||
{
|
||||
return v_.ref_count();
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline Array1D<T> Array1D<T>::subarray(int i0, int i1)
|
||||
{
|
||||
if (((i0 >= 0) && (i1 < n_)) || (i0 <= i1))
|
||||
{
|
||||
Array1D<T> X(*this); /* create a new instance of this array. */
|
||||
X.n_ = i1-i0+1;
|
||||
X.data_ += i0;
|
||||
|
||||
return X;
|
||||
}
|
||||
else
|
||||
{
|
||||
return Array1D<T>();
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* private internal functions */
|
||||
|
||||
|
||||
template <class T>
|
||||
void Array1D<T>::set_(T* begin, T* end, const T& a)
|
||||
{
|
||||
for (T* p=begin; p<end; p++)
|
||||
*p = a;
|
||||
|
||||
}
|
||||
|
||||
template <class T>
|
||||
void Array1D<T>::copy_(T* p, const T* q, int len) const
|
||||
{
|
||||
T *end = p + len;
|
||||
while (p<end )
|
||||
*p++ = *q++;
|
||||
|
||||
}
|
||||
|
||||
|
||||
} /* namespace TNT */
|
||||
|
||||
#endif
|
||||
/* TNT_ARRAY1D_H */
|
||||
|
||||
|
|
@ -1,451 +0,0 @@
|
|||
/*
|
||||
*
|
||||
* 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_ARRAY2D_H
|
||||
#define TNT_ARRAY2D_H
|
||||
|
||||
#include <cstdlib>
|
||||
#include <iostream>
|
||||
#ifdef TNT_BOUNDS_CHECK
|
||||
#include <assert.h>
|
||||
#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..
|
||||
|
||||
<p>
|
||||
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.
|
||||
|
||||
<p>
|
||||
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 T>
|
||||
class Array2D
|
||||
{
|
||||
|
||||
|
||||
private:
|
||||
|
||||
|
||||
|
||||
Array1D<T> data_;
|
||||
Array1D<T*> 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,
|
||||
<pre>
|
||||
template < class ArrayTwoD >
|
||||
void foo(ArrayTwoD &A)
|
||||
{
|
||||
A::value_type first_entry = A[0][0];
|
||||
...
|
||||
}
|
||||
</pre>
|
||||
*/
|
||||
typedef T value_type;
|
||||
|
||||
/**
|
||||
Create a null array. This is <b>not</b> 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<dobule>.
|
||||
|
||||
*/
|
||||
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<dobule>&".
|
||||
|
||||
*/
|
||||
|
||||
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).
|
||||
|
||||
<p>
|
||||
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 <class T>
|
||||
Array2D<T>::Array2D() : data_(), v_(), m_(0), n_(0) {}
|
||||
|
||||
|
||||
template <class T>
|
||||
Array2D<T>::Array2D(const Array2D<T> &A) : data_(A.data_), v_(A.v_),
|
||||
m_(A.m_), n_(A.n_) {}
|
||||
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
Array2D<T>::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<m; i++)
|
||||
{
|
||||
v_[i] = p;
|
||||
p += n;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
Array2D<T>::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<m; i++)
|
||||
{
|
||||
v_[i] = p;
|
||||
p += n;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <class T>
|
||||
Array2D<T>::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<m; i++)
|
||||
{
|
||||
v_[i] = p;
|
||||
p += n;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template <class T>
|
||||
inline T* Array2D<T>::operator[](int i)
|
||||
{
|
||||
#ifdef TNT_BOUNDS_CHECK
|
||||
assert(i >= 0);
|
||||
assert(i < m_);
|
||||
#endif
|
||||
|
||||
return v_[i];
|
||||
|
||||
}
|
||||
|
||||
|
||||
template <class T>
|
||||
inline const T* Array2D<T>::operator[](int i) const
|
||||
{
|
||||
#ifdef TNT_BOUNDS_CHECK
|
||||
assert(i >= 0);
|
||||
assert(i < m_);
|
||||
#endif
|
||||
|
||||
return v_[i];
|
||||
|
||||
}
|
||||
|
||||
template <class T>
|
||||
Array2D<T> & Array2D<T>::operator=(const T &a)
|
||||
{
|
||||
/* non-optimzied, but will work with subarrays in future verions */
|
||||
|
||||
for (int i=0; i<m_; i++)
|
||||
for (int j=0; j<n_; j++)
|
||||
v_[i][j] = a;
|
||||
return *this;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
Array2D<T> Array2D<T>::copy() const
|
||||
{
|
||||
Array2D A(m_, n_);
|
||||
|
||||
for (int i=0; i<m_; i++)
|
||||
for (int j=0; j<n_; j++)
|
||||
A[i][j] = v_[i][j];
|
||||
|
||||
|
||||
return A;
|
||||
}
|
||||
|
||||
|
||||
template <class T>
|
||||
Array2D<T> & Array2D<T>::inject(const Array2D &A)
|
||||
{
|
||||
if (A.m_ == m_ && A.n_ == n_)
|
||||
{
|
||||
for (int i=0; i<m_; i++)
|
||||
for (int j=0; j<n_; j++)
|
||||
v_[i][j] = A[i][j];
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
Array2D<T> & Array2D<T>::ref(const Array2D<T> &A)
|
||||
{
|
||||
if (this != &A)
|
||||
{
|
||||
v_ = A.v_;
|
||||
data_ = A.data_;
|
||||
m_ = A.m_;
|
||||
n_ = A.n_;
|
||||
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
Array2D<T> & Array2D<T>::operator=(const Array2D<T> &A)
|
||||
{
|
||||
return ref(A);
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline int Array2D<T>::dim1() const { return m_; }
|
||||
|
||||
template <class T>
|
||||
inline int Array2D<T>::dim2() const { return n_; }
|
||||
|
||||
|
||||
template <class T>
|
||||
Array2D<T>::~Array2D() {}
|
||||
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
inline Array2D<T>::operator T**()
|
||||
{
|
||||
return &(v_[0]);
|
||||
}
|
||||
template <class T>
|
||||
inline Array2D<T>::operator const T**() const
|
||||
{
|
||||
return static_cast<const T**>(&(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 <class T>
|
||||
Array2D<T> Array2D<T>::subarray(int i0, int i1, int j0, int j1)
|
||||
{
|
||||
Array2D<T> 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<T*>(m);
|
||||
T* p = &(data_[0]) + i0 * n_ + j0;
|
||||
for (int i=0; i<m; i++)
|
||||
{
|
||||
A.v_[i] = p + i*n_;
|
||||
|
||||
}
|
||||
return A;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline int Array2D<T>::ref_count()
|
||||
{
|
||||
return ref_count_data();
|
||||
}
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
inline int Array2D<T>::ref_count_data()
|
||||
{
|
||||
return data_.ref_count();
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline int Array2D<T>::ref_count_dim1()
|
||||
{
|
||||
return v_.ref_count();
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
} /* namespace TNT */
|
||||
|
||||
#endif
|
||||
/* TNT_ARRAY2D_H */
|
||||
|
||||
|
|
@ -1,243 +0,0 @@
|
|||
/*
|
||||
*
|
||||
* 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 <cstdlib>
|
||||
#include <iostream>
|
||||
|
||||
#ifdef TNT_BOUNDS_CHECK
|
||||
#include <assert.h>
|
||||
#endif
|
||||
|
||||
#ifndef NULL
|
||||
#define NULL 0
|
||||
#endif
|
||||
|
||||
namespace TNT
|
||||
{
|
||||
/*
|
||||
Internal representation of ref-counted array. The TNT
|
||||
arrays all use this building block.
|
||||
|
||||
<p>
|
||||
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 T>
|
||||
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<T> & operator=(const i_refvec<T> &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 <class T>
|
||||
void i_refvec<T>::copy_(T* p, const T* q, const T* e)
|
||||
{
|
||||
for (T* t=p; q<e; t++, q++)
|
||||
*t= *q;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
i_refvec<T>::i_refvec() : data_(NULL), ref_count_(NULL) {}
|
||||
|
||||
/**
|
||||
In case n is 0 or negative, it does NOT call new.
|
||||
*/
|
||||
template <class T>
|
||||
i_refvec<T>::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 <class T>
|
||||
inline i_refvec<T>::i_refvec(const i_refvec<T> &V): data_(V.data_),
|
||||
ref_count_(V.ref_count_)
|
||||
{
|
||||
if (V.ref_count_ != NULL)
|
||||
(*(V.ref_count_))++;
|
||||
}
|
||||
|
||||
|
||||
template <class T>
|
||||
i_refvec<T>::i_refvec(T* data) : data_(data), ref_count_(NULL) {}
|
||||
|
||||
template <class T>
|
||||
inline T* i_refvec<T>::begin()
|
||||
{
|
||||
return data_;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline const T& i_refvec<T>::operator[](int i) const
|
||||
{
|
||||
return data_[i];
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline T& i_refvec<T>::operator[](int i)
|
||||
{
|
||||
return data_[i];
|
||||
}
|
||||
|
||||
|
||||
template <class T>
|
||||
inline const T* i_refvec<T>::begin() const
|
||||
{
|
||||
return data_;
|
||||
}
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
i_refvec<T> & i_refvec<T>::operator=(const i_refvec<T> &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 <class T>
|
||||
void i_refvec<T>::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<class T>
|
||||
int i_refvec<T>::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 <class T>
|
||||
int i_refvec<T>::ref_count() const
|
||||
{
|
||||
if (data_ == NULL)
|
||||
return 0;
|
||||
else
|
||||
return (ref_count_ != NULL ? *ref_count_ : -1) ;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
i_refvec<T>::~i_refvec()
|
||||
{
|
||||
if (ref_count_ != NULL)
|
||||
{
|
||||
(*ref_count_)--;
|
||||
|
||||
if (*ref_count_ == 0)
|
||||
destroy();
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
} /* namespace TNT */
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
#endif
|
||||
/* TNT_I_REFVEC_H */
|
||||
|
||||
|
|
@ -1,37 +0,0 @@
|
|||
#ifndef MATH_UTILS_H
|
||||
#define MATH_UTILS_H
|
||||
|
||||
|
||||
/* needed for abs(), sqrt() below */
|
||||
#include <cmath>
|
||||
|
||||
|
||||
|
||||
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 <class Real>
|
||||
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 */
|
||||
|
|
@ -1,49 +0,0 @@
|
|||
// contour_tree.cxx -- routines for building a contour tree showing
|
||||
// which contours are inside if which contours
|
||||
//
|
||||
// Written by Curtis Olson, started June 2000.
|
||||
//
|
||||
// Copyright (C) 2000 Curtis L. Olson - http://www.flightgear.org/~curt
|
||||
//
|
||||
// This program is free software; you can redistribute it and/or
|
||||
// modify it under the terms of the GNU General Public License as
|
||||
// published by the Free Software Foundation; either version 2 of the
|
||||
// License, or (at your option) any later version.
|
||||
//
|
||||
// This program is distributed in the hope that it will be useful, but
|
||||
// WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
// General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
|
||||
//
|
||||
// $Id: contour_tree.cxx,v 1.3 2004-11-19 22:25:50 curt Exp $
|
||||
|
||||
|
||||
#include "contour_tree.hxx"
|
||||
|
||||
|
||||
// Constructor
|
||||
TGContourNode::TGContourNode() {
|
||||
kids.clear();
|
||||
}
|
||||
|
||||
|
||||
// Constructor
|
||||
TGContourNode::TGContourNode( int n ) {
|
||||
contour_num = n;
|
||||
kids.clear();
|
||||
}
|
||||
|
||||
|
||||
// Destructor
|
||||
TGContourNode::~TGContourNode() {
|
||||
for ( unsigned int i = 0; i < kids.size(); ++i ) {
|
||||
if ( kids[i] != NULL ) {
|
||||
delete kids[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -1,81 +0,0 @@
|
|||
// contour_tree.hxx -- routines for building a contour tree showing
|
||||
// which contours are inside if which contours
|
||||
//
|
||||
// Written by Curtis Olson, started June 2000.
|
||||
//
|
||||
// Copyright (C) 2000 Curtis L. Olson - http://www.flightgear.org/~curt
|
||||
//
|
||||
// This program is free software; you can redistribute it and/or
|
||||
// modify it under the terms of the GNU General Public License as
|
||||
// published by the Free Software Foundation; either version 2 of the
|
||||
// License, or (at your option) any later version.
|
||||
//
|
||||
// This program is distributed in the hope that it will be useful, but
|
||||
// WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
// General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
|
||||
//
|
||||
// $Id: contour_tree.hxx,v 1.5 2004-11-19 22:25:50 curt Exp $
|
||||
|
||||
|
||||
#ifndef _CONTOUR_TREE_HXX
|
||||
#define _CONTOUR_TREE_HXX
|
||||
|
||||
|
||||
#ifndef __cplusplus
|
||||
# error This library requires C++
|
||||
#endif
|
||||
|
||||
|
||||
#include <simgear/compiler.h>
|
||||
|
||||
#include <vector>
|
||||
#include <cstddef>
|
||||
|
||||
// forward declaration
|
||||
class TGContourNode;
|
||||
|
||||
typedef std::vector < TGContourNode * > contour_kids;
|
||||
typedef contour_kids::iterator contour_kids_iterator;
|
||||
typedef contour_kids::const_iterator const_contour_kids_iterator;
|
||||
|
||||
|
||||
// a simple class for building a contour tree for a polygon. The
|
||||
// contour tree shows the hierarchy of which contour is inside which
|
||||
// contour.
|
||||
|
||||
class TGContourNode {
|
||||
|
||||
private:
|
||||
|
||||
int contour_num; // -1 for the root node
|
||||
contour_kids kids;
|
||||
|
||||
public:
|
||||
|
||||
TGContourNode();
|
||||
TGContourNode( int n );
|
||||
|
||||
~TGContourNode();
|
||||
|
||||
inline int get_contour_num() const { return contour_num; }
|
||||
inline void set_contour_num( int n ) { contour_num = n; }
|
||||
|
||||
inline int get_num_kids() const { return kids.size(); }
|
||||
inline TGContourNode *get_kid( int n ) const { return kids[n]; }
|
||||
inline void add_kid( TGContourNode *kid ) { kids.push_back( kid ); }
|
||||
inline void remove_kid( int n ) {
|
||||
// cout << "kids[" << n << "] = " << kids[n] << endl;
|
||||
delete kids[n];
|
||||
kids[n] = NULL;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
#endif // _CONTOUR_TREE_HXX
|
||||
|
||||
|
||||
|
|
@ -1,78 +0,0 @@
|
|||
// line.cxx - a simple multi-segment line class.
|
||||
//
|
||||
// Started by David Megginson, July 2002
|
||||
//
|
||||
// This file is in the Public Domain and comes with NO WARRANTY OF ANY KIND.
|
||||
|
||||
|
||||
#include "line.hxx"
|
||||
|
||||
namespace tg {
|
||||
|
||||
Line::Line ()
|
||||
{
|
||||
}
|
||||
|
||||
Line::Line (const Line &l)
|
||||
{
|
||||
int nPoints = l.getPointCount();
|
||||
for (int i = 0; i < nPoints; i++)
|
||||
addPoint(l.getPoint(i));
|
||||
}
|
||||
|
||||
Line::~Line ()
|
||||
{
|
||||
}
|
||||
|
||||
int
|
||||
Line::getPointCount () const
|
||||
{
|
||||
return _points.size();
|
||||
}
|
||||
|
||||
const Point3D &
|
||||
Line::getPoint (int index) const
|
||||
{
|
||||
return _points[index];
|
||||
}
|
||||
|
||||
Point3D &
|
||||
Line::getPoint (int index)
|
||||
{
|
||||
return _points[index];
|
||||
}
|
||||
|
||||
void
|
||||
Line::addPoint (const Point3D &point)
|
||||
{
|
||||
_points.push_back(point);
|
||||
}
|
||||
|
||||
Rectangle
|
||||
Line::getBounds () const
|
||||
{
|
||||
SGGeod min;
|
||||
SGGeod max;
|
||||
|
||||
int nPoints = _points.size();
|
||||
for (int i = 0; i < nPoints; i++) {
|
||||
if (i == 0) {
|
||||
min = max = _points[i].toSGGeod();
|
||||
} else {
|
||||
if (_points[i].x() < min.getLongitudeDeg())
|
||||
min.setLongitudeDeg(_points[i].x());
|
||||
if (_points[i].x() > max.getLongitudeDeg())
|
||||
max.setLongitudeDeg(_points[i].x());
|
||||
if (_points[i].y() < min.getLatitudeDeg())
|
||||
min.setLatitudeDeg(_points[i].y());
|
||||
if (_points[i].y() > max.getLatitudeDeg())
|
||||
max.setLatitudeDeg(_points[i].y());
|
||||
}
|
||||
}
|
||||
|
||||
return Rectangle(min, max);
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
// end of line.cxx
|
||||
|
|
@ -1,92 +0,0 @@
|
|||
// line.hxx - a simple multi-segment line class.
|
||||
//
|
||||
// Started by David Megginson, July 2002
|
||||
//
|
||||
// This file is in the Public Domain and comes with NO WARRANTY OF ANY KIND.
|
||||
|
||||
#ifndef __LINE_HXX
|
||||
#define __LINE_HXX
|
||||
|
||||
#ifndef __cplusplus
|
||||
# error This library requires C++
|
||||
#endif
|
||||
|
||||
#include <simgear/compiler.h>
|
||||
#include <Geometry/point3d.hxx>
|
||||
|
||||
#include <vector>
|
||||
|
||||
#include "rectangle.hxx"
|
||||
|
||||
namespace tg {
|
||||
|
||||
/**
|
||||
* A simple multi-segment line class.
|
||||
*
|
||||
* A line segment is a growable list of points. I will add more
|
||||
* functionality if or when it is needed.
|
||||
*/
|
||||
class Line
|
||||
{
|
||||
public:
|
||||
/**
|
||||
* Create a new line with no points.
|
||||
*/
|
||||
Line ();
|
||||
|
||||
/**
|
||||
* Copy an existing line.
|
||||
*
|
||||
* @param l The line to copy.
|
||||
*/
|
||||
Line (const Line &l);
|
||||
|
||||
/**
|
||||
* Destructor.
|
||||
*/
|
||||
virtual ~Line ();
|
||||
|
||||
/**
|
||||
* Get the number of points currently in the line.
|
||||
*
|
||||
* @return The point count.
|
||||
*/
|
||||
virtual int getPointCount () const;
|
||||
|
||||
/**
|
||||
* Get a point in the line (const).
|
||||
*
|
||||
* @param index The index of the point, zero-based.
|
||||
* @return The point at the index specified.
|
||||
*/
|
||||
virtual const Point3D &getPoint (int index) const;
|
||||
|
||||
/**
|
||||
* Get a point in the line (non-const).
|
||||
*
|
||||
* @param index The index of the point, zero-based.
|
||||
* @return The point at the index specified.
|
||||
*/
|
||||
virtual Point3D &getPoint (int index);
|
||||
|
||||
/**
|
||||
* Add a new point to the end of the line.
|
||||
*
|
||||
* @param point The point to add.
|
||||
*/
|
||||
virtual void addPoint (const Point3D &point);
|
||||
|
||||
/**
|
||||
* Get the bounding rectangle for the line.
|
||||
*
|
||||
* @return The bounding rectangle.
|
||||
*/
|
||||
virtual Rectangle getBounds () const;
|
||||
|
||||
private:
|
||||
std::vector<Point3D> _points;
|
||||
};
|
||||
|
||||
};
|
||||
|
||||
#endif // __LINE_HXX
|
||||
|
|
@ -1,165 +0,0 @@
|
|||
#include <float.h>
|
||||
#include <stdio.h>
|
||||
|
||||
#include <simgear/compiler.h>
|
||||
#include <simgear/constants.h>
|
||||
#include <Geometry/point3d.hxx>
|
||||
#include <simgear/math/sg_types.hxx>
|
||||
#include <simgear/debug/logstream.hxx>
|
||||
#include <simgear/structure/exception.hxx>
|
||||
|
||||
#include <Polygon/polygon.hxx>
|
||||
#include "poly_support.hxx"
|
||||
|
||||
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
|
||||
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
|
||||
#include <CGAL/Triangulation_face_base_with_info_2.h>
|
||||
#include <CGAL/Polygon_2.h>
|
||||
#include <CGAL/Triangle_2.h>
|
||||
#include <iostream>
|
||||
|
||||
/* determining if a face is within the reulting poly */
|
||||
struct FaceInfo2
|
||||
{
|
||||
FaceInfo2() {}
|
||||
int nesting_level;
|
||||
|
||||
bool in_domain(){
|
||||
return nesting_level%2 == 1;
|
||||
}
|
||||
};
|
||||
|
||||
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
|
||||
typedef CGAL::Triangulation_vertex_base_2<K> Vb;
|
||||
typedef CGAL::Triangulation_face_base_with_info_2<FaceInfo2,K> Fbb;
|
||||
typedef CGAL::Constrained_triangulation_face_base_2<K,Fbb> Fb;
|
||||
typedef CGAL::Triangulation_data_structure_2<Vb,Fb> TDS;
|
||||
typedef CGAL::Exact_predicates_tag Itag;
|
||||
typedef CGAL::Constrained_Delaunay_triangulation_2<K, TDS, Itag> CDT;
|
||||
typedef CDT::Point Point;
|
||||
typedef CGAL::Polygon_2<K> Polygon_2;
|
||||
typedef CGAL::Triangle_2<K> Triangle_2;
|
||||
|
||||
void mark_domains(CDT& ct, CDT::Face_handle start, int index, std::list<CDT::Edge>& border )
|
||||
{
|
||||
if(start->info().nesting_level != -1) {
|
||||
return;
|
||||
}
|
||||
|
||||
std::list<CDT::Face_handle> queue;
|
||||
queue.push_back(start);
|
||||
|
||||
while( !queue.empty() ){
|
||||
CDT::Face_handle fh = queue.front();
|
||||
queue.pop_front();
|
||||
if(fh->info().nesting_level == -1) {
|
||||
fh->info().nesting_level = index;
|
||||
for(int i = 0; i < 3; i++) {
|
||||
CDT::Edge e(fh,i);
|
||||
CDT::Face_handle n = fh->neighbor(i);
|
||||
if(n->info().nesting_level == -1) {
|
||||
if(ct.is_constrained(e)) border.push_back(e);
|
||||
else queue.push_back(n);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//explore set of facets connected with non constrained edges,
|
||||
//and attribute to each such set a nesting level.
|
||||
//We start from facets incident to the infinite vertex, with a nesting
|
||||
//level of 0. Then we recursively consider the non-explored facets incident
|
||||
//to constrained edges bounding the former set and increase the nesting level by 1.
|
||||
//Facets in the domain are those with an odd nesting level.
|
||||
void mark_domains(CDT& cdt)
|
||||
{
|
||||
for(CDT::All_faces_iterator it = cdt.all_faces_begin(); it != cdt.all_faces_end(); ++it){
|
||||
it->info().nesting_level = -1;
|
||||
}
|
||||
|
||||
int index = 0;
|
||||
std::list<CDT::Edge> border;
|
||||
mark_domains(cdt, cdt.infinite_face(), index++, border);
|
||||
while(! border.empty()) {
|
||||
CDT::Edge e = border.front();
|
||||
border.pop_front();
|
||||
CDT::Face_handle n = e.first->neighbor(e.second);
|
||||
if(n->info().nesting_level == -1) {
|
||||
mark_domains(cdt, n, e.first->info().nesting_level+1, border);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void insert_polygon(CDT& cdt,const Polygon_2& polygon)
|
||||
{
|
||||
if ( polygon.is_empty() ) return;
|
||||
|
||||
CDT::Vertex_handle v_prev=cdt.insert(*CGAL::cpp0x::prev(polygon.vertices_end()));
|
||||
for (Polygon_2::Vertex_iterator vit=polygon.vertices_begin(); vit!=polygon.vertices_end();++vit) {
|
||||
CDT::Vertex_handle vh=cdt.insert(*vit);
|
||||
cdt.insert_constraint(vh,v_prev);
|
||||
v_prev=vh;
|
||||
}
|
||||
}
|
||||
|
||||
TGPolygon polygon_tesselate_alt_with_extra_cgal( TGPolygon &p, const std::vector<SGGeod>& extra_nodes, bool verbose ) {
|
||||
TGPolygon result;
|
||||
CDT cdt;
|
||||
|
||||
result.erase();
|
||||
|
||||
// Bail right away if polygon is empty
|
||||
if ( p.contours() == 0 ) {
|
||||
return result;
|
||||
}
|
||||
|
||||
// First, insert the extra points
|
||||
std::vector<Point> points;
|
||||
points.reserve(extra_nodes.size());
|
||||
for (unsigned int n = 0; n < extra_nodes.size(); n++) {
|
||||
points.push_back( Point(extra_nodes[n].getLongitudeDeg(), extra_nodes[n].getLatitudeDeg()) );
|
||||
}
|
||||
cdt.insert(points.begin(), points.end());
|
||||
|
||||
// then insert each polygon as a constraint into the triangulation
|
||||
for (int c = 0; c < p.contours(); c++) {
|
||||
point_list contour = p.get_contour( c );
|
||||
Polygon_2 poly;
|
||||
|
||||
for (unsigned int n = 0; n < contour.size(); n++ ) {
|
||||
Point3D node = contour[n];
|
||||
poly.push_back( Point( node.x(), node.y()) );
|
||||
}
|
||||
|
||||
insert_polygon(cdt, poly);
|
||||
}
|
||||
|
||||
mark_domains( cdt );
|
||||
|
||||
int count=0;
|
||||
for (CDT::Finite_faces_iterator fit=cdt.finite_faces_begin(); fit!=cdt.finite_faces_end();++fit) {
|
||||
if ( fit->info().in_domain() ) {
|
||||
Triangle_2 tri = cdt.triangle(fit);
|
||||
|
||||
Point3D p0 = Point3D( tri.vertex(0).x(), tri.vertex(0).y(), 0.0f );
|
||||
Point3D p1 = Point3D( tri.vertex(1).x(), tri.vertex(1).y(), 0.0f );
|
||||
Point3D p2 = Point3D( tri.vertex(2).x(), tri.vertex(2).y(), 0.0f );
|
||||
|
||||
result.add_node( count, p0 );
|
||||
result.add_node( count, p1 );
|
||||
result.add_node( count, p2 );
|
||||
|
||||
++count;
|
||||
// create a contour in result with this face
|
||||
}
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
TGPolygon polygon_tesselate_alt_cgal( TGPolygon &p, bool verbose ) {
|
||||
|
||||
std::vector<SGGeod> pl; pl.clear();
|
||||
return ( polygon_tesselate_alt_with_extra_cgal(p, pl, verbose) );
|
||||
}
|
||||
|
|
@ -1,177 +0,0 @@
|
|||
// poly_extra.cxx -- Extra polygon manipulation routines
|
||||
//
|
||||
// Written by Curtis Olson, started February 2002.
|
||||
//
|
||||
// Copyright (C) 2002 Curtis L. Olson - http://www.flightgear.org/~curt
|
||||
//
|
||||
// This program is free software; you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation; either version 2 of the License, or
|
||||
// (at your option) any later version.
|
||||
//
|
||||
// This program is distributed in the hope that it will be useful,
|
||||
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
// GNU General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
|
||||
//
|
||||
// $Id: poly_extra.cxx,v 1.9 2004-11-19 22:25:49 curt Exp $
|
||||
//
|
||||
|
||||
#include <stdio.h>
|
||||
|
||||
#include <simgear/compiler.h>
|
||||
#include <simgear/debug/logstream.hxx>
|
||||
|
||||
#include <Geometry/poly_support.hxx>
|
||||
|
||||
#include "poly_extra.hxx"
|
||||
|
||||
|
||||
// Divide segment if there are other existing points on it, return the
|
||||
// new polygon
|
||||
void add_intermediate_nodes( int contour, const Point3D& start,
|
||||
const Point3D& end, point_list& tmp_nodes,
|
||||
TGPolygon *result,
|
||||
const double bbEpsilon,
|
||||
const double errEpsilon
|
||||
)
|
||||
{
|
||||
Point3D new_pt;
|
||||
|
||||
SG_LOG(SG_GENERAL, SG_BULK, " " << start << " <==> " << end );
|
||||
|
||||
bool found_extra = find_intermediate_node( start, end, tmp_nodes, &new_pt, bbEpsilon, errEpsilon );
|
||||
|
||||
if ( found_extra ) {
|
||||
// recurse with two sub segments
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, "dividing " << p0 << " " << nodes[extra_index]
|
||||
// << " " << p1);
|
||||
add_intermediate_nodes( contour, start, new_pt, tmp_nodes, result, bbEpsilon, errEpsilon );
|
||||
|
||||
result->add_node( contour, new_pt );
|
||||
SG_LOG(SG_GENERAL, SG_BULK, " adding = " << new_pt);
|
||||
|
||||
add_intermediate_nodes( contour, new_pt, end, tmp_nodes, result, bbEpsilon, errEpsilon );
|
||||
}
|
||||
}
|
||||
|
||||
void add_intermediate_nodes( int contour, const Point3D& start,
|
||||
const Point3D& end, std::vector<SGGeod>&tmp_nodes,
|
||||
TGPolygon *result,
|
||||
const double bbEpsilon,
|
||||
const double errEpsilon
|
||||
)
|
||||
{
|
||||
Point3D new_pt;
|
||||
|
||||
SG_LOG(SG_GENERAL, SG_BULK, " " << start << " <==> " << end );
|
||||
|
||||
bool found_extra = find_intermediate_node( start, end, tmp_nodes, &new_pt, bbEpsilon, errEpsilon );
|
||||
|
||||
if ( found_extra ) {
|
||||
// recurse with two sub segments
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, "dividing " << p0 << " " << nodes[extra_index]
|
||||
// << " " << p1);
|
||||
add_intermediate_nodes( contour, start, new_pt, tmp_nodes, result, bbEpsilon, errEpsilon );
|
||||
|
||||
result->add_node( contour, new_pt );
|
||||
SG_LOG(SG_GENERAL, SG_BULK, " adding = " << new_pt);
|
||||
|
||||
add_intermediate_nodes( contour, new_pt, end, tmp_nodes, result, bbEpsilon, errEpsilon );
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Search each segment for additional vertex points that may have been
|
||||
// created elsewhere that lie on the segment and split it there to
|
||||
// avoid "T" intersections.
|
||||
|
||||
TGPolygon add_nodes_to_poly( const TGPolygon& poly,
|
||||
const TGTriNodes& nodes ) {
|
||||
int i, j;
|
||||
TGPolygon result; result.erase();
|
||||
point_list tmp_nodes = nodes.get_node_list();
|
||||
Point3D p0, p1;
|
||||
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, "add_nodes_to_poly");
|
||||
|
||||
for ( i = 0; i < poly.contours(); ++i ) {
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, "contour = " << i);
|
||||
for ( j = 0; j < poly.contour_size(i) - 1; ++j ) {
|
||||
p0 = poly.get_pt( i, j );
|
||||
p1 = poly.get_pt( i, j + 1 );
|
||||
|
||||
// add start of segment
|
||||
result.add_node( i, p0 );
|
||||
|
||||
// add intermediate points
|
||||
add_intermediate_nodes( i, p0, p1, tmp_nodes, &result, SG_EPSILON*10, SG_EPSILON*4 );
|
||||
|
||||
// end of segment is beginning of next segment
|
||||
}
|
||||
p0 = poly.get_pt( i, poly.contour_size(i) - 1 );
|
||||
p1 = poly.get_pt( i, 0 );
|
||||
|
||||
// add start of segment
|
||||
result.add_node( i, p0 );
|
||||
|
||||
// add intermediate points
|
||||
add_intermediate_nodes( i, p0, p1, tmp_nodes, &result, SG_EPSILON*10, SG_EPSILON*4 );
|
||||
|
||||
// end of segment is beginning of next segment
|
||||
// 5/9/2000 CLO - this results in duplicating the last point
|
||||
// of a contour so I have removed this line.
|
||||
// result.add_node( i, p1 );
|
||||
|
||||
// maintain original hole flag setting
|
||||
result.set_hole_flag( i, poly.get_hole_flag( i ) );
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
TGPolygon add_tgnodes_to_poly( const TGPolygon& poly,
|
||||
const TGNodes* nodes ) {
|
||||
TGPolygon result; result.erase();
|
||||
SGGeod min, max;
|
||||
Point3D p0, p1;
|
||||
std::vector<SGGeod> poly_points;
|
||||
|
||||
poly.get_bounding_box(min, max);
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "add_tgnodes_to_poly : min " << min << " max " << max );
|
||||
|
||||
nodes->get_geod_inside( min, max, poly_points );
|
||||
|
||||
for ( int i = 0; i < poly.contours(); ++i ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "contour = " << i);
|
||||
for ( int j = 0; j < poly.contour_size(i) - 1; ++j ) {
|
||||
p0 = poly.get_pt( i, j );
|
||||
p1 = poly.get_pt( i, j + 1 );
|
||||
|
||||
// add start of segment
|
||||
result.add_node( i, p0 );
|
||||
|
||||
// add intermediate points
|
||||
add_intermediate_nodes( i, p0, p1, poly_points, &result, SG_EPSILON*10, SG_EPSILON*4 );
|
||||
|
||||
// end of segment is beginning of next segment
|
||||
}
|
||||
p0 = poly.get_pt( i, poly.contour_size(i) - 1 );
|
||||
p1 = poly.get_pt( i, 0 );
|
||||
|
||||
// add start of segment
|
||||
result.add_node( i, p0 );
|
||||
|
||||
// add intermediate points
|
||||
add_intermediate_nodes( i, p0, p1, poly_points, &result, SG_EPSILON*10, SG_EPSILON*4 );
|
||||
|
||||
// maintain original hole flag setting
|
||||
result.set_hole_flag( i, poly.get_hole_flag( i ) );
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
|
@ -1,66 +0,0 @@
|
|||
// poly_extra.hxx -- Extra polygon manipulation routines
|
||||
//
|
||||
// Written by Curtis Olson, started February 2002.
|
||||
//
|
||||
// Copyright (C) 2002 Curtis L. Olson - http://www.flightgear.org/~curt
|
||||
//
|
||||
// This program is free software; you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation; either version 2 of the License, or
|
||||
// (at your option) any later version.
|
||||
//
|
||||
// This program is distributed in the hope that it will be useful,
|
||||
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
// GNU General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
|
||||
//
|
||||
// $Id: poly_extra.hxx,v 1.5 2004-11-19 22:25:49 curt Exp $
|
||||
//
|
||||
|
||||
|
||||
#ifndef _POLY_EXTRA_HXX
|
||||
#define _POLY_EXTRA_HXX
|
||||
|
||||
|
||||
#include <Geometry/point3d.hxx>
|
||||
|
||||
#include <Geometry/trinodes.hxx>
|
||||
#include <Geometry/tg_nodes.hxx>
|
||||
|
||||
#include <Polygon/polygon.hxx>
|
||||
|
||||
|
||||
// Divide segment if there are other existing points on it, return the
|
||||
// new polygon
|
||||
void add_intermediate_nodes( int contour, const Point3D& start,
|
||||
const Point3D& end, const point_list& tmp_nodes,
|
||||
TGPolygon *result,
|
||||
const double bbEpsilon = SG_EPSILON*10,
|
||||
const double errEpsilon = SG_EPSILON*4
|
||||
);
|
||||
|
||||
|
||||
// TEMP - converting tgnodes to SGGeod/SGVec3d
|
||||
void add_intermediate_nodes( int contour, const Point3D& start,
|
||||
const Point3D& end, const std::vector<SGGeod>& tmp_nodes,
|
||||
TGPolygon *result,
|
||||
const double bbEpsilon = SG_EPSILON*10,
|
||||
const double errEpsilon = SG_EPSILON*4
|
||||
);
|
||||
|
||||
// Search each segment for additional vertex points that may have been
|
||||
// created elsewhere that lie on the segment and split it there to
|
||||
// avoid "T" intersections.
|
||||
|
||||
TGPolygon add_nodes_to_poly( const TGPolygon& poly,
|
||||
const TGTriNodes& tmp_nodes );
|
||||
|
||||
TGPolygon add_tgnodes_to_poly( const TGPolygon& poly,
|
||||
const TGNodes* nodes );
|
||||
|
||||
|
||||
#endif // _POLY_EXTRA_HXX
|
||||
File diff suppressed because it is too large
Load diff
|
|
@ -1,122 +0,0 @@
|
|||
// poly_support.hxx -- additional supporting routines for the TGPolygon class
|
||||
// specific to the object building process.
|
||||
//
|
||||
// Written by Curtis Olson, started October 1999.
|
||||
//
|
||||
// Copyright (C) 1999 Curtis L. Olson - http://www.flightgear.org/~curt
|
||||
//
|
||||
// This program is free software; you can redistribute it and/or
|
||||
// modify it under the terms of the GNU General Public License as
|
||||
// published by the Free Software Foundation; either version 2 of the
|
||||
// License, or (at your option) any later version.
|
||||
//
|
||||
// This program is distributed in the hope that it will be useful, but
|
||||
// WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
// General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
|
||||
//
|
||||
// $Id: poly_support.hxx,v 1.14 2004-11-19 22:25:50 curt Exp $
|
||||
|
||||
|
||||
#ifndef _POLY_SUPPORT_HXX
|
||||
#define _POLY_SUPPORT_HXX
|
||||
|
||||
|
||||
#ifndef __cplusplus
|
||||
# error This library requires C++
|
||||
#endif
|
||||
|
||||
|
||||
#include <simgear/compiler.h>
|
||||
#include <simgear/math/sg_types.hxx>
|
||||
|
||||
#include <Polygon/polygon.hxx>
|
||||
|
||||
#include "trieles.hxx"
|
||||
#include "trinodes.hxx"
|
||||
|
||||
|
||||
// Calculate the area of a triangle
|
||||
inline double triangle_area( const SGGeod& p1, const SGGeod& p2, const SGGeod& p3 )
|
||||
{
|
||||
return fabs(0.5 * ( p1.getLongitudeDeg() * p2.getLatitudeDeg() - p2.getLongitudeDeg() * p1.getLatitudeDeg() +
|
||||
p2.getLongitudeDeg() * p3.getLatitudeDeg() - p3.getLongitudeDeg() * p2.getLatitudeDeg() +
|
||||
p3.getLongitudeDeg() * p1.getLatitudeDeg() - p1.getLongitudeDeg() * p3.getLatitudeDeg() ));
|
||||
}
|
||||
|
||||
|
||||
// Alternate basic triangulation of a polygon with out adding points
|
||||
// or splitting edges and without regard for holes. Returns a polygon
|
||||
// with one contour per tesselated triangle.
|
||||
TGPolygon polygon_tesselate_alt_with_extra_cgal( TGPolygon &p,
|
||||
const std::vector<SGGeod>& extra_nodes, bool verbose );
|
||||
|
||||
TGPolygon polygon_tesselate_alt_cgal( TGPolygon &p, bool verbose );
|
||||
|
||||
// calculate some "arbitrary" point inside each of the polygons contours
|
||||
void calc_points_inside( TGPolygon& p );
|
||||
|
||||
// snap all points to a grid
|
||||
TGPolygon snap( const TGPolygon &poly, double grid_size );
|
||||
|
||||
// remove duplicate nodes in a polygon should they exist. Returns the
|
||||
// fixed polygon
|
||||
TGPolygon remove_dups( const TGPolygon &poly );
|
||||
|
||||
|
||||
// Search each segment of each contour for degenerate points (i.e. out
|
||||
// of order points that lie coincident on other segments
|
||||
TGPolygon reduce_degeneracy( const TGPolygon& poly );
|
||||
|
||||
|
||||
// Occasionally the outline of the clipped polygon can take a side
|
||||
// track, then double back on return to the start of the side track
|
||||
// branch and continue normally. Attempt to detect and clear this
|
||||
// extraneous nonsense.
|
||||
TGPolygon remove_cycles( const TGPolygon& poly );
|
||||
|
||||
// Occasionally the outline of the clipped polygon can have long spikes
|
||||
// that come close to doubling back on the same segment - this kills
|
||||
// triangulation
|
||||
TGPolygon remove_spikes( const TGPolygon& poly );
|
||||
|
||||
|
||||
// Find a point in the given node list that lies between start and
|
||||
// end, return true if something found, false if nothing found.
|
||||
bool find_intermediate_node( const Point3D& start, const Point3D& end,
|
||||
const point_list& nodes, Point3D *result,
|
||||
const double bbEpsilon = SG_EPSILON*10,
|
||||
const double errEpsilon = SG_EPSILON*4
|
||||
);
|
||||
|
||||
// TEMP
|
||||
bool find_intermediate_node( const Point3D& start, const Point3D& end,
|
||||
const std::vector<SGGeod>& nodes, Point3D *result,
|
||||
const double bbEpsilon = SG_EPSILON*10,
|
||||
const double errEpsilon = SG_EPSILON*4
|
||||
);
|
||||
|
||||
// remove any degenerate contours
|
||||
TGPolygon remove_bad_contours( const TGPolygon &poly );
|
||||
|
||||
// remove any too small contours
|
||||
TGPolygon remove_tiny_contours( const TGPolygon &poly );
|
||||
TGPolygon remove_small_contours( const TGPolygon &poly );
|
||||
|
||||
|
||||
// Write Polygons to Shapefile Support
|
||||
void tgShapefileInit( void );
|
||||
void* tgShapefileOpenDatasource( const char* datasource_name );
|
||||
void* tgShapefileOpenLayer( void* ds_id, const char* layer_name );
|
||||
void tgShapefileCreateFeature( void* ds_id, void* l_id, const TGPolygon &poly, const char* feature_name, bool has_point_inside = false );
|
||||
void tgShapefileCloseLayer( void* l_id );
|
||||
void* tgShapefileCloseDatasource( void* ds_id );
|
||||
|
||||
|
||||
#endif // _POLY_SUPPORT_HXX
|
||||
|
||||
|
||||
|
|
@ -1,72 +0,0 @@
|
|||
// trieles.hxx -- "Triangle" element management class
|
||||
//
|
||||
// Written by Curtis Olson, started March 1999.
|
||||
//
|
||||
// Copyright (C) 1999 Curtis L. Olson - http://www.flightgear.org/~curt
|
||||
//
|
||||
// This program is free software; you can redistribute it and/or
|
||||
// modify it under the terms of the GNU General Public License as
|
||||
// published by the Free Software Foundation; either version 2 of the
|
||||
// License, or (at your option) any later version.
|
||||
//
|
||||
// This program is distributed in the hope that it will be useful, but
|
||||
// WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
// General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
|
||||
//
|
||||
// $Id: trieles.hxx,v 1.4 2004-11-19 22:25:50 curt Exp $
|
||||
|
||||
|
||||
#ifndef _TRIELES_HXX
|
||||
#define _TRIELES_HXX
|
||||
|
||||
|
||||
#ifndef __cplusplus
|
||||
# error This library requires C++
|
||||
#endif
|
||||
|
||||
|
||||
#include <simgear/compiler.h>
|
||||
|
||||
#include <vector>
|
||||
|
||||
// a segment is two integer pointers into the node list
|
||||
class TGTriEle {
|
||||
int n1, n2, n3;
|
||||
|
||||
double attribute;
|
||||
|
||||
public:
|
||||
|
||||
// Constructor and destructor
|
||||
inline TGTriEle( void ) { };
|
||||
inline TGTriEle( int i1, int i2, int i3, double a ) {
|
||||
n1 = i1; n2 = i2; n3 = i3; attribute = a;
|
||||
}
|
||||
|
||||
inline ~TGTriEle( void ) { };
|
||||
|
||||
inline int get_n1() const { return n1; }
|
||||
inline void set_n1( int i ) { n1 = i; }
|
||||
inline int get_n2() const { return n2; }
|
||||
inline void set_n2( int i ) { n2 = i; }
|
||||
inline int get_n3() const { return n3; }
|
||||
inline void set_n3( int i ) { n3 = i; }
|
||||
|
||||
inline double get_attribute() const { return attribute; }
|
||||
inline void set_attribute( double a ) { attribute = a; }
|
||||
};
|
||||
|
||||
|
||||
typedef std::vector < TGTriEle > triele_list;
|
||||
typedef triele_list::iterator triele_list_iterator;
|
||||
typedef triele_list::const_iterator const_triele_list_iterator;
|
||||
|
||||
|
||||
#endif // _TRIELES_HXX
|
||||
|
||||
|
||||
|
|
@ -1,209 +0,0 @@
|
|||
// trisegs.cxx -- "Triangle" segment management class
|
||||
//
|
||||
// Written by Curtis Olson, started March 1999.
|
||||
//
|
||||
// Copyright (C) 1999 Curtis L. Olson - http://www.flightgear.org/~curt
|
||||
//
|
||||
// This program is free software; you can redistribute it and/or
|
||||
// modify it under the terms of the GNU General Public License as
|
||||
// published by the Free Software Foundation; either version 2 of the
|
||||
// License, or (at your option) any later version.
|
||||
//
|
||||
// This program is distributed in the hope that it will be useful, but
|
||||
// WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
// General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
|
||||
//
|
||||
// $Id: trisegs.cxx,v 1.10 2004-11-19 22:25:50 curt Exp $
|
||||
|
||||
#include <simgear/compiler.h>
|
||||
#include <simgear/constants.h>
|
||||
#include <simgear/debug/logstream.hxx>
|
||||
|
||||
#include <Geometry/point3d.hxx>
|
||||
|
||||
#include <iostream>
|
||||
|
||||
#include "trinodes.hxx"
|
||||
|
||||
#include "trisegs.hxx"
|
||||
|
||||
// Constructor
|
||||
TGTriSegments::TGTriSegments( void ) {
|
||||
}
|
||||
|
||||
|
||||
// Destructor
|
||||
TGTriSegments::~TGTriSegments( void ) {
|
||||
}
|
||||
|
||||
|
||||
// Add a segment to the segment list if it doesn't already exist.
|
||||
// Returns the index (starting at zero) of the segment in the list.
|
||||
int TGTriSegments::unique_add( const TGTriSeg& s )
|
||||
{
|
||||
triseg_list_iterator current, last;
|
||||
int counter = 0;
|
||||
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, s.get_n1() << "," << s.get_n2() );
|
||||
|
||||
// check if segment has duplicated endpoints
|
||||
if ( s.get_n1() == s.get_n2() ) {
|
||||
SG_LOG(SG_GENERAL, SG_ALERT, "WARNING: ignoring null segment with the same point for both endpoints" );
|
||||
return -1;
|
||||
}
|
||||
|
||||
// check if segment already exists
|
||||
current = seg_list.begin();
|
||||
last = seg_list.end();
|
||||
for ( ; current != last; ++current ) {
|
||||
if ( s == *current ) {
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, "found an existing segment match" );
|
||||
return counter;
|
||||
}
|
||||
|
||||
++counter;
|
||||
}
|
||||
|
||||
// add to list
|
||||
seg_list.push_back( s );
|
||||
|
||||
return counter;
|
||||
}
|
||||
|
||||
|
||||
// Divide segment if there are other points on it, return the divided
|
||||
// list of segments
|
||||
void TGTriSegments::unique_divide_and_add( const point_list& nodes,
|
||||
const TGTriSeg& s )
|
||||
{
|
||||
Point3D p0 = nodes[ s.get_n1() ];
|
||||
Point3D p1 = nodes[ s.get_n2() ];
|
||||
|
||||
bool found_extra = false;
|
||||
int extra_index = 0;
|
||||
int counter;
|
||||
double m, m1, b, b1, y_err, x_err, y_err_min, x_err_min;
|
||||
const_point_list_iterator current, last;
|
||||
|
||||
// bool temp = false;
|
||||
// if ( s == TGTriSeg( 170, 206 ) ) {
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, "this is it!" );
|
||||
// temp = true;
|
||||
// }
|
||||
|
||||
double xdist = fabs(p0.x() - p1.x());
|
||||
double ydist = fabs(p0.y() - p1.y());
|
||||
x_err_min = xdist + 1.0;
|
||||
y_err_min = ydist + 1.0;
|
||||
|
||||
if ( xdist > ydist ) {
|
||||
// use y = mx + b
|
||||
|
||||
// sort these in a sensible order
|
||||
if ( p0.x() > p1.x() ) {
|
||||
Point3D tmp = p0;
|
||||
p0 = p1;
|
||||
p1 = tmp;
|
||||
}
|
||||
|
||||
m = (p0.y() - p1.y()) / (p0.x() - p1.x());
|
||||
b = p1.y() - m * p1.x();
|
||||
|
||||
// if ( temp ) {
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, "m = " << m << " b = " << b );
|
||||
// }
|
||||
|
||||
current = nodes.begin();
|
||||
last = nodes.end();
|
||||
counter = 0;
|
||||
for ( ; current != last; ++current ) {
|
||||
if ( (current->x() > (p0.x() + SG_EPSILON))
|
||||
&& (current->x() < (p1.x() - SG_EPSILON)) ) {
|
||||
|
||||
// if ( temp ) {
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, counter );
|
||||
// }
|
||||
|
||||
y_err = fabs(current->y() - (m * current->x() + b));
|
||||
|
||||
if ( y_err < FG_PROXIMITY_EPSILON ) {
|
||||
//SG_LOG(SG_GENERAL, SG_DEBUG, "FOUND EXTRA SEGMENT NODE (Y)" );
|
||||
//SG_LOG(SG_GENERAL, SG_DEBUG, p0 << " < " << *current << " < " << p1 );
|
||||
found_extra = true;
|
||||
if ( y_err < y_err_min ) {
|
||||
extra_index = counter;
|
||||
y_err_min = y_err;
|
||||
}
|
||||
}
|
||||
}
|
||||
++counter;
|
||||
}
|
||||
} else {
|
||||
// use x = m1 * y + b1
|
||||
|
||||
// sort these in a sensible order
|
||||
if ( p0.y() > p1.y() ) {
|
||||
Point3D tmp = p0;
|
||||
p0 = p1;
|
||||
p1 = tmp;
|
||||
}
|
||||
|
||||
m1 = (p0.x() - p1.x()) / (p0.y() - p1.y());
|
||||
b1 = p1.x() - m1 * p1.y();
|
||||
|
||||
// bool temp = true;
|
||||
// if ( temp ) {
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, "xdist = " << xdist << " ydist = " << ydist );
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, " p0 = " << p0 << " p1 = " << p1 );
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, " m1 = " << m1 << " b1 = " << b1 );
|
||||
// }
|
||||
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, " should = 0 = " << fabs(p0.x() - (m1 * p0.y() + b1)) );
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, " should = 0 = " << fabs(p1.x() - (m1 * p1.y() + b1)) );
|
||||
|
||||
current = nodes.begin();
|
||||
last = nodes.end();
|
||||
counter = 0;
|
||||
for ( ; current != last; ++current ) {
|
||||
if ( (current->y() > (p0.y() + SG_EPSILON))
|
||||
&& (current->y() < (p1.y() - SG_EPSILON)) ) {
|
||||
|
||||
x_err = fabs(current->x() - (m1 * current->y() + b1));
|
||||
|
||||
// if ( temp ) {
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, " (" << counter << ") x_err = " << x_err );
|
||||
// }
|
||||
|
||||
if ( x_err < FG_PROXIMITY_EPSILON ) {
|
||||
//SG_LOG(SG_GENERAL, SG_DEBUG, "FOUND EXTRA SEGMENT NODE (X)" );
|
||||
//SG_LOG(SG_GENERAL, SG_DEBUG, p0 << " < " << *current << " < " << p1 );
|
||||
found_extra = true;
|
||||
if ( x_err < x_err_min ) {
|
||||
extra_index = counter;
|
||||
x_err_min = x_err;
|
||||
}
|
||||
}
|
||||
}
|
||||
++counter;
|
||||
}
|
||||
}
|
||||
|
||||
if ( found_extra ) {
|
||||
// recurse with two sub segments
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "dividing " << s.get_n1() << " " << extra_index << " " << s.get_n2() );
|
||||
unique_divide_and_add( nodes, TGTriSeg( s.get_n1(), extra_index,
|
||||
s.get_boundary_marker() ) );
|
||||
unique_divide_and_add( nodes, TGTriSeg( extra_index, s.get_n2(),
|
||||
s.get_boundary_marker() ) );
|
||||
} else {
|
||||
// this segment does not need to be divided, lets add it
|
||||
unique_add( s );
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -1,122 +0,0 @@
|
|||
// trisegs.hxx -- "Triangle" segment management class
|
||||
//
|
||||
// Written by Curtis Olson, started March 1999.
|
||||
//
|
||||
// Copyright (C) 1999 Curtis L. Olson - http://www.flightgear.org/~curt
|
||||
//
|
||||
// This program is free software; you can redistribute it and/or
|
||||
// modify it under the terms of the GNU General Public License as
|
||||
// published by the Free Software Foundation; either version 2 of the
|
||||
// License, or (at your option) any later version.
|
||||
//
|
||||
// This program is distributed in the hope that it will be useful, but
|
||||
// WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
// General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
|
||||
//
|
||||
// $Id: trisegs.hxx,v 1.4 2004/11/19 22:25:50 curt Exp $
|
||||
|
||||
|
||||
#ifndef _TRISEGS_HXX
|
||||
#define _TRISEGS_HXX
|
||||
|
||||
|
||||
#ifndef __cplusplus
|
||||
# error This library requires C++
|
||||
#endif
|
||||
|
||||
|
||||
#include <simgear/compiler.h>
|
||||
|
||||
#include <vector>
|
||||
|
||||
#include "trinodes.hxx"
|
||||
|
||||
// a segment is two integer pointers into the node list
|
||||
class TGTriSeg {
|
||||
int n1, n2; // indices into point list
|
||||
int boundary_marker; // flag if segment is a boundary
|
||||
// (i.e. shouldn't get split)
|
||||
|
||||
public:
|
||||
|
||||
// Constructor and destructor
|
||||
inline TGTriSeg( void ) { };
|
||||
inline TGTriSeg( int i1, int i2, int b ) {
|
||||
n1 = i1;
|
||||
n2 = i2;
|
||||
boundary_marker = b;
|
||||
}
|
||||
|
||||
inline ~TGTriSeg( void ) { };
|
||||
|
||||
inline int get_n1() const { return n1; }
|
||||
inline void set_n1( int i ) { n1 = i; }
|
||||
inline int get_n2() const { return n2; }
|
||||
inline void set_n2( int i ) { n2 = i; }
|
||||
inline int get_boundary_marker() const { return boundary_marker; }
|
||||
inline void set_boundary_marker( int b ) { boundary_marker = b; }
|
||||
|
||||
friend bool operator == (const TGTriSeg& a, const TGTriSeg& b);
|
||||
|
||||
};
|
||||
|
||||
inline bool operator == (const TGTriSeg& a, const TGTriSeg& b)
|
||||
{
|
||||
return ((a.n1 == b.n1) && (a.n2 == b.n2))
|
||||
|| ((a.n1 == b.n2) && (a.n2 == b.n1));
|
||||
}
|
||||
|
||||
|
||||
typedef std::vector < TGTriSeg > triseg_list;
|
||||
typedef triseg_list::iterator triseg_list_iterator;
|
||||
typedef triseg_list::const_iterator const_triseg_list_iterator;
|
||||
|
||||
|
||||
class TGTriSegments {
|
||||
|
||||
private:
|
||||
|
||||
triseg_list seg_list;
|
||||
|
||||
// Divide segment if there are other points on it, return the
|
||||
// divided list of segments
|
||||
triseg_list divide_segment( const point_list& nodes,
|
||||
const TGTriSeg& s );
|
||||
|
||||
public:
|
||||
|
||||
// Constructor and destructor
|
||||
TGTriSegments( void );
|
||||
~TGTriSegments( void );
|
||||
|
||||
// delete all the data out of seg_list
|
||||
inline void clear() { seg_list.clear(); }
|
||||
|
||||
// Add a segment to the segment list if it doesn't already exist.
|
||||
// Returns the index (starting at zero) of the segment in the
|
||||
// list.
|
||||
int unique_add( const TGTriSeg& s );
|
||||
|
||||
// Add a segment to the segment list if it doesn't already exist.
|
||||
// Returns the index (starting at zero) of the segment in the list.
|
||||
void unique_divide_and_add( const point_list& node_list,
|
||||
const TGTriSeg& s );
|
||||
|
||||
// return the master segment list
|
||||
inline triseg_list& get_seg_list() { return seg_list; }
|
||||
inline const triseg_list& get_seg_list() const { return seg_list; }
|
||||
|
||||
// return the ith segment
|
||||
inline TGTriSeg& get_seg( int i ) { return seg_list[i]; }
|
||||
inline const TGTriSeg& get_seg( int i ) const { return seg_list[i]; }
|
||||
};
|
||||
|
||||
|
||||
#endif // _TRISEGS_HXX
|
||||
|
||||
|
||||
|
|
@ -3,15 +3,25 @@ include_directories(${GDAL_INCLUDE_DIR})
|
|||
include( ${CGAL_USE_FILE} )
|
||||
|
||||
add_library(Polygon STATIC
|
||||
chop-bin.cxx
|
||||
index.cxx
|
||||
polygon.cxx
|
||||
simple_clip.cxx
|
||||
chop.hxx
|
||||
index.hxx
|
||||
names.hxx
|
||||
polygon.hxx
|
||||
simple_clip.hxx
|
||||
clipper.cpp
|
||||
clipper.hpp
|
||||
)
|
||||
polygon.cxx
|
||||
polygon.hxx
|
||||
rectangle.cxx
|
||||
rectangle.hxx
|
||||
tg_nodes.cxx
|
||||
tg_nodes.hxx
|
||||
tg_unique_geod.hxx
|
||||
tg_unique_tgnode.hxx
|
||||
tg_unique_vec2f.hxx
|
||||
tg_unique_vec3d.hxx
|
||||
tg_unique_vec3f.hxx
|
||||
trinodes.cxx
|
||||
trinodes.hxx
|
||||
TNT/jama_qr.h
|
||||
TNT/tnt_array1d.h
|
||||
TNT/tnt_array2d.h
|
||||
TNT/tnt_i_refvec.h
|
||||
TNT/tnt_math_utils.h
|
||||
point3d.hxx
|
||||
)
|
||||
File diff suppressed because it is too large
Load diff
|
|
@ -1,50 +0,0 @@
|
|||
// chop.hxx -- routine to chop a polygon up along tile boundaries and
|
||||
// write the individual pieces to the TG working polygon
|
||||
// file format.
|
||||
//
|
||||
// Written by Curtis Olson, started February 1999.
|
||||
//
|
||||
// Copyright (C) 1999-2004 Curtis L. Olson - http://www.flightgear.org/~curt
|
||||
//
|
||||
// This program is free software; you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation; either version 2 of the License, or
|
||||
// (at your option) any later version.
|
||||
//
|
||||
// This program is distributed in the hope that it will be useful,
|
||||
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
// GNU General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
|
||||
//
|
||||
// $Id: chop.hxx,v 1.3 2004-11-19 22:25:50 curt Exp $
|
||||
|
||||
|
||||
#ifndef _TG_CHOP_HXX
|
||||
#define _TG_CHOP_HXX
|
||||
|
||||
|
||||
#include <string>
|
||||
#include "polygon.hxx"
|
||||
|
||||
// process polygon shape (chop up along tile boundaries and write each
|
||||
// polygon piece to a file)
|
||||
void tgChopNormalPolygon( const std::string& path, const std::string& poly_type,
|
||||
const TGPolygon& shape, bool preserve3d );
|
||||
|
||||
void tgChopNormalPolygonsWithMask(const std::string& path, const std::string& poly_type,
|
||||
const poly_list& segments, bool preserve3d );
|
||||
|
||||
// process polygon shape (chop up along tile boundaries and write each
|
||||
// polygon piece to a file) This has a front end to a crude clipper
|
||||
// that doesn't handle holes so beware. This routine is appropriate
|
||||
// for breaking down really huge structures if needed.
|
||||
void tgChopBigSimplePolygon( const std::string& path, const std::string& poly_type,
|
||||
const TGPolygon& shape, bool preserve3d );
|
||||
|
||||
#endif // _TG_CHOP_HXX
|
||||
|
||||
|
||||
|
|
@ -1,76 +0,0 @@
|
|||
// index.cxx -- routines to handle a unique/persistant integer polygon index
|
||||
//
|
||||
// Written by Curtis Olson, started February 1999.
|
||||
//
|
||||
// Copyright (C) 1999 Curtis L. Olson - http://www.flightgear.org/~curt
|
||||
//
|
||||
// This program is free software; you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation; either version 2 of the License, or
|
||||
// (at your option) any later version.
|
||||
//
|
||||
// This program is distributed in the hope that it will be useful,
|
||||
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
// GNU General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
|
||||
//
|
||||
// $Id: index.cxx,v 1.7 2004-11-19 22:25:50 curt Exp $
|
||||
|
||||
#include <simgear/compiler.h>
|
||||
#include <simgear/debug/logstream.hxx>
|
||||
|
||||
#include <string>
|
||||
|
||||
#include <stdio.h>
|
||||
|
||||
#include "index.hxx"
|
||||
|
||||
using std::string;
|
||||
|
||||
static long int poly_index;
|
||||
static string poly_path;
|
||||
|
||||
// initialize the unique polygon index counter stored in path
|
||||
bool poly_index_init( string path )
|
||||
{
|
||||
poly_path = path;
|
||||
|
||||
FILE *fp = fopen( poly_path.c_str(), "r" );
|
||||
|
||||
if ( fp == NULL ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "Warning: cannot open " << poly_path);
|
||||
poly_index = 0;
|
||||
return false;
|
||||
}
|
||||
|
||||
fscanf( fp, "%ld", &poly_index );
|
||||
|
||||
fclose( fp );
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
// increment the persistant counter and return the next poly_index
|
||||
long int poly_index_next()
|
||||
{
|
||||
++poly_index;
|
||||
|
||||
FILE *fp = fopen( poly_path.c_str(), "w" );
|
||||
|
||||
if ( fp == NULL )
|
||||
SG_LOG(SG_GENERAL, SG_ALERT,
|
||||
"Error cannot open " << poly_path << " for writing");
|
||||
|
||||
fprintf( fp, "%ld\n", poly_index );
|
||||
|
||||
fclose( fp );
|
||||
|
||||
return poly_index;
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -1,42 +0,0 @@
|
|||
// index.cxx -- routines to handle a unique/persistant integer polygon index
|
||||
//
|
||||
// Written by Curtis Olson, started February 1999.
|
||||
//
|
||||
// Copyright (C) 1999 Curtis L. Olson - http://www.flightgear.org/~curt
|
||||
//
|
||||
// This program is free software; you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation; either version 2 of the License, or
|
||||
// (at your option) any later version.
|
||||
//
|
||||
// This program is distributed in the hope that it will be useful,
|
||||
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
// GNU General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
|
||||
//
|
||||
// $Id: index.hxx,v 1.4 2004-11-19 22:25:50 curt Exp $
|
||||
|
||||
|
||||
#ifndef _INDEX_HXX
|
||||
#define _INDEX_HXX
|
||||
|
||||
|
||||
#include <simgear/compiler.h>
|
||||
|
||||
#include <string>
|
||||
|
||||
// initialize the unique polygon index counter stored in path
|
||||
bool poly_index_init( std::string path );
|
||||
|
||||
// increment the persistant counter and return the next poly_index
|
||||
long int poly_index_next();
|
||||
|
||||
|
||||
|
||||
#endif // _INDEX_HXX
|
||||
|
||||
|
||||
|
|
@ -1 +0,0 @@
|
|||
|
||||
|
|
@ -1,48 +0,0 @@
|
|||
// names.hxx -- process shapefiles names
|
||||
//
|
||||
// Written by Curtis Olson, started February 1999.
|
||||
//
|
||||
// Copyright (C) 1999 Curtis L. Olson - http://www.flightgear.org/~curt
|
||||
//
|
||||
// This program is free software; you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation; either version 2 of the License, or
|
||||
// (at your option) any later version.
|
||||
//
|
||||
// This program is distributed in the hope that it will be useful,
|
||||
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
// GNU General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
|
||||
//
|
||||
// $Id: names.hxx,v 1.16 2007-10-31 15:05:13 curt Exp $
|
||||
|
||||
|
||||
#ifndef _NAMES_HXX
|
||||
#define _NAMES_HXX
|
||||
|
||||
|
||||
#include <simgear/compiler.h>
|
||||
|
||||
#include <string>
|
||||
|
||||
inline static bool is_ocean_area( const std::string &area )
|
||||
{
|
||||
return area == "Ocean" || area == "Bay Estuary or Ocean";
|
||||
}
|
||||
|
||||
inline static bool is_void_area( const std::string &area )
|
||||
{
|
||||
return area == "Void Area";
|
||||
}
|
||||
|
||||
inline static bool is_null_area( const std::string& area )
|
||||
{
|
||||
return area == "Null";
|
||||
}
|
||||
|
||||
#endif // _NAMES_HXX
|
||||
|
||||
|
|
@ -30,10 +30,9 @@
|
|||
#include <simgear/misc/texcoord.hxx>
|
||||
#include <simgear/structure/exception.hxx>
|
||||
#include <simgear/debug/logstream.hxx>
|
||||
#include <Geometry/point3d.hxx>
|
||||
#include <Geometry/poly_support.hxx>
|
||||
|
||||
#include <Geometry/trinodes.hxx>
|
||||
#include <Polygon/point3d.hxx>
|
||||
#include <Polygon/trinodes.hxx>
|
||||
|
||||
#include "polygon.hxx"
|
||||
|
||||
|
|
@ -76,7 +75,7 @@ void TGPolygon::get_bounding_box( SGGeod& min, SGGeod& max ) const
|
|||
max = SGGeod::fromDeg( maxx, maxy );
|
||||
}
|
||||
|
||||
|
||||
#if 0
|
||||
// Set the elevations of points in the current polgyon based on the
|
||||
// elevations of points in source. For points that are not found in
|
||||
// source, propogate the value from the nearest matching point.
|
||||
|
|
@ -142,7 +141,7 @@ void TGPolygon::inherit_elevations( const TGPolygon &source ) {
|
|||
}
|
||||
}
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
// Set the elevations of all points to the specified values
|
||||
void TGPolygon::set_elevations( double elev ) {
|
||||
|
|
@ -532,7 +531,7 @@ void clipper_to_shapefile( ClipperLib::Polygons polys, char* ds )
|
|||
TGPolygon tgcontour;
|
||||
char layer[32];
|
||||
|
||||
void* ds_id = tgShapefileOpenDatasource( ds );
|
||||
void* ds_id = tgShapefile::OpenDatasource( ds );
|
||||
|
||||
for (unsigned int i = 0; i < polys.size(); ++i) {
|
||||
if ( Orientation( polys[i] ) ) {
|
||||
|
|
@ -541,18 +540,18 @@ void clipper_to_shapefile( ClipperLib::Polygons polys, char* ds )
|
|||
sprintf( layer, "%04d_hole", i );
|
||||
}
|
||||
|
||||
void* l_id = tgShapefileOpenLayer( ds_id, layer );
|
||||
void* l_id = tgShapefile::OpenLayer( ds_id, layer );
|
||||
contour.clear();
|
||||
contour.push_back( polys[i] );
|
||||
|
||||
tgcontour.erase();
|
||||
make_tg_poly_from_clipper( contour, &tgcontour );
|
||||
|
||||
tgShapefileCreateFeature( ds_id, l_id, tgcontour, "contour" );
|
||||
tgShapefile::CreateFeature( ds_id, l_id, tgcontour, "contour" );
|
||||
}
|
||||
|
||||
// close after each write
|
||||
ds_id = tgShapefileCloseDatasource( ds_id );
|
||||
ds_id = tgShapefile::CloseDatasource( ds_id );
|
||||
}
|
||||
|
||||
TGPolygon polygon_clip_clipper( clip_op poly_op, const TGPolygon& subject, const TGPolygon& clip )
|
||||
|
|
@ -652,17 +651,17 @@ void tgPolygonFreeClipperAccumulator( void )
|
|||
|
||||
void tgPolygonDumpAccumulator( char* ds, char* layer, char* name )
|
||||
{
|
||||
void* ds_id = tgShapefileOpenDatasource( ds );
|
||||
void* l_id = tgShapefileOpenLayer( ds_id, layer );
|
||||
void* ds_id = tgShapefile::OpenDatasource( ds );
|
||||
void* l_id = tgShapefile::OpenLayer( ds_id, layer );
|
||||
TGPolygon accum;
|
||||
|
||||
for (unsigned int i=0; i<clipper_accumulator.size(); i++) {
|
||||
make_tg_poly_from_clipper( clipper_accumulator[i], &accum );
|
||||
tgShapefileCreateFeature( ds_id, l_id, accum, name );
|
||||
tgShapefile::CreateFeature( ds_id, l_id, accum, name );
|
||||
}
|
||||
|
||||
// close after each write
|
||||
ds_id = tgShapefileCloseDatasource( ds_id );
|
||||
ds_id = tgShapefile::CloseDatasource( ds_id );
|
||||
}
|
||||
|
||||
void tgPolygonAddToClipperAccumulator( const TGPolygon& subject, bool dump )
|
||||
|
|
@ -3053,11 +3052,11 @@ void tgPolygon::Texture( void )
|
|||
|
||||
void tgPolygon::ToShapefile( const tgPolygon& subject, const std::string& datasource, const std::string& layer, const std::string& description )
|
||||
{
|
||||
void* ds_id = tgShapefileOpenDatasource( datasource.c_str() );
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "tgShapefileOpenDatasource returned " << (unsigned long)ds_id);
|
||||
void* ds_id = tgShapefile::OpenDatasource( datasource.c_str() );
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "tgShapefile::OpenDatasource returned " << (unsigned long)ds_id);
|
||||
|
||||
OGRLayer* l_id = (OGRLayer *)tgShapefileOpenLayer( ds_id, layer.c_str() );
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "tgShapefileOpenLayer returned " << (unsigned long)l_id);
|
||||
OGRLayer* l_id = (OGRLayer *)tgShapefile::OpenLayer( ds_id, layer.c_str() );
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "tgShapefile::OpenLayer returned " << (unsigned long)l_id);
|
||||
|
||||
OGRPolygon* polygon = new OGRPolygon();
|
||||
|
||||
|
|
@ -3102,7 +3101,7 @@ void tgPolygon::ToShapefile( const tgPolygon& subject, const std::string& dataso
|
|||
}
|
||||
|
||||
// close after each write
|
||||
ds_id = tgShapefileCloseDatasource( ds_id );
|
||||
ds_id = tgShapefile::CloseDatasource( ds_id );
|
||||
}
|
||||
|
||||
tgPolygon tgPolygon::FromOGR( const OGRPolygon* subject )
|
||||
|
|
@ -3224,25 +3223,6 @@ bool tgPolygon::ChopIdxInit( const std::string& path )
|
|||
return true;
|
||||
}
|
||||
|
||||
SGMutex index_next_mutex;
|
||||
static long int tgPolygon_index_next()
|
||||
{
|
||||
SGGuard<SGMutex> g(index_next_mutex);
|
||||
|
||||
++poly_index;
|
||||
|
||||
FILE *fp = fopen( poly_path.c_str(), "w" );
|
||||
|
||||
if ( fp == NULL ) {
|
||||
SG_LOG(SG_GENERAL, SG_ALERT, "Error cannot open BLAH BLAG BLAH " << poly_path << " for writing");
|
||||
}
|
||||
|
||||
fprintf( fp, "%ld\n", poly_index );
|
||||
fclose( fp );
|
||||
|
||||
return poly_index;
|
||||
}
|
||||
|
||||
/************************ TESSELATION ***********************************/
|
||||
|
||||
void tg_mark_domains(CDT& ct, CDT::Face_handle start, int index, std::list<CDT::Edge>& border )
|
||||
|
|
@ -3993,3 +3973,121 @@ void tgChopper::Save( void )
|
|||
gzclose( fp );
|
||||
}
|
||||
}
|
||||
|
||||
const char* format_name="ESRI Shapefile";
|
||||
|
||||
void tgShapefile::Init( void )
|
||||
{
|
||||
OGRRegisterAll();
|
||||
}
|
||||
|
||||
void* tgShapefile::OpenDatasource( const char* datasource_name )
|
||||
{
|
||||
OGRDataSource *datasource;
|
||||
OGRSFDriver *ogrdriver;
|
||||
|
||||
ogrdriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(format_name);
|
||||
if (!ogrdriver) {
|
||||
SG_LOG(SG_GENERAL, SG_ALERT, "Unknown datasource format driver: " << format_name);
|
||||
exit(1);
|
||||
}
|
||||
|
||||
datasource = ogrdriver->Open(datasource_name, TRUE);
|
||||
|
||||
if (!datasource) {
|
||||
datasource = ogrdriver->CreateDataSource(datasource_name, NULL);
|
||||
}
|
||||
|
||||
if (!datasource) {
|
||||
SG_LOG(SG_GENERAL, SG_ALERT, "Unable to open or create datasource: " << datasource_name);
|
||||
exit(1);
|
||||
}
|
||||
|
||||
return (void*)datasource;
|
||||
}
|
||||
|
||||
void* tgShapefile::OpenLayer( void* ds_id, const char* layer_name ) {
|
||||
OGRDataSource* datasource = (OGRDataSource *)ds_id;
|
||||
OGRLayer* layer;
|
||||
|
||||
OGRSpatialReference srs;
|
||||
srs.SetWellKnownGeogCS("WGS84");
|
||||
|
||||
layer = datasource->GetLayerByName(layer_name);
|
||||
|
||||
if (!layer) {
|
||||
layer = datasource->CreateLayer( layer_name, &srs, wkbPolygon25D, NULL);
|
||||
|
||||
OGRFieldDefn descriptionField("ID", OFTString);
|
||||
descriptionField.SetWidth(128);
|
||||
|
||||
if( layer->CreateField( &descriptionField ) != OGRERR_NONE ) {
|
||||
SG_LOG(SG_GENERAL, SG_ALERT, "Creation of field 'Description' failed");
|
||||
}
|
||||
}
|
||||
|
||||
if (!layer) {
|
||||
SG_LOG(SG_GENERAL, SG_ALERT, "Creation of layer '" << layer_name << "' failed");
|
||||
return NULL;
|
||||
}
|
||||
|
||||
return (void*)layer;
|
||||
}
|
||||
|
||||
void tgShapefile::CreateFeature( void* ds_id, void* l_id, const TGPolygon &poly, const char* description )
|
||||
{
|
||||
OGRLayer* layer = (OGRLayer*)l_id;
|
||||
OGRPolygon* polygon = new OGRPolygon();
|
||||
|
||||
for ( int i = 0; i < poly.contours(); i++ ) {
|
||||
bool skip_ring=false;
|
||||
point_list contour = poly.get_contour( i );
|
||||
|
||||
if (contour.size()<3) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "Polygon with less than 3 points");
|
||||
skip_ring=true;
|
||||
}
|
||||
|
||||
// FIXME: Current we ignore the hole-flag and instead assume
|
||||
// that the first ring is not a hole and the rest
|
||||
// are holes
|
||||
OGRLinearRing *ring=new OGRLinearRing();
|
||||
for (unsigned int pt = 0; pt < contour.size(); pt++) {
|
||||
OGRPoint *point=new OGRPoint();
|
||||
|
||||
point->setX( contour[pt].x() );
|
||||
point->setY( contour[pt].y() );
|
||||
point->setZ( 0.0 );
|
||||
ring->addPoint(point);
|
||||
}
|
||||
ring->closeRings();
|
||||
|
||||
if (!skip_ring) {
|
||||
polygon->addRingDirectly(ring);
|
||||
}
|
||||
|
||||
OGRFeature* feature = NULL;
|
||||
feature = new OGRFeature( layer->GetLayerDefn() );
|
||||
feature->SetField("ID", description);
|
||||
feature->SetGeometry(polygon);
|
||||
if( layer->CreateFeature( feature ) != OGRERR_NONE )
|
||||
{
|
||||
SG_LOG(SG_GENERAL, SG_ALERT, "Failed to create feature in shapefile");
|
||||
}
|
||||
OGRFeature::DestroyFeature(feature);
|
||||
}
|
||||
}
|
||||
|
||||
void tgShapefile::CloseLayer( void* l_id )
|
||||
{
|
||||
//OGRLayer::DestroyLayer( layer );
|
||||
}
|
||||
|
||||
void* tgShapefile::CloseDatasource( void* ds_id )
|
||||
{
|
||||
OGRDataSource* datasource = (OGRDataSource *)ds_id;
|
||||
|
||||
OGRDataSource::DestroyDataSource( datasource );
|
||||
|
||||
return (void *)-1;
|
||||
}
|
||||
|
|
@ -29,16 +29,15 @@
|
|||
# error This library requires C++
|
||||
#endif
|
||||
|
||||
|
||||
#include <simgear/compiler.h>
|
||||
#include <simgear/math/sg_types.hxx>
|
||||
#include <Geometry/point3d.hxx>
|
||||
|
||||
|
||||
#include <iostream>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
#include <simgear/compiler.h>
|
||||
#include <simgear/math/sg_types.hxx>
|
||||
|
||||
#include <Polygon/point3d.hxx>
|
||||
|
||||
// forward declaration
|
||||
class TGPolygon;
|
||||
|
||||
|
|
@ -316,8 +315,8 @@ std::ostream &operator<<(std::ostream &output, const TGPolygon &poly);
|
|||
#include <simgear/bucket/newbucket.hxx>
|
||||
#include <simgear/threads/SGThread.hxx>
|
||||
|
||||
#include <Geometry/trinodes.hxx>
|
||||
#include <Geometry/rectangle.hxx>
|
||||
#include <Polygon/trinodes.hxx>
|
||||
#include <Polygon/rectangle.hxx>
|
||||
|
||||
#include "tg_unique_geod.hxx"
|
||||
|
||||
|
|
@ -486,6 +485,12 @@ public:
|
|||
return face_area;
|
||||
}
|
||||
|
||||
static double area( const SGGeod& p1, const SGGeod& p2, const SGGeod& p3 ) {
|
||||
return fabs(0.5 * ( p1.getLongitudeDeg() * p2.getLatitudeDeg() - p2.getLongitudeDeg() * p1.getLatitudeDeg() +
|
||||
p2.getLongitudeDeg() * p3.getLatitudeDeg() - p3.getLongitudeDeg() * p2.getLatitudeDeg() +
|
||||
p3.getLongitudeDeg() * p1.getLatitudeDeg() - p1.getLongitudeDeg() * p3.getLatitudeDeg() ));
|
||||
}
|
||||
|
||||
void SaveToGzFile( gzFile& fp ) const;
|
||||
void LoadFromGzFile( gzFile& fp );
|
||||
|
||||
|
|
@ -917,6 +922,17 @@ private:
|
|||
clipper_polygons_list accum;
|
||||
};
|
||||
|
||||
class tgShapefile
|
||||
{
|
||||
public:
|
||||
static void Init( void );
|
||||
static void* OpenDatasource( const char* datasource_name );
|
||||
static void* OpenLayer( void* ds_id, const char* layer_name );
|
||||
static void CreateFeature( void* ds_id, void* l_id, const TGPolygon &poly, const char* description );
|
||||
static void CloseLayer( void* l_id );
|
||||
static void* CloseDatasource( void* ds_id );
|
||||
};
|
||||
|
||||
|
||||
#endif // _POLYGON_HXX
|
||||
|
||||
|
|
|
|||
|
|
@ -1,550 +0,0 @@
|
|||
// simple_clip.cxx -- simplistic polygon clipping routine. Returns
|
||||
// the portion of a polygon that is above or below
|
||||
// a horizontal line of y = a. Only works with
|
||||
// single contour polygons (no holes.)
|
||||
//
|
||||
// Written by Curtis Olson, started August 1999.
|
||||
//
|
||||
// Copyright (C) 1999 Curtis L. Olson - http://www.flightgear.org/~curt
|
||||
//
|
||||
// This program is free software; you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation; either version 2 of the License, or
|
||||
// (at your option) any later version.
|
||||
//
|
||||
// This program is distributed in the hope that it will be useful,
|
||||
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
// GNU General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
|
||||
//
|
||||
// $Id: simple_clip.cxx,v 1.7 2004-11-19 22:25:50 curt Exp $
|
||||
|
||||
|
||||
#include <simgear/constants.h>
|
||||
#include <simgear/debug/logstream.hxx>
|
||||
#include <simgear/structure/exception.hxx>
|
||||
|
||||
#include "simple_clip.hxx"
|
||||
|
||||
#define CLIP_EPSILON 0.000000000001
|
||||
|
||||
|
||||
// Given a line segment specified by two endpoints p1 and p2, return
|
||||
// the x value of a point on the line that intersects with the
|
||||
// horizontal line through y. Return true if an intersection is found,
|
||||
// false otherwise.
|
||||
static bool intersects_y( Point3D p0, Point3D p1, double y, Point3D *result )
|
||||
{
|
||||
// sort the end points
|
||||
if ( p0.y() > p1.y() ) {
|
||||
Point3D tmp = p0;
|
||||
p0 = p1;
|
||||
p1 = tmp;
|
||||
}
|
||||
|
||||
if ( (y < p0.y()) || (y > p1.y()) )
|
||||
// out of range of line segment, bail right away
|
||||
return false;
|
||||
|
||||
// equation of a line through (x0,y0) and (x1,y1):
|
||||
//
|
||||
// y = y1 + (x - x1) * (y0 - y1) / (x0 - x1)
|
||||
// x = x1 + (y - y1) * (x0 - x1) / (y0 - y1)
|
||||
|
||||
double x;
|
||||
|
||||
if ( fabs(p0.y() - p1.y()) > CLIP_EPSILON )
|
||||
x = p1.x() + (y - p1.y()) * (p0.x() - p1.x()) / (p0.y() - p1.y());
|
||||
else
|
||||
return false;
|
||||
result->setx(x);
|
||||
result->sety(y);
|
||||
|
||||
if ( p0.x() <= p1.x() ) {
|
||||
if ( (p0.x() <= x) && (x <= p1.x()) )
|
||||
return true;
|
||||
} else if ( (p0.x() >= x) && (x >= p1.x()) )
|
||||
return true;
|
||||
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
// find the smallest point in contour 0 of poly such that x > min_x
|
||||
// and y = y. Returns index of the found point, -1 if no match found.
|
||||
static int find_point( const TGPolygon& poly, double min_x, double y )
|
||||
{
|
||||
Point3D p, save;
|
||||
int index = -1;
|
||||
|
||||
save.setx( 361.0 );
|
||||
|
||||
for ( int i = 0; i < poly.contour_size( 0 ); ++i ) {
|
||||
p = poly.get_pt( 0, i );
|
||||
if ( p.y() == y ) {
|
||||
// printf("(%d) p.y() = %.12f y = %.12f\n", i, p.y(), y);
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, " " << p);
|
||||
if ( p.x() > min_x ) {
|
||||
if ( p.x() < save.x() ) {
|
||||
save = p;
|
||||
index = i;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return index;
|
||||
}
|
||||
|
||||
|
||||
// return if interesection is valid (not in the ignore list)
|
||||
static bool valid_intersection( int intersection, const int_list& ignore_ints )
|
||||
{
|
||||
for ( int i = 0; i < (int)ignore_ints.size(); ++i )
|
||||
if ( intersection == ignore_ints[i] )
|
||||
return false;
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
// return index of next valid intersection
|
||||
static int next_intersection( const int_list& keep_ints,
|
||||
const int_list& ignore_ints,
|
||||
const int beginning_at )
|
||||
{
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, "[ni] start_int = " << beginning_at);
|
||||
int i = beginning_at;
|
||||
|
||||
if ( i < 0 ) i = 0; while ( i < (int)keep_ints.size() ) {
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, " i = " << i);
|
||||
if ( keep_ints[i] != -1 )
|
||||
if ( valid_intersection(keep_ints[i], ignore_ints) )
|
||||
return i;
|
||||
++i;
|
||||
}
|
||||
|
||||
return -1;
|
||||
}
|
||||
|
||||
|
||||
// return true if point.y() touches or is inside of line, else return false
|
||||
inline bool is_on_or_inside( double line, Point3D p, fgSideType side )
|
||||
{
|
||||
if ( side == Above ) {
|
||||
if ( p.y() >= line )
|
||||
return true;
|
||||
} else if ( side == Below )
|
||||
if ( p.y() <= line )
|
||||
return true;
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
// return true if point.y() is inside of line, else return false
|
||||
inline bool is_inside( double line, Point3D p, fgSideType side )
|
||||
{
|
||||
if ( side == Above ) {
|
||||
if ( p.y() > line )
|
||||
return true;
|
||||
} else if ( side == Below )
|
||||
if ( p.y() < line )
|
||||
return true;
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
// Walk through the input polygon and split it into the
|
||||
// portion that is inside the clip region
|
||||
static bool simple_clip( const TGPolygon& in, const double y,
|
||||
const fgSideType side,
|
||||
TGPolygon& result )
|
||||
{
|
||||
Point3D p, p_last, p_int;
|
||||
int i;
|
||||
|
||||
result.erase();
|
||||
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "input poly size = " << in.total_size());
|
||||
|
||||
p_last = in.get_pt( 0, in.contour_size(0) - 1 );
|
||||
|
||||
for ( i = 0; i < (int)in.contour_size(0); ++i ) {
|
||||
p = in.get_pt( 0, i );
|
||||
|
||||
if ( (fabs(p.x() - p_last.x()) < CLIP_EPSILON) &&
|
||||
(fabs(p.y() - p_last.y()) < CLIP_EPSILON) &&
|
||||
(i > 0) ) {
|
||||
// SG_LOG(SG_GENERAL, SG_WARN,
|
||||
// "WARNING: p and p_last are identical at index = " << i);
|
||||
}
|
||||
|
||||
if ( is_on_or_inside(y, p, side) ) {
|
||||
if ( is_on_or_inside(y, p_last, side) )
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, "inside & inside " << i << " " << p);
|
||||
result.add_node( 0, p );
|
||||
else {
|
||||
if ( !intersects_y(p, p_last, y, &p_int) ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "Huh, this should have intersected!");
|
||||
return false;
|
||||
} else {
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG,
|
||||
// "intersection outside to inside " << i << " "
|
||||
// << p_int " i - 1 = " << in.get_pt( 0, i-1 ));
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, " i = " << in.get_pt( 0, i ));
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, " i + 1 = " << in.get_pt( 0, i+1 ));
|
||||
result.add_node( 0, p_int );
|
||||
if ( (fabs(p.x() - p_int.x()) < CLIP_EPSILON) &&
|
||||
(fabs(p.y() - p_int.y()) < CLIP_EPSILON) ) {
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG,
|
||||
// "WARNING: p and p_int are identical, "
|
||||
// << "omitting p");
|
||||
} else {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "adding intersection" << i << " " << p);
|
||||
result.add_node( 0, p );
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
if ( is_inside(y, p_last, side) ) {
|
||||
if ( !intersects_y(p, p_last, y, &p_int) ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "Huh, this should have intersected!");
|
||||
return false;
|
||||
} else {
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG,
|
||||
// "intersection inside to outside " << i << " "
|
||||
// << p_int);
|
||||
if ( (fabs(p.x() - p_int.x()) < CLIP_EPSILON) &&
|
||||
(fabs(p.y() - p_int.y()) < CLIP_EPSILON) ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG,
|
||||
"WARNING: p and p_int are identical, "
|
||||
<< "omitting p");
|
||||
} else
|
||||
result.add_node( 0, p_int );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
p_last = p;
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
// build the list of intersections
|
||||
static bool build_intersections( const TGPolygon& arcs, double line,
|
||||
fgSideType side,
|
||||
int_list& keep_ints,
|
||||
int_list& ignore_ints )
|
||||
{
|
||||
keep_ints.clear();
|
||||
ignore_ints.clear();
|
||||
|
||||
int index = 0;
|
||||
double current_x = -181.0;
|
||||
|
||||
while ( index >= 0 ) {
|
||||
index = find_point( arcs, current_x, line );
|
||||
if ( index >= 0 ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "intersection at " << index << " = "
|
||||
<< arcs.get_pt( 0, index ));
|
||||
keep_ints.push_back( index );
|
||||
current_x = arcs.get_pt( 0, index ).x();
|
||||
|
||||
int before = index - 1;
|
||||
if ( before < 0 ) before += arcs.contour_size(0); int after = (index + 1) % arcs.contour_size(0);
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, " before = "
|
||||
<< arcs.get_pt(0, before));
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, " int = "
|
||||
<< arcs.get_pt(0, index));
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, " after = "
|
||||
<< arcs.get_pt(0, after));
|
||||
if ( side == Above ) {
|
||||
if ( (arcs.get_pt(0, before).y() > line) &&
|
||||
(arcs.get_pt(0, after).y() > line) ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "side = above");
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG,
|
||||
"V intersection with clip line from above");
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG,
|
||||
"Adding intersection to ignore_ints");
|
||||
ignore_ints.push_back( index );
|
||||
}
|
||||
if ( (arcs.get_pt(0, before).y() <= line) &&
|
||||
(arcs.get_pt(0, after).y() <= line) ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "side = above");
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG,
|
||||
"V intersection with clip line from BELOW\n"
|
||||
<< "or an extra in-clip-line intersection.\n"
|
||||
<< "Adding intersection to ignore_ints");
|
||||
ignore_ints.push_back( index );
|
||||
}
|
||||
} else if ( side == Below ) {
|
||||
if ( (arcs.get_pt(0, before).y() >= line) &&
|
||||
(arcs.get_pt(0, after).y() >= line) ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "side = below");
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "V intersection with clip line from above");
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "Adding intersection to ignore_ints");
|
||||
ignore_ints.push_back( index );
|
||||
}
|
||||
if ( (arcs.get_pt(0, before).y() < line) &&
|
||||
(arcs.get_pt(0, after).y() < line) ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "side = below");
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "V intersection with clip line from BELOW");
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "or an extra in-clip-line intersection.");
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "Adding intersection to ignore_ints");
|
||||
ignore_ints.push_back( index );
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
// test for duplicate nodes
|
||||
TGPolygon fix_dups( TGPolygon& in )
|
||||
{
|
||||
TGPolygon result;
|
||||
|
||||
double x_last = -20000.0;
|
||||
double y_last = -20000.0;
|
||||
double x, y;
|
||||
|
||||
for ( int i = 0; i < (int)in.contour_size(0); ++i ) {
|
||||
x = in.get_pt(0, i).x();
|
||||
y = in.get_pt(0, i).y();
|
||||
if ( (x == x_last) && (y == y_last) ) {
|
||||
// ignore
|
||||
} else
|
||||
result.add_node(0, in.get_pt(0, i));
|
||||
x_last = x;
|
||||
y_last = y;
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
|
||||
// simple polygon clipping routine. Returns the portion of a polygon
|
||||
// that is above and below a horizontal line of y = a. Only works
|
||||
// with single contour polygons (no holes.) Returns true if routine
|
||||
// thinks it was successful.
|
||||
|
||||
static bool clip_contour( const TGPolygon& in, const double y,
|
||||
const fgSideType side,
|
||||
TGPolygon& result )
|
||||
{
|
||||
TGPolygon result_arcs, arcs;
|
||||
int i, i1, i2, index;
|
||||
|
||||
|
||||
// Step 1: sanity checks
|
||||
if ( (int)in.contours() != 1 ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "we only handle single contour polygons");
|
||||
return false;
|
||||
}
|
||||
|
||||
if ( (int)in.contour_size( 0 ) < 3 ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "we must have at least three vertices to work");
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
// Step 2: walk through the input polygon and split it into the
|
||||
// portion that is on or inside the clip line
|
||||
|
||||
if ( simple_clip( in, y, side, result_arcs ) ) {
|
||||
if ( result_arcs.contours() > 0 )
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "result_arcs size = "
|
||||
<< result_arcs.total_size());
|
||||
else
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "empty result");
|
||||
} else
|
||||
throw sg_exception("simple_clip_above() failed!");
|
||||
|
||||
|
||||
// Step 3: check for trivial cases
|
||||
|
||||
result.erase();
|
||||
|
||||
// trivial -- nothing inside of clip line
|
||||
if ( result_arcs.contours() == 0 ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "trivially empty");
|
||||
return true;
|
||||
}
|
||||
|
||||
// trivial -- everything inside of clip line
|
||||
i1 = find_point( result_arcs, -181.0, y );
|
||||
if ( i1 < 0 ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "trivially full");
|
||||
result = result_arcs;
|
||||
return true;
|
||||
}
|
||||
|
||||
// trivial -- single clip line intersection (polygon just nicks
|
||||
// it) -- everything inside
|
||||
i2 = find_point( result_arcs, result_arcs.get_pt(0, i1).x(), y );
|
||||
if ( i2 < 0 ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "trivially full (clip line nicks edge)");
|
||||
result = result_arcs;
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
// Step 4: If we are finding the "below" clip, reverse the points
|
||||
// before extracting the cycles. (and remove duplicates)
|
||||
|
||||
TGPolygon tmp;
|
||||
tmp.erase();
|
||||
|
||||
if ( side == Below ) {
|
||||
for ( i = result_arcs.contour_size(0) - 1; i >= 0; --i ) {
|
||||
Point3D p = result_arcs.get_pt( 0, i );
|
||||
tmp.add_node( 0, p );
|
||||
}
|
||||
} else
|
||||
tmp = result_arcs;
|
||||
|
||||
arcs = fix_dups( tmp );
|
||||
|
||||
// Step 5: Build the intersections lists
|
||||
|
||||
int_list keep_ints, ignore_ints;
|
||||
build_intersections( arcs, y, side, keep_ints, ignore_ints );
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "total keep_ints = " << keep_ints.size());
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "total ignore_ints = " << ignore_ints.size());
|
||||
|
||||
|
||||
// Step 6: Walk through the result_arcs and extract the cycles (or
|
||||
// individual contours.)
|
||||
|
||||
int start_int = next_intersection( keep_ints, ignore_ints, 0 );
|
||||
int next_int = next_intersection( keep_ints, ignore_ints, start_int + 1 );
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "start_int = " << start_int);
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "next_int = " << next_int);
|
||||
|
||||
int count = 0;
|
||||
|
||||
// while we have keep_ints left to process
|
||||
while ( start_int >= 0 ) {
|
||||
point_list contour;
|
||||
contour.clear();
|
||||
|
||||
index = keep_ints[next_int];
|
||||
keep_ints[next_int] = -1;
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "\nstarting at point = "
|
||||
<< arcs.get_pt(0, index));
|
||||
|
||||
while ( index != keep_ints[start_int] ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "index = " << index
|
||||
<< " start_int = " << start_int
|
||||
<< " keep_ints[start_int] = " << keep_ints[start_int]);
|
||||
|
||||
// start with the 2nd item in the intersection list and
|
||||
// traverse until we find another intersection
|
||||
contour.push_back( arcs.get_pt(0, index) );
|
||||
index = (index + 1) % arcs.contour_size(0);
|
||||
|
||||
while ( (arcs.get_pt(0, index).y() != y) ||
|
||||
!valid_intersection(index, ignore_ints) ) {
|
||||
contour.push_back( arcs.get_pt(0, index) );
|
||||
index = (index + 1) % arcs.contour_size(0);
|
||||
}
|
||||
contour.push_back( arcs.get_pt(0, index) );
|
||||
SG_LOG(SG_GENERAL, SG_WARN, "exited at poly index = "
|
||||
<< index << " " << arcs.get_pt(0, index));
|
||||
|
||||
// find which intersection we came out on in our list
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "finding exit intersection for " << index);
|
||||
i = 0;
|
||||
while ( i < (int)keep_ints.size() ) {
|
||||
// SG_LOG(SG_GENERAL, SG_DEBUG, " keep_int[" << i << "] = " << keep_ints[i]);
|
||||
if ( index == keep_ints[i] ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, " intersection index = " << i);
|
||||
if ( index != keep_ints[start_int] ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, " not start index so keep going");
|
||||
keep_ints[i] = -1;
|
||||
next_int = next_intersection( keep_ints, ignore_ints,
|
||||
i + 1 );
|
||||
index = keep_ints[next_int];
|
||||
keep_ints[next_int] = -1;
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG,
|
||||
" next_int = " << next_int << " index = "
|
||||
<< index);
|
||||
}
|
||||
break;
|
||||
}
|
||||
++i;
|
||||
}
|
||||
if ( i == (int)keep_ints.size() )
|
||||
throw sg_exception("oops, didn't find that intersection, you are screwed");
|
||||
}
|
||||
keep_ints[start_int] = -1;
|
||||
result.add_contour( contour, count );
|
||||
count++;
|
||||
|
||||
// find next keep_ints
|
||||
start_int = next_intersection( keep_ints, ignore_ints, -1 );
|
||||
next_int = next_intersection( keep_ints, ignore_ints, start_int + 1 );
|
||||
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "start_int = " << start_int);
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "next_int = " << next_int);
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
// simple polygon clipping routine. Returns the portion of a polygon
|
||||
// that is above and below a horizontal line of y = a. Clips
|
||||
// multi-contour polygons individually and then reassembles the
|
||||
// results. Doesn't work with holes. Returns true if routine thinks
|
||||
// it was successful.
|
||||
|
||||
TGPolygon horizontal_clip( const TGPolygon& in, const double y,
|
||||
const fgSideType side )
|
||||
{
|
||||
TGPolygon result;
|
||||
|
||||
result.erase();
|
||||
|
||||
// Step 1: sanity checks
|
||||
if ( (int)in.contours() == 0 ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "Error: 0 contour polygon");
|
||||
return result;
|
||||
}
|
||||
|
||||
// clip each contour individually
|
||||
TGPolygon tmp_poly, clip_poly;
|
||||
point_list contour;
|
||||
|
||||
for ( int i = 0; i < in.contours(); ++i ) {
|
||||
if ( (int)in.contour_size( i ) < 3 ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "we must have at least three vertices to work");
|
||||
return result;
|
||||
}
|
||||
|
||||
tmp_poly.erase();
|
||||
|
||||
contour = in.get_contour( i );
|
||||
tmp_poly.add_contour( contour, 0 );
|
||||
|
||||
clip_contour( tmp_poly, y, side, clip_poly );
|
||||
|
||||
// add each clip_poly contour to the result poly
|
||||
for ( int j = 0; j < clip_poly.contours(); ++j ) {
|
||||
contour = clip_poly.get_contour( j );
|
||||
result.add_contour( contour, 0 );
|
||||
}
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
|
@ -1,53 +0,0 @@
|
|||
// simple_clip.hxx -- simplistic polygon clipping routine. Returns
|
||||
// the portion of a polygon that is above or below
|
||||
// a horizontal line of y = a. Only works with
|
||||
// single contour polygons (no holes.)
|
||||
//
|
||||
// Written by Curtis Olson, started August 1999.
|
||||
//
|
||||
// Copyright (C) 1999 Curtis L. Olson - http://www.flightgear.org/~curt
|
||||
//
|
||||
// This program is free software; you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation; either version 2 of the License, or
|
||||
// (at your option) any later version.
|
||||
//
|
||||
// This program is distributed in the hope that it will be useful,
|
||||
// but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
// GNU General Public License for more details.
|
||||
//
|
||||
// You should have received a copy of the GNU General Public License
|
||||
// along with this program; if not, write to the Free Software
|
||||
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
|
||||
//
|
||||
// $Id: simple_clip.hxx,v 1.4 2004-11-19 22:25:50 curt Exp $
|
||||
|
||||
|
||||
#ifndef _SIMPLE_CLIP_HXX
|
||||
#define _SIMPLE_CLIP_HXX
|
||||
|
||||
|
||||
#include "polygon.hxx"
|
||||
|
||||
|
||||
// side type
|
||||
enum fgSideType {
|
||||
Above = 0,
|
||||
Below = 1
|
||||
};
|
||||
|
||||
|
||||
// simple polygon clipping routine. Returns the portion of a polygon
|
||||
// that is above and below a horizontal line of y = a. Clips
|
||||
// multi-contour polygons individually and then reassembles the
|
||||
// results. Doesn't work with holes. Returns true if routine thinks
|
||||
// it was successful.
|
||||
|
||||
TGPolygon horizontal_clip( const TGPolygon& in, const double y,
|
||||
const fgSideType side );
|
||||
|
||||
|
||||
#endif // _SIMPLE_CLIP_HXX
|
||||
|
||||
|
||||
|
|
@ -31,9 +31,9 @@
|
|||
|
||||
|
||||
#include <simgear/compiler.h>
|
||||
#include <Geometry/point3d.hxx>
|
||||
#include <simgear/math/sg_types.hxx>
|
||||
|
||||
#include <Polygon/point3d.hxx>
|
||||
|
||||
#define FG_PROXIMITY_EPSILON 0.000001
|
||||
//#define FG_COURSE_EPSILON 0.0003
|
||||
|
|
@ -3,9 +3,9 @@ add_executable(ogr-decode ogr-decode.cxx)
|
|||
|
||||
target_link_libraries(ogr-decode
|
||||
${GDAL_LIBRARY}
|
||||
Polygon Geometry
|
||||
${SIMGEAR_CORE_LIBRARIES}
|
||||
${SIMGEAR_CORE_LIBRARY_DEPENDENCIES}
|
||||
Polygon
|
||||
${SIMGEAR_CORE_LIBRARIES}
|
||||
${SIMGEAR_CORE_LIBRARY_DEPENDENCIES}
|
||||
)
|
||||
|
||||
install(TARGETS ogr-decode RUNTIME DESTINATION bin)
|
||||
|
|
|
|||
|
|
@ -38,9 +38,6 @@
|
|||
#include <simgear/misc/sg_path.hxx>
|
||||
|
||||
#include <Include/version.h>
|
||||
#include <Polygon/chop.hxx>
|
||||
#include <Polygon/index.hxx>
|
||||
#include <Polygon/names.hxx>
|
||||
#include <Polygon/polygon.hxx>
|
||||
|
||||
/* stretch endpoints to reduce slivers in linear data ~.1 meters */
|
||||
|
|
@ -68,6 +65,19 @@ int num_threads = 1;
|
|||
|
||||
SGLockedQueue<OGRFeature *> global_workQueue;
|
||||
|
||||
/* very GDAL specific here... */
|
||||
inline static bool is_ocean_area( const std::string &area ) {
|
||||
return area == "Ocean" || area == "Bay Estuary or Ocean";
|
||||
}
|
||||
|
||||
inline static bool is_void_area( const std::string &area ) {
|
||||
return area == "Void Area";
|
||||
}
|
||||
|
||||
inline static bool is_null_area( const std::string& area ) {
|
||||
return area == "Null";
|
||||
}
|
||||
|
||||
class Decoder : public SGThread
|
||||
{
|
||||
public:
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue