fly/terragear
1
0
Fork 0
Code Activity

Initial revision.

Browse source
This commit is contained in:
curt 2001-04-24 21:09:49 +00:00
parent f6e553a1b6
commit 522bff9911
8 changed files with 997 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

11
src/BuildTiles/Osgb36/Makefile.am Normal file
View file

@ -0,0 +1,11 @@
noinst_LIBRARIES = libOsgb36.a
libOsgb36_a_SOURCES = osgb36.cxx osgb36.hxx osgbtc.cxx osgbtc.hxx uk.cxx uk.hxx
noinst_PROGRAMS = testosgb36
testosgb36_SOURCES = testosgb36.cxx
INCLUDES += \
-I$(top_srcdir) \
-I$(top_srcdir)/src/Lib

511
src/BuildTiles/Osgb36/osgb36.cxx Normal file
View file

@ -0,0 +1,511 @@
// osgb36.cxx -- Routines to convert UK OS coordinates to and
// from world WGS84 coordinates.
//
// Written by David Luff, started December 2000.
//
// Copyright (C) 2000 David C. Luff - david.luff@nottingham.ac.uk
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#include "osgb36.hxx"
double DEG_TO_RAD = 2.0 * 3.14159265358979323846264338327950288419716939967511 / 360.0;
double RAD_TO_DEG = 360.0 / (2.0 * 3.14159265358979323846264338327950288419716939967511);
//Convert OSGB36 Lat and Lon to OSGB36 Eastings and Northings (Grid Reference)
Point3D ConvertLatLonToEastingsNorthings(Point3D LatLon)
{
// ellipsoid constants
double a; //semi-major axis (m)
double b; //semi-minor axis (m)
double e_squared; //ellipsoid squared eccentricity constant
// Airy 1830
a = 6377563.396;
b = 6356256.910;
// GRS80 (aka WGS84)
// a = 6378137.000;
// b = 6356752.3141;
e_squared = ((a*a)-(b*b))/(a*a);
// Projection constants
double N_zero; //northing of true origin (m)
double E_zero; //easting of true origin (m)
double F_zero; //scale factor on central meridian
double lat_zero; //latitude of true origin (rad)
double lon_zero; //longitude of true origin and central meridian (rad)
// OSGB36 (UK National Grid)
N_zero = -100000;
E_zero = 400000;
F_zero = 0.9996012717;
lat_zero = 49 * DEG_TO_RAD;
lon_zero = -2 * DEG_TO_RAD;
double n;
double v;
double rho;
double neta_squared;
double M;
double I;
double II;
double III;
double IIIA;
double IV;
double V;
double VI;
double lon = LatLon.lon() * DEG_TO_RAD;
double lat = LatLon.lat() * DEG_TO_RAD;
n = (a-b)/(a+b);
v = a*F_zero*pow( (1-e_squared*sin(lat)*sin(lat)) , -0.5 );
rho = a*F_zero*(1-e_squared)*pow((1-e_squared*sin(lat)*sin(lat)) , -1.5 );
neta_squared = v/rho - 1.0;
// terms in the M equation to make it more manageable (eqn A11)
double term1;
double term2;
double term3;
double term4;
term1 = (1.0 + n + (5.0/4.0)*n*n + (5.0/4.0)*n*n*n)*(lat-lat_zero);
term2 = ((3.0*n) + (3.0*n*n) + ((21.0/8.0)*n*n*n))*(sin(lat-lat_zero))*(cos(lat+lat_zero));
term3 = (((15.0/8.0)*n*n) + ((15.0/8.0)*n*n*n))*(sin(2.0*(lat-lat_zero)))*(cos(2.0*(lat+lat_zero)));
term4 = ((35.0/24.0)*n*n*n)*(sin(3.0*(lat-lat_zero)))*(cos(3.0*(lat+lat_zero)));
M = b * F_zero * (term1 - term2 + term3 - term4); //Eqn A11
I = M + N_zero;
II = (v/2.0)*sin(lat)*cos(lat);
III = (v/24.0)*sin(lat)*pow(cos(lat),3)*(5.0 - pow(tan(lat),2) + 9.0*neta_squared);
IIIA = (v/720.0)*sin(lat)*pow(cos(lat),5)*(61.0 - 58.0*pow(tan(lat),2) + pow(tan(lat),4));
IV = v*cos(lat);
V = (v/6.0)*pow(cos(lat),3)*(v/rho - pow(tan(lat),2));
VI = (v/120.0)*pow(cos(lat),5)*(5.0 - 18.0*pow(tan(lat),2) + pow(tan(lat),4) + 14.0*neta_squared - 58.0*pow(tan(lat),2)*neta_squared);
double N; //calculated Northing
double E; //calculated Easting
N = I + II*(lon-lon_zero)*(lon-lon_zero) + III*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero) + IIIA*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero); //eqn A12
E = E_zero + IV*(lon-lon_zero) + V*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero) + VI*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero); //eqn A13
Point3D OSref;
OSref.setx(E);
OSref.sety(N);
OSref.setz(0.0); // Not used at present but could be used to return ODN height in future
return(OSref);
}
Point3D ConvertEastingsNorthingsToLatLon(Point3D GridRef)
{
// ellipsoid constants
double a; //semi-major axis (m)
double b; //semi-minor axis (m)
double e_squared; //ellipsoid squared eccentricity constant
// Airy 1830
a = 6377563.396;
b = 6356256.910;
// GRS80 (aka WGS84)
// a = 6378137.000;
// b = 6356752.3141;
e_squared = ((a*a)-(b*b))/(a*a);
// Projection constants
double N_zero; //northing of true origin (m)
double E_zero; //easting of true origin (m)
double F_zero; //scale factor on central meridian
double lat_zero; //latitude of true origin (rad)
double lon_zero; //longitude of true origin and central meridian (rad)
// OSGB36 (UK National Grid)
N_zero = -100000;
E_zero = 400000;
F_zero = 0.9996012717;
lat_zero = 49 * DEG_TO_RAD;
lon_zero = -2 * DEG_TO_RAD;
double n;
double M;
n = (a-b)/(a+b);
double E = GridRef.x();
double N = GridRef.y();
// terms in the M equation to make it more manageable (eqn A11)
double term1;
double term2;
double term3;
double term4;
double lat_prime;
//calculate initial lat_prime
lat_prime = ((N-N_zero)/(a*F_zero)) + lat_zero;
int i;
for(i=0;i<100;i++)
{
term1 = (1.0 + n + ((5.0/4.0)*n*n) + ((5.0/4.0)*n*n*n))*(lat_prime-lat_zero);
term2 = ((3.0*n) + (3.0*n*n) + ((21.0/8.0)*n*n*n))*(sin(lat_prime-lat_zero))*(cos(lat_prime+lat_zero));
term3 = (((15.0/8.0)*n*n) + ((15.0/8.0)*n*n*n))*(sin(2.0*(lat_prime-lat_zero)))*(cos(2.0*(lat_prime+lat_zero)));
term4 = ((35.0/24.0)*n*n*n)*(sin(3.0*(lat_prime-lat_zero)))*(cos(3.0*(lat_prime+lat_zero)));
M = b * F_zero * (term1 - term2 + term3 - term4); //Eqn A11
if((N - N_zero - M) < 0.001)
break; //N - N_zero - M is less than 1mm
else //recompute lat_prime and keep iterating until it is
lat_prime += ((N - N_zero - M)/(a*F_zero));
if(i == 99)
cout << "WARNING: Iteration failed to converge in ConvertEastingsNorthingsToLatLon\n";
}
double v;
double rho;
double neta_squared;
double VII;
double VIII;
double IX;
double X;
double XI;
double XII;
double XIIA;
v = a*F_zero*pow( (1-e_squared*sin(lat_prime)*sin(lat_prime)) , -0.5 );
rho = a*F_zero*(1-e_squared)*pow((1-e_squared*sin(lat_prime)*sin(lat_prime)) , -1.5 );
neta_squared = v/rho - 1.0;
VII = tan(lat_prime) / (2.0*rho*v);
VIII = (tan(lat_prime)/(24.0*rho*v*v*v)) * (5.0 + 3.0*pow(tan(lat_prime),2.0) + neta_squared - 9.0*pow(tan(lat_prime),2.0)*neta_squared);
IX = (tan(lat_prime)/(720.0*rho*v*v*v*v*v)) * (61.0 + 90.0*pow(tan(lat_prime),2.0) + 45.0*pow(tan(lat_prime),4.0));
X = (1/cos(lat_prime)) / v;
XI = ((1/cos(lat_prime)) / (6.0 * pow(v,3))) * (v/rho + 2*pow(tan(lat_prime),2));
XII = ((1/cos(lat_prime)) / (120.0 * pow(v,5))) * (5.0 + 28.0*pow(tan(lat_prime),2) + 24.0*pow(tan(lat_prime),4));
XIIA = ((1/cos(lat_prime)) / (5040.0 * pow(v,7))) * (61.0 + 662.0*pow(tan(lat_prime),2) + 1320.0*pow(tan(lat_prime),4) + 720.0*pow(tan(lat_prime),6));
double lat;
double lon;
lat = lat_prime - VII*pow((E-E_zero),2) + VIII*pow((E-E_zero),4) + IX*pow((E-E_zero),6);
lon = lon_zero + X*(E-E_zero) - XI*pow((E-E_zero),3) + XII*pow((E-E_zero),5) - XIIA*pow((E-E_zero),7);
Point3D LatLon;
LatLon.setlon(lon * RAD_TO_DEG);
LatLon.setlat(lat * RAD_TO_DEG);
LatLon.setelev(0.0);
return(LatLon);
}
Point3D ConvertAiry1830PolarToCartesian(Point3D LatLon)
{
double lon = LatLon.lon() * DEG_TO_RAD;
double lat = LatLon.lat() * DEG_TO_RAD;
double H = LatLon.elev() * DEG_TO_RAD;
// ellipsoid constants
double a; //semi-major axis (m)
double b; //semi-minor axis (m)
double e_squared; //ellipsoid squared eccentricity constant
// Airy 1830
a = 6377563.396;
b = 6356256.910;
e_squared = ((a*a)-(b*b))/(a*a);
double v;
v = a*pow( (1-e_squared*sin(lat)*sin(lat)) , -0.5 );
double x;
double y;
double z;
Point3D Cartesian;
x = (v+H)*cos(lat)*cos(lon);
y = (v+H)*cos(lat)*sin(lon);
z = ((1-e_squared)*v + H)*sin(lat);
Cartesian.setx(x);
Cartesian.sety(y);
Cartesian.setz(z);
return(Cartesian);
}
Point3D ConvertGRS80PolarToCartesian(Point3D LatLon)
{
double lon = LatLon.lon() * DEG_TO_RAD;
double lat = LatLon.lat() * DEG_TO_RAD;
double H = LatLon.elev() * DEG_TO_RAD;
// ellipsoid constants
double a; //semi-major axis (m)
double b; //semi-minor axis (m)
double e_squared; //ellipsoid squared eccentricity constant
// GRS80 (aka WGS84)
a = 6378137.000;
b = 6356752.3141;
e_squared = ((a*a)-(b*b))/(a*a);
double v;
v = a*pow( (1-e_squared*sin(lat)*sin(lat)) , -0.5 );
double x;
double y;
double z;
Point3D Cartesian;
x = (v+H)*cos(lat)*cos(lon);
y = (v+H)*cos(lat)*sin(lon);
z = ((1-e_squared)*v + H)*sin(lat);
Cartesian.setx(x);
Cartesian.sety(y);
Cartesian.setz(z);
return(Cartesian);
}
Point3D ConvertAiry1830CartesianToPolar(Point3D XYZCoord)
{
double x = XYZCoord.x();
double y = XYZCoord.y();
double z = XYZCoord.z();
double p;
double lat1;
double lat2;
double lat;
double lon;
double H;
Point3D LatLon;
// ellipsoid constants
double a; //semi-major axis (m)
double b; //semi-minor axis (m)
double e_squared; //ellipsoid squared eccentricity constant
// Airy 1830
a = 6377563.396;
b = 6356256.910;
e_squared = ((a*a)-(b*b))/(a*a);
double v;
lon = atan(y/x);
p = sqrt(x*x + y*y);
//lat is obtained iteratively
//obtain initial value of lat
lat1 = atan(z/(p*(1-e_squared)));
// cout << "Initial value of lat = " << lat1 << '\n';
int i;
for(i=0;i<100;i++)
{
v = a*pow( (1-e_squared*sin(lat1)*sin(lat1)) , -0.5 );
lat2 = atan((z + e_squared*v*sin(lat1))/p);
// cout << "lat2 = " << lat2 << '\n';
if(fabs(lat2-lat1) < 0.00000001)
{
lat = lat2;
// cout << i << " iterations required in ConvertAiry1830CartesianToPolar\n";
break;
}
else
{
lat1 = lat2;
}
if(i == 99)
cout << "WARNING: Iteration failed to converge in ConvertAiry1830CartesianToPolar\n";
}
H = p/cos(lat) - v;
LatLon.setlon(lon * RAD_TO_DEG);
LatLon.setlat(lat * RAD_TO_DEG);
LatLon.setelev(H);
return(LatLon);
}
Point3D ConvertGRS80CartesianToPolar(Point3D XYZCoord)
{
double x = XYZCoord.x();
double y = XYZCoord.y();
double z = XYZCoord.z();
double p;
double lat1;
double lat2;
double lat;
double lon;
double H;
Point3D LatLon;
// ellipsoid constants
double a; //semi-major axis (m)
double b; //semi-minor axis (m)
double e_squared; //ellipsoid squared eccentricity constant
// GRS80 (aka WGS84)
a = 6378137.000;
b = 6356752.3141;
e_squared = ((a*a)-(b*b))/(a*a);
double v;
lon = atan(y/x);
p = sqrt(x*x + y*y);
//lat is obtained iteratively
//obtain initial value of lat
lat1 = atan(z/(p*(1-e_squared)));
int i;
for(i=0;i<100;i++)
{
v = a*pow( (1-e_squared*sin(lat1)*sin(lat1)) , -0.5 );
lat2 = atan((z + e_squared*v*sin(lat1))/p);
if(fabs(lat2-lat1) < 0.00000001)
{
lat = lat2;
// cout << i << " iterations required in ConvertGRS80CartesianToPolar\n";
break;
}
else
{
lat1 = lat2;
}
if(i == 99)
cout << "WARNING: Iteration failed to converge in ConvertGRS80CartesianToPolar\n";
}
H = p/cos(lat) - v;
LatLon.setlon(lon * RAD_TO_DEG);
LatLon.setlat(lat * RAD_TO_DEG);
LatLon.setelev(H);
return(LatLon);
}
//Transform a point in WGS84 cartesian coordinates to OSGB36 cartesian coordinates
//Uses the Helmert Transformation with coefficients from www.gps.gov.org
//Only accurate to around 5m since OSGB36 is an imhomogenous TRF
Point3D ConvertWGS84ToOSGB36(Point3D WGS84)
{
Point3D OSGB36;
double tx = -446.448; // (m)
double ty = 125.157;
double tz = -542.060;
double rx = -0.1502; // (sec)
double ry = -0.2470;
double rz = -0.8421;
double s = 20.4894; // scale factor in ppm
// Convert rotations from seconds of arc to radians
rx = rx / 3600.0 * DEG_TO_RAD;
ry = ry / 3600.0 * DEG_TO_RAD;
rz = rz / 3600.0 * DEG_TO_RAD;
s /= 1000000;
// cout << "(1+s) = " << (1+s) << '\n';
OSGB36.setx(tx + (1+s)*WGS84.x() - rz*WGS84.y() + ry*WGS84.z());
OSGB36.sety(ty + rz*WGS84.x() + (1+s)*WGS84.y() - rx*WGS84.z());
OSGB36.setz(tz - ry*WGS84.x() + rx*WGS84.y() + (1+s)*WGS84.z());
return(OSGB36);
}
//Reverse transformation achieved by reversing the signs of coefficients for the above.
//Only valid for small angles, and has same accuracy limitations
Point3D ConvertOSGB36ToWGS84(Point3D OSGB36)
{
Point3D WGS84;
double tx = 446.448; // (m)
double ty = -125.157;
double tz = 542.060;
double rx = 0.1502; // (sec)
double ry = 0.2470;
double rz = 0.8421;
double s = -20.4894; // scale factor in ppm
// Convert rotations from seconds to radians
rx = rx / 3600.0 * DEG_TO_RAD;
ry = ry / 3600.0 * DEG_TO_RAD;
rz = rz / 3600.0 * DEG_TO_RAD;
s /= 1000000;
WGS84.setx(tx + (1+s)*OSGB36.x() - rz*OSGB36.y() + ry*OSGB36.z());
WGS84.sety(ty + rz*OSGB36.x() + (1+s)*OSGB36.y() - rx*OSGB36.z());
WGS84.setz(tz - ry*OSGB36.x() + rx*OSGB36.y() + (1+s)*OSGB36.z());
return(WGS84);
}

64
src/BuildTiles/Osgb36/osgb36.hxx Normal file
View file

@ -0,0 +1,64 @@
// osgb36.hxx -- Routines to convert UK OS coordinates to and
// from world WGS84 coordinates.
//
// Written by David Luff, started December 2000.
//
// Copyright (C) 2000 David C. Luff - david.luff@nottingham.ac.uk
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#include <iostream.h>
#include <math.h>
#include <simgear/math/point3d.hxx>
//*******************************************************************
//
// NOTE All functions which take or return Latitude and Longitude
// take and return degrees. Conversion to and from radians
// is performed internaly. XYZ coordinates and ellipsoid
// heights are returned and required in meters.
//
//*******************************************************************
//Convert OSGB36 Lat and Lon to OSGB36 Eastings and Northings (Grid Reference)
Point3D ConvertLatLonToEastingsNorthings(Point3D LatLon);
//Convert OSGB36 Eastings and Northings to OSGB36 Lat and Lon
Point3D ConvertEastingsNorthingsToLatLon(Point3D GridRef);
//Convert OSGB36 coordinates from polar to cartesian (x,y,z) form
Point3D ConvertAiry1830PolarToCartesian(Point3D LatLon);
//Convert WGS84 coordinates from polar to cartesian (x,y,z) form
Point3D ConvertGRS80PolarToCartesian(Point3D LatLon);
//Convert OSGB36 coordinates from cartesian (x,y,z) to polar form
Point3D ConvertAiry1830CartesianToPolar(Point3D XYZCoord);
//Convert WGS84 coordinates from cartesian (x,y,z) to polar form
Point3D ConvertGRS80CartesianToPolar(Point3D XYZCoord);
//Transform a point in WGS84 cartesian coordinates to OSGB36 cartesian coordinates
//Uses the Helmert Transformation with coefficients from www.gps.gov.org
//Only accurate to around 5m since OSGB36 is an inhomogenous TRF
Point3D ConvertWGS84ToOSGB36(Point3D WGS84);
//Transform a point in OSGB36 cartesian coordinates to WGS84 cartesian coordinates
//Uses the Helmert Transformation with coefficients from www.gps.gov.org
//Only accurate to around 5m since OSGB36 is an inhomogenous TRF
Point3D ConvertOSGB36ToWGS84(Point3D OSGB36);

220
src/BuildTiles/Osgb36/osgbtc.cxx Normal file
View file

@ -0,0 +1,220 @@
#include "osgb36.hxx"
#include "uk.hxx"
#include "osgbtc.hxx"
// traverse the specified fan/strip/list of vertices and attempt to
// calculate "none stretching" texture coordinates
point_list UK_calc_tex_coords( const SGBucket& b, const point_list& geod_nodes,
const int_list& fan, double scale )
{
// cout << "calculating texture coordinates for a specific fan of size = "
// << fan.size() << endl;
// calculate perimeter based on center of this degree (not center
// of bucket)
/* double clat = (int)b.get_center_lat();
if ( clat > 0 ) {
clat = (int)clat + 0.5;
} else {
clat = (int)clat - 0.5;
}
*/
// double clat_rad = clat * DEG_TO_RAD;
// double cos_lat = cos( clat_rad );
// double local_radius = cos_lat * EQUATORIAL_RADIUS_M;
// double local_perimeter = 2.0 * local_radius * FG_PI;
// double degree_width = local_perimeter / 360.0;
// cout << "clat = " << clat << endl;
// cout << "clat (radians) = " << clat_rad << endl;
// cout << "cos(lat) = " << cos_lat << endl;
// cout << "local_radius = " << local_radius << endl;
// cout << "local_perimeter = " << local_perimeter << endl;
// cout << "degree_width = " << degree_width << endl;
// double perimeter = 2.0 * EQUATORIAL_RADIUS_M * FG_PI;
// double degree_height = perimeter / 360.0;
// cout << "degree_height = " << degree_height << endl;
// find min/max of fan
Point3D tmin, tmax, WGS84LatLon, t;
Point3D WGS84xyz, OSGB36xyz, OSGB36LatLon;
Point3D OSGridRef;
bool first = true;
int i;
for ( i = 0; i < (int)fan.size(); ++i ) {
WGS84LatLon = geod_nodes[ fan[i] ];
// cout << "point WGS84LatLon = " << WGS84LatLon << endl;
WGS84xyz = ConvertGRS80PolarToCartesian(WGS84LatLon);
// cout << "WGS84XYZ = " << WGS84XYZ.x << ", " << WGS84XYZ.y << ", " << WGS84XYZ.z << '\n';
OSGB36xyz = ConvertWGS84ToOSGB36(WGS84xyz);
// cout << "OSXYZ = " << OSXYZ.x << ", " << OSXYZ.y << ", " << OSXYZ.z << '\n';
OSGB36LatLon = ConvertAiry1830CartesianToPolar(OSGB36xyz);
// cout << "OSGB36LatLon = " << OSGB36LatLon.x() << ", " << OSGB36LatLon.y() << ", " << OSGB36LatLon.z() << '\n';
OSGridRef = ConvertLatLonToEastingsNorthings(OSGB36LatLon);
// cout << "OS Eastings and Northings = " << OSGridRef << '\n';
t = UK_basic_tex_coord( OSGridRef );
// cout << "basic_tex_coord = " << t << endl;
// cout << "p = " << p << '\n';
// cout << "t = " << t << '\n';
// cout << "Degree width = " << degree_width << '\n';
// cout << "Degree height = " << degree_height << '\n';
// char dcl_pause;
// cin >> dcl_pause;
if ( first ) {
tmin = tmax = t;
first = false;
} else {
if ( t.x() < tmin.x() ) {
tmin.setx( t.x() );
}
if ( t.y() < tmin.y() ) {
tmin.sety( t.y() );
}
if ( t.x() > tmax.x() ) {
tmax.setx( t.x() );
}
if ( t.y() > tmax.y() ) {
tmax.sety( t.y() );
}
}
// cout << "tmin.x = " << tmin.x() << " tmax.x = " << tmax.x() << " tmin.y = " << tmin.y() << " tmax.y = " << tmax.y() << '\n';
}
double dx = fabs( tmax.x() - tmin.x() );
double dy = fabs( tmax.y() - tmin.y() );
// cout << "dx = " << dx << " dy = " << dy << endl;
bool do_shift = false;
// Point3D mod_shift;
if ( (dx > HALF_MAX_TEX_COORD) || (dy > HALF_MAX_TEX_COORD) ) {
// structure is too big, we'll just have to shift it so that
// tmin = (0,0). This messes up subsequent texture scaling,
// but is the best we can do.
// cout << "SHIFTING" << endl;
do_shift = true;
if ( tmin.x() < 0 ) {
tmin.setx( (double)( (int)tmin.x() - 1 ) );
} else {
tmin.setx( (int)tmin.x() );
}
if ( tmin.y() < 0 ) {
tmin.sety( (double)( (int)tmin.y() - 1 ) );
} else {
tmin.sety( (int)tmin.y() );
}
// cout << "found tmin = " << tmin << endl;
} else {
if ( tmin.x() < 0 ) {
tmin.setx( ( (int)(tmin.x() / HALF_MAX_TEX_COORD) - 1 )
* HALF_MAX_TEX_COORD );
} else {
tmin.setx( ( (int)(tmin.x() / HALF_MAX_TEX_COORD) )
* HALF_MAX_TEX_COORD );
}
if ( tmin.y() < 0 ) {
tmin.sety( ( (int)(tmin.y() / HALF_MAX_TEX_COORD) - 1 )
* HALF_MAX_TEX_COORD );
} else {
tmin.sety( ( (int)(tmin.y() / HALF_MAX_TEX_COORD) )
* HALF_MAX_TEX_COORD );
}
#if 0
// structure is small enough ... we can mod it so we can
// properly scale the texture coordinates later.
// cout << "MODDING" << endl;
double x1 = fmod(tmin.x(), MAX_TEX_COORD);
while ( x1 < 0 ) { x1 += MAX_TEX_COORD; }
double y1 = fmod(tmin.y(), MAX_TEX_COORD);
while ( y1 < 0 ) { y1 += MAX_TEX_COORD; }
double x2 = fmod(tmax.x(), MAX_TEX_COORD);
while ( x2 < 0 ) { x2 += MAX_TEX_COORD; }
double y2 = fmod(tmax.y(), MAX_TEX_COORD);
while ( y2 < 0 ) { y2 += MAX_TEX_COORD; }
// At this point we know that the object is < 16 wide in
// texture coordinate space. If the modulo of the tmin is >
// the mod of the tmax at this point, then we know that the
// starting tex coordinate for the tmax > 16 so we can shift
// everything down by 16 and get it within the 0-32 range.
if ( x1 > x2 ) {
mod_shift.setx( HALF_MAX_TEX_COORD );
} else {
mod_shift.setx( 0.0 );
}
if ( y1 > y2 ) {
mod_shift.sety( HALF_MAX_TEX_COORD );
} else {
mod_shift.sety( 0.0 );
}
#endif
// cout << "mod_shift = " << mod_shift << endl;
}
// generate tex_list
Point3D adjusted_t;
point_list tex;
tex.clear();
for ( i = 0; i < (int)fan.size(); ++i ) {
WGS84LatLon = geod_nodes[ fan[i] ];
// cout << "About to do second run through tex coords\n";
// cout << "point WGS84LatLon = " << WGS84LatLon << endl;
WGS84xyz = ConvertGRS80PolarToCartesian(WGS84LatLon);
// cout << "WGS84XYZ = " << WGS84XYZ.x << ", " << WGS84XYZ.y << ", " << WGS84XYZ.z << '\n';
OSGB36xyz = ConvertWGS84ToOSGB36(WGS84xyz);
// cout << "OSXYZ = " << OSXYZ.x << ", " << OSXYZ.y << ", " << OSXYZ.z << '\n';
OSGB36LatLon = ConvertAiry1830CartesianToPolar(OSGB36xyz);
// cout << "OSGB36LatLon = " << OSGB36LatLon.x() << ", " << OSGB36LatLon.y() << ", " << OSGB36LatLon.z() << '\n';
OSGridRef = ConvertLatLonToEastingsNorthings(OSGB36LatLon);
// cout << "OS Eastings and Northings = " << OSGridRef << '\n';
t = UK_basic_tex_coord(OSGridRef);
// cout << "second t = " << t << endl;
// if ( do_shift ) {
adjusted_t = t - tmin;
// cout << "adjusted_t = " << adjusted_t << '\n';
#if 0
} else {
adjusted_t.setx( fmod(t.x() + mod_shift.x(), MAX_TEX_COORD) );
while ( adjusted_t.x() < 0 ) {
adjusted_t.setx( adjusted_t.x() + MAX_TEX_COORD );
}
adjusted_t.sety( fmod(t.y() + mod_shift.y(), MAX_TEX_COORD) );
while ( adjusted_t.y() < 0 ) {
adjusted_t.sety( adjusted_t.y() + MAX_TEX_COORD );
}
// cout << "adjusted_t " << adjusted_t << endl;
}
#endif
if ( adjusted_t.x() < SG_EPSILON ) {
adjusted_t.setx( 0.0 );
}
if ( adjusted_t.y() < SG_EPSILON ) {
adjusted_t.sety( 0.0 );
}
adjusted_t.setz( 0.0 );
cout << "adjusted_t after check for SG_EPSILON = " << adjusted_t << endl;
tex.push_back( adjusted_t );
}
// char dcl_pause2;
// cin >> dcl_pause2;
return tex;
}

31
src/BuildTiles/Osgb36/osgbtc.hxx Normal file
View file

@ -0,0 +1,31 @@
#ifndef OSBTEXCOORD_H
#define OSBTEXCOORD_H
#include <simgear/bucket/newbucket.hxx>
#include <simgear/math/point3d.hxx>
#include <simgear/math/sg_types.hxx>
#define FG_STANDARD_TEXTURE_DIMENSION 1000.0 // meters
#define MAX_TEX_COORD 8.0
#define HALF_MAX_TEX_COORD ( MAX_TEX_COORD / 2.0 )
// return the basic unshifted/unmoded texture coordinate for a Easting/Northing (UK) grid reference
inline Point3D UK_basic_tex_coord( const Point3D& p )
{
// cout << "Point in dcl_basic_tex_coord is " << p << '\n';
return Point3D( p.x() / FG_STANDARD_TEXTURE_DIMENSION ,
p.y() / FG_STANDARD_TEXTURE_DIMENSION ,
0.0 );
}
// traverse the specified fan/strip/list of vertices and attempt to
// calculate "none stretching" texture coordinates
point_list UK_calc_tex_coords( const SGBucket& b, const point_list& geod_nodes,
const int_list& fan, double scale );
#endif // OSBTEXCOORD_H

74
src/BuildTiles/Osgb36/testosgb36.cxx Normal file
View file

@ -0,0 +1,74 @@
#include <simgear/math/point3d.hxx>
#include "osgb36.hxx"
//Test harness for various WGS84 and OSGB36 conversion functions
int main()
{
cout << "Test OSGB36\n";
/* Point3D initial_position;
Point3D final_position;
initial_position.setlat(52.0 + ((39.0 + 27.2531/60.0)/60.0));
initial_position.setlon(1.0 + ((43.0 + 4.5177/60.0)/60.0));
initial_position.setelev(0.0);
final_position = ConvertLatLonToEastingsNorthings(initial_position);
cout << "Eastings = " << final_position.x() << '\n';
cout << "Northings = " << final_position.y() << '\n';
Point3D check_position;
check_position = ConvertEastingsNorthingsToLatLon(final_position);
cout << "Lat = " << check_position.lat() << '\n';
cout << "Lon = " << check_position.lon() << '\n';
cout << "Elipsiodal height = " << check_position.elev() << '\n';
check_position.setelev(24.7);
Point3D cartesian_position = ConvertAiry1830PolarToCartesian(check_position);
cout << "X = " << cartesian_position.x() << '\n';
cout << "Y = " << cartesian_position.y() << '\n';
cout << "Z = " << cartesian_position.z() << '\n';
cout << "\nTesting ConvertAiry1830CartesianToPolar.....\n";
Point3D XYZPosition;
XYZPosition.setx(3874938.849);
XYZPosition.sety(116218.624);
XYZPosition.setz(5047168.208);
Point3D LatLonPosition = ConvertAiry1830CartesianToPolar(XYZPosition);
cout << "Lat = " << LatLonPosition.lat() << '\n';
cout << "Lon = " << LatLonPosition.lon() << '\n';
cout << "Elipsiodal height = " << LatLonPosition.elev() << '\n';
cout << "\nNow attempting to convert WGS84 lat & lon to OS grid reference .....\n";
Point3D WGS84LatLon;
Point3D WGS84XYZ;
Point3D OSXYZ;
Point3D OSLatLon;
Point3D OSGridRef;
WGS84LatLon.setlat(52.0);
WGS84LatLon.setlon(-2.0);
WGS84LatLon.setelev(0.0);
WGS84XYZ = ConvertGRS80PolarToCartesian(WGS84LatLon);
cout << "WGS84XYZ = " << WGS84XYZ.x() << ", " << WGS84XYZ.y() << ", " << WGS84XYZ.z() << '\n';
OSXYZ = ConvertWGS84ToOSGB36(WGS84XYZ);
cout << "OSXYZ = " << OSXYZ.x() << ", " << OSXYZ.y() << ", " << OSXYZ.z() << '\n';
OSLatLon = ConvertAiry1830CartesianToPolar(OSXYZ);
cout << "OSLatLon = " << OSLatLon.lat() << ", " << OSLatLon.lon() << ", " << OSLatLon.elev() << '\n';
OSGridRef = ConvertLatLonToEastingsNorthings(OSLatLon);
cout << "OS Grid Reference = " << OSGridRef.x() << ", " << OSGridRef.y() << "\n\n";
*/
return(0);
}

61
src/BuildTiles/Osgb36/uk.cxx Normal file
View file

@ -0,0 +1,61 @@
// uk.cxx -- Routines to determine whether a point is within the UK.
//
// Written by David Luff, started Janurary 2001.
//
// Copyright (C) 2001 David C. Luff - david.luff@nottingham.ac.uk
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#include "uk.hxx"
//Returns true if a point is within the mainland UK, excluding
//Northern Ireland. Requires lat and lon to be passed in degrees.
bool isInUK(Point3D p)
{
bool inUK = false;
double lat = p.lat();
double lon = p.lon();
//The entire UK (excepting a small portion of Northern Ireland which we arn't doing anyway)
//falls within the box of -8 -> 2 deg lon, and 49 -> 60 deg lat
//We'll check for within this box first, and if so then we'll refine
//further to avoid bits of France and Ireland.
if((lat >= 49.0) && (lat <= 60.0) && (lon >= -8.0) && (lon <= 2.0))
{
//we might be within the UK
//set inUK = true and then test for the exceptions
inUK = true;
//check for France (Normandy/Calais)
if((lon >= -2.0) && (lat <=50.0))
//Normandy
inUK = false;
if((lon > 1.0) && (lat <= 51.0))
//Calais area
inUK = false;
//Check for Ireland
if((lon <= -6.0) && (lat > 51.0) && (lat <= 55.5))
//Ireland
inUK = false;
//Check for last bit of NI
if((lat > 54.0) && (lat <= 55.1) && (lon <= -5.3333))
inUK = false;
}
else
inUK = false;
return(inUK);
}

25
src/BuildTiles/Osgb36/uk.hxx Normal file
View file

@ -0,0 +1,25 @@
// uk.hxx -- Routines to determine whether a point is within the UK.
//
// Written by David Luff, started Janurary 2001.
//
// Copyright (C) 2001 David C. Luff - david.luff@nottingham.ac.uk
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#include <simgear/math/point3d.hxx>
//Returns true if a point is within the mainland UK, excluding
//Northern Ireland. Requires lat and lon to be passed in degrees.
bool isInUK(Point3D p);