/*
This projection contributed by Ian Turner <vectro@vectro.org>. His
comments:
I have created a new projection for XPlanet, entitled
"icosagnomonic". It is very similar to the Fuller projection, but I
have used the gnomonic projection instead of fuller's method to
perform the transformation from spherical triangle to planar
triangle (since the exact reverse transformation equation for
Fuller's method is unknown). I have also changed Fuller's
orientation of some of the triangles, so that the antarctic circle
is not broken up.
*/
#include <cmath>
#include <sstream>
#include <vector>
using namespace std;
#include "Options.h"
#include "ProjectionIcosagnomonic.h"
#include "xpUtil.h"
// Aspect ratio of map.
#define TWIDE 5.25
#define THIGH 3
#define RATIO (((THIGH)*sqrt(3.0)/2)/(TWIDE))
/* signbit() appaired only in FreeBSD 5.1, so we need this hack */
#ifndef signbit
static int
signbit(double x)
{
if ((x < 0.0) || (x = -0.0))
return 1;
else
return 0;
}
#endif /* signbit */
bool
ProjectionIcosagnomonic::PointXY::sameSide(const PointXY& p1,
const PointXY& p2,
const PointXY& l1,
const PointXY& l2)
{
PointXY linec = l2 - l1;
PointXY pointc = p1 - l1;
double cp1 = linec.x*pointc.y - linec.y*pointc.x;
pointc = p2 - l1;
double cp2 = linec.x*pointc.y - linec.y*pointc.x;
return ((cp1 < 0 && cp2 < 0) || (cp1 > 0 && cp2 > 0));
}
bool
ProjectionIcosagnomonic::PointXY::inTriangle(const PointXY& a,
const PointXY& b,
const PointXY& c)
{
return (sameSide(*this, a, b, c) &&
sameSide(*this, b, a, c) &&
sameSide(*this, c, a, b));
}
double
ProjectionIcosagnomonic::PointXY::distance(const PointXY& p) const {
PointXY q = p-*this;
return sqrt(q.x*q.x+q.y*q.y);
}
// Rotate counterclockwise around the origin by the specified angle.
void
ProjectionIcosagnomonic::PointXY::rotate(double phi) {
// We convert to a vector, add the angle, and go back.
double mag = sqrt(x*x + y*y);
double theta = atan2(y, x);
theta += phi;
x = mag * cos(theta);
y = mag * sin(theta);
}
int
ProjectionIcosagnomonic::PointXY::cmpLine(const PointXY& a, const PointXY& b) {
double v = (b.x - a.x)*(y - a.y) - (x - a.x)*(b.y - a.y);
if (v < 0) return -1;
if (v > 0) return 1;
return 0;
}
ProjectionIcosagnomonic::PointXYZ
ProjectionIcosagnomonic::PointXYZ::crossP(const PointXYZ& a,
const PointXYZ& b) {
return PointXYZ(a.y*b.z - b.y*a.z, b.x*a.z - a.x*b.z, a.x*b.y - b.x*a.y);
}
double
ProjectionIcosagnomonic::PointXYZ::dotP(const PointXYZ& a,
const PointXYZ& b) {
return a.x*b.x + a.y*b.y + a.z*b.z;
}
bool
ProjectionIcosagnomonic::PointLL::sameSide(const PointLL& p1,
const PointLL& p2,
const PointLL& l1,
const PointLL& l2)
{
PointXYZ cp = PointXYZ::crossP(l1, l2);
double dp1 = PointXYZ::dotP(cp, p1);
double dp2 = PointXYZ::dotP(cp, p2);
return (signbit(dp1) == signbit(dp2) || fabs(dp1) < 1e-10);
}
bool
ProjectionIcosagnomonic::PointLL::inTriangle(const PointLL& a,
const PointLL& b,
const PointLL& c)
{
return (sameSide(*this, a, b, c) &&
sameSide(*this, b, a, c) &&
sameSide(*this, c, a, b));
}
double
ProjectionIcosagnomonic::PointLL::bearing(const PointLL& b) const {
const PointLL & a = *this;
double y = sin(b.lon-a.lon)*cos(b.lat);
double x = cos(a.lat)*sin(b.lat) - sin(a.lat)*cos(b.lat)*cos(b.lon-a.lon);
return atan2(y, x);
}
void
ProjectionIcosagnomonic::Polygon::scale(double x, double y) {
for (PointList::iterator i = V.begin(); i != V.end(); i ++) {
i->x *= x;
i->y *= y;
}
}
bool
ProjectionIcosagnomonic::Polygon::contains(PointXY p) const {
// Winding number algorithm from
// http://geometryalgorithms.com/Archive/algorithm_0103/algorithm_0103.htm
int wn = 0;
// loop through all edges of the polygon
for (PointList::const_iterator i = V.begin(); i != V.end(); i ++) {
PointList::const_iterator j = i;
j ++;
if (i->y <= p.y) { // start y <= P.y
if (j->y > p.y) // an upward crossing
if (p.cmpLine(*i, *j) > 0) // P left of edge
++wn; // have a valid up intersect
}
else { // start y > P.y (no test needed)
if (j->y <= p.y) // a downward crossing
if (p.cmpLine(*i, *j) < 0) // P right of edge
--wn; // have a valid down intersect
}
}
return wn;
}
ProjectionIcosagnomonic::Triangle::Triangle(PointXY a, PointXY b, PointXY c,
PointLL la, PointLL lb, PointLL lc)
: cA(a), cB(b), cC(c), sA(la), sB(lb), sC(lc) {
centerXY = PointXY((a.x + b.x + c.x)/3, (a.y + b.y + c.y)/3);
centerLL = sCentroid(la, lb, lc);
Options o;
o.Latitude(centerLL.lat);
o.Longitude(centerLL.lon);
double dim = cB.distance(cA);
o.CenterX(0);
o.CenterY(0);
// Reversed because (0,0) in upper left.
PointXY dA(cA.x-centerXY.x, centerXY.y-cA.y);
double d1 = atan2(dA.y, dA.x) - M_PI/2;// Angle to make map N-S.
double d2 = centerLL.bearing(la); // Angle to make point A go north.
rotation = d1+d2;
o.AddProjectionParameter(33.5*deg_to_rad);
P = new ProjectionGnomonic(1, (int)dim, (int)dim, &o);
}
ProjectionIcosagnomonic::Triangle::Triangle(const Triangle& T)
: cA(T.cA), cB(T.cB), cC(T.cC),sA(T.sA), sB(T.sB), sC(T.sC),
centerXY(T.centerXY), centerLL(T.centerLL),
rotation(T.rotation), P(new ProjectionGnomonic(*(T.P))) {}
ProjectionIcosagnomonic::Triangle::~Triangle() {
delete P;
}
/* Finds the spherical centroid of the triangle, by converting to Cartesian
* coordinates, finding the Cartesian centroid, normalizing to radius 1,
* and converting back to lat/lon.
*
* This algorithm is inexact in general, but I think it works for
* equilateral spherical triangles, which is what we're dealing with. */
ProjectionIcosagnomonic::PointLL
ProjectionIcosagnomonic::Triangle::sCentroid(PointLL a, PointLL b, PointLL c) {
PointXYZ ca(a);
PointXYZ cb(b);
PointXYZ cc(c);
PointXYZ center1 = PointXYZ((ca.x+cb.x+cc.x)/3,
(ca.y+cb.y+cc.y)/3,
(ca.z+cb.z+cc.z)/3);
double scale = sqrt(center1.x*center1.x+
center1.y*center1.y+
center1.z*center1.z);
PointXYZ center2 = PointXYZ(center1.x / scale,
center1.y / scale,
center1.z / scale);
return PointLL(center2);
}
bool
ProjectionIcosagnomonic::Triangle::contains(PointXY p) const {
return p.inTriangle(cA, cB, cC);
}
bool
ProjectionIcosagnomonic::Triangle::contains(PointLL p) const {
return p.inTriangle(sA, sB, sC);
}
bool
ProjectionIcosagnomonic::Triangle::pixelToSpherical(PointXY p,
double &lon,
double &lat) const {
if (!contains(p))
return false;
PointXY q(p - centerXY);
q.rotate(rotation);
return P->pixelToSpherical(q.x, q.y, lon, lat);
}
bool
ProjectionIcosagnomonic::Triangle::sphericalToPixel(double lon,
double lat,
double &x,
double &y) const {
if (!Triangle::contains(PointLL(lat, lon)))
return false;
PointXY p;
P->sphericalToPixel(lon, lat, p.x, p.y);
p.rotate(-rotation);
p += centerXY;
x = p.x;
y = p.y;
return contains(PointXY(x, y));
}
bool
ProjectionIcosagnomonic::ClippedTriangle::contains(PointXY p) const {
return Triangle::contains(p) && clip.contains(p);
}
bool
ProjectionIcosagnomonic::ClippedTriangle::contains(PointLL p) const {
double x, y;
bool success = Triangle::sphericalToPixel(p.lon, p.lat, x, y);
return success && contains(PointXY(x, y));
}
ProjectionIcosagnomonic::ProjectionIcosagnomonic(const int f, const int w,
const int h)
: ProjectionBase(f, w, h)
{
isWrapAround_ = false;
if (w * RATIO > h) {
_w = h / RATIO;
_h = h;
} else {
_w = w;
_h = w * RATIO;
}
baseh = _h / THIGH;
basew = _w / TWIDE;
offset = PointXY(centerX_ - _w/2, centerY_ - _h/2);
makeTriangles();
}
ProjectionIcosagnomonic::~ProjectionIcosagnomonic() {
for (TList::iterator i = T.begin(); i != T.end(); i ++) {
delete *i;
}
}
static int posmin(int a, int b, int c) {
if (a < b && a < c && a >= 0) {
return a;
} else if (b < c && b >= 0) {
return b;
} else if (c >= 0) {
return c;
} else
return 0;
}
void
ProjectionIcosagnomonic::makeTriangles() {
// Relates icosahedron vertices to points on the sphere.
static const PointLL icosa[] =
{ PointLL(+64.700000 * deg_to_rad, 10.536200 * deg_to_rad), // 0
PointLL( +2.300882 * deg_to_rad, -5.245390 * deg_to_rad), // 1
PointLL(+10.447378 * deg_to_rad, 58.157710 * deg_to_rad), // 2
PointLL(+39.100000 * deg_to_rad, 122.300000 * deg_to_rad), // 3
PointLL(+50.103201 * deg_to_rad, -143.478490 * deg_to_rad), // 4
PointLL(+23.717925 * deg_to_rad, -67.132330 * deg_to_rad), // 5
PointLL(-39.100000 * deg_to_rad, -57.700000 * deg_to_rad), // 6
PointLL(-50.103200 * deg_to_rad, 36.521500 * deg_to_rad), // 7
PointLL(-23.717925 * deg_to_rad, 112.867670 * deg_to_rad), // 8
PointLL( -2.300900 * deg_to_rad, 174.754600 * deg_to_rad), // 9
PointLL(-10.447345 * deg_to_rad, -121.842290 * deg_to_rad), // 10
PointLL(-64.700000 * deg_to_rad, -169.463800 * deg_to_rad), // 11
PointLL(+70.750433 * deg_to_rad, 26.020416 * deg_to_rad), // 12
PointLL(+60.846374 * deg_to_rad, 43.157248 * deg_to_rad), // 13
PointLL( +0.000000 * deg_to_rad, 0.000000 * deg_to_rad) };// 14
// Relates plane vertices to icosahedron vertices and to map coordinates.
static const struct { PointXY p; int v; } rect[] =
{ { PointXY(0.5, 0), 7 }, // 0
{ PointXY(1.5, 0), 1 }, // 1
{ PointXY(2.5, 0), 5 }, // 2
{ PointXY(3.5, 0), 1 }, // 3
{ PointXY(4.5, 0), 1 }, // 4
{ PointXY(0 , 1), 7 }, // 5
{ PointXY(1 , 1), 2 }, // 6
{ PointXY(2 , 1), 0 }, // 7
{ PointXY(3 , 1), 5 }, // 8
{ PointXY(4 , 1), 6 }, // 9
{ PointXY(5 , 1), 7 }, // 10
{ PointXY(-.5, 2), 7 }, // 11
{ PointXY(0.5, 2), 8 }, // 12
{ PointXY(1.5, 2), 3 }, // 13
{ PointXY(2.5, 2), 4 }, // 14
{ PointXY(3.5, 2), 10 }, // 15
{ PointXY(4.5, 2), 11 }, // 16
{ PointXY(5.5, 2), 8 }, // 17
{ PointXY(0 , 3), 11 }, // 18
{ PointXY(1 , 3), 9 }, // 19
{ PointXY(2 , 3), 9 }, // 20
{ PointXY(3 , 3), 9 }, // 21
{ PointXY(4 , 3), 9 }, // 22
{ PointXY(5 , 3), 9 }, // 23
{ PointXY(1 , 3), 8 } }; // 24 -- special duplicate of 19.
// Relates triangles to their planar vertices.
const int num_triangles = 23;
// clip is clipping polygon number or -1 for whole triangle.
static const struct { int v1, v2, v3, clip; } tri[] =
{ { 0, 1, 6, -1 }, // 0
{ 1, 6, 7, -1 }, // 1
{ 1, 2, 7, -1 }, // 2
{ 3, 8, 9, -1 }, // 3
{ 4, 9, 10, -1 }, // 4
{ 5, 6, 12, -1 }, // 5
{ 6, 12, 13, -1 }, // 6
{ 6, 7, 13, -1 }, // 7
{ 7, 13, 14, -1 }, // 8
{ 7, 8, 14, -1 }, // 9
{ 8, 14, 15, -1 }, // 10
{ 8, 9, 15, -1 }, // 11
{ 9, 10, 16, -1 }, // 12
{ 9, 15, 16, -1 }, // 13
{ 10, 16, 17, 3 }, // 14
{ 11, 12, 18, 2 }, // 15
{ 12, 18, 19, 4 }, // 16
{ 12, 13, 19, 0 }, // 17
{ 13, 14, 20, -1 }, // 18
{ 13, 24, 20, 1 }, // 19
{ 14, 15, 21, -1 }, // 20
{ 15, 16, 22, -1 }, // 21
{ 16, 17, 23, 5 } }; // 22
for (int i = 0; i < num_triangles; i ++) {
int pv1 = tri[i].v1;
int pv2 = tri[i].v2;
int pv3 = tri[i].v3;
PointXY pp1 = PointXY(rect[pv1].p.x*basew, rect[pv1].p.y*baseh);
PointXY pp2 = PointXY(rect[pv2].p.x*basew, rect[pv2].p.y*baseh);
PointXY pp3 = PointXY(rect[pv3].p.x*basew, rect[pv3].p.y*baseh);
PointLL ip1 = icosa[rect[pv1].v];
PointLL ip2 = icosa[rect[pv2].v];
PointLL ip3 = icosa[rect[pv3].v];
if (tri[i].clip >= 0) {
Polygon clip = makeClippingPath(tri[i].clip);
clip.scale(basew, baseh);
T.push_back(new ClippedTriangle(pp1, pp2, pp3,
ip1, ip2, ip3, clip));
} else {
T.push_back(new Triangle(pp1, pp2, pp3, ip1, ip2, ip3));
}
int findx = posmin((int)(rect[pv1].p.x*2),
(int)(rect[pv2].p.x*2),
(int)(rect[pv3].p.x*2));
int findy = posmin((int)rect[pv1].p.y,
(int)rect[pv2].p.y,
(int)rect[pv3].p.y);
// Triangles always span exactly two 'x' and one 'y'
// (not including clipping).
TFinder[findx][findy].push_back(T.back());
TFinder[findx+1][findy].push_back(T.back());
}
}
ProjectionIcosagnomonic::Polygon
ProjectionIcosagnomonic::makeClippingPath(unsigned int n) {
Polygon::PointList pp;
switch(n) {
case 0:
pp.push_back(PointXY(0.5, 2));
pp.push_back(PointXY(1, 3));
pp.push_back(PointXY(1, 2 + 1/3.));
pp.push_back(PointXY(1.5, 2));
break;
case 1:
pp.push_back(PointXY(1.5, 2));
pp.push_back(PointXY(1.5, 2 + 2/3.));
pp.push_back(PointXY(2, 3));
break;
case 2:
pp.push_back(PointXY(0, 2));
pp.push_back(PointXY(0.25, 2.5));
pp.push_back(PointXY(0.5, 2));
break;
case 3:
pp.push_back(PointXY(5, 1));
pp.push_back(PointXY(5.25, 1.5));
pp.push_back(PointXY(5, 2));
pp.push_back(PointXY(4.5, 2));
break;
case 4:
pp.push_back(PointXY(0.5, 2));
pp.push_back(PointXY(0.25, 2.5));
pp.push_back(PointXY(0.3, 3));
pp.push_back(PointXY(1, 3));
break;
case 5:
pp.push_back(PointXY(4.5, 2));
pp.push_back(PointXY(5, 2));
pp.push_back(PointXY(4.65, 2.3));
break;
default: break;
}
Polygon p(pp);
return p;
}
bool
ProjectionIcosagnomonic::pixelToSpherical(const double x, const double y,
double &lon, double &lat)
{
PointXY p(x, y);
p -= offset;
int fx = (int)(p.x / (basew/2));
int fy = (int)(p.y / baseh);
if (fx < 0 || fy < 0 || fx >= FINDCOLS || fy >= FINDROWS)
return false;
const TList& box = TFinder[fx][fy];
for (TList::const_iterator i = box.begin(); i != box.end(); i ++) {
if ((*i)->pixelToSpherical(p, lon, lat)) {
if (rotate_) RotateXYZ(lat, lon);
if (lon > M_PI) lon -= TWO_PI;
else if (lon < -M_PI) lon += TWO_PI;
return true;
}
}
return(false);
}
bool
ProjectionIcosagnomonic::sphericalToPixel(double lon, double lat,
double &x, double &y) const
{
if (rotate_) RotateZYX(lat, lon);
PointLL p(lat, lon);
for (TList::const_iterator i = T.begin(); i != T.end(); i ++) {
if ((*i)->contains(p)) {
if ((*i)->sphericalToPixel(lon, lat, x, y)) {
x += offset.x;
y += offset.y;
return true;
}
}
}
return(false);
}
syntax highlighted by Code2HTML, v. 0.9.1