I'm a General! Wheeeeeeeeee!

I implemented the Lee Conformal Projection, and let me tell you, it took
a lot of work. Oh, also, I tried another family of pseudo-azimuthal-ish
projections that just kind of sucked. I think I'm going to try
optimizing some tetrahedral projections next. Hence the Lee projection.
I used a 28th-order McLaurin series for the Lee projection, so it
doesn't actually tesselate perfectly. I think I might have entered a
number wrong, because seriously?! The 28th order with six decimal places
of precision wasn't good enough? Anyway, the series was a suggestion
from Lee himself in the paper that I read on JStor, that I had to make a
JStor account for, because I couldn't get the equations literally
anywhere else on the internet. NGAHHH!
This commit is contained in:
Justin Kunimune 2017-05-27 20:27:38 -10:00
parent 2fcfa325e7
commit 6853216990
6 changed files with 157 additions and 60 deletions

Binary file not shown.

Before

Width:  |  Height:  |  Size: 22 KiB

After

Width:  |  Height:  |  Size: 25 KiB

View File

@ -334,9 +334,10 @@ public class MapAnalyzer extends Application {
final double s1ms2 = Math.hypot((p1[0]-p0[0])-(p2[1]-p0[1]), (p1[1]-p0[1])+(p2[0]-p0[0]));
final double factor = Math.abs((s1ps2-s1ms2)/(s1ps2+s1ms2)); //there's some linear algebra behind this formula. Don't worry about it.
output[1] = Math.asin(1-factor);
if (output[1] >= Math.PI/2)
output[1] = Double.NaN;
//output[1] = Math.asin(1-factor);
output[1] = (1-factor)/Math.sqrt(factor);
//if (output[1] >= Math.PI/2)
// output[1] = Double.NaN;
return output;
}
@ -356,14 +357,14 @@ public class MapAnalyzer extends Application {
final int r, g, b;
if (sizeDistort < 0) { //if compressing
r = (int)(256*(1-shapeDistort/(Math.PI/2)));
g = (int)(256*(1-shapeDistort/(Math.PI/2))*(1-Math.min(1, -sizeDistort/3)));
r = (int)(255.9*Math.exp(-shapeDistort/1.5));
g = (int)(255.9*Math.exp(-shapeDistort/1.5)*Math.exp(sizeDistort/1.5));
b = g;
}
else { //if dilating
r = (int)(256*(1-shapeDistort/(Math.PI/2))*(1-Math.min(1, sizeDistort/3)));
r = (int)(255.9*Math.exp(-shapeDistort/1.5)*Math.exp(-sizeDistort/1.5));
g = r;
b = (int)(256*(1-shapeDistort/(Math.PI/2)));
b = (int)(255.9*Math.exp(-shapeDistort/1.5));
}
final int argb = ((((((0xFF)<<8)+r)<<8)+g)<<8)+b;

View File

@ -22,8 +22,8 @@ public class MapOptimizer extends Application {
private static final String[] EXISTING_PROJECTIONS = { "Hobo-Dyer", "Winkel Tripel", "Robinson", "Van der Grinten",
"Mercator" };
private static final double[] WEIGHTS = {.22, .50, .86, 1.3, 2.0, 2.67, 3.0, 4.7, 8.0, 18.0 };
"Mercator" }; //"Pierce Quincuncial" };
private static final double[] WEIGHTS = { .25, .43, .67, 1.0, 1.5, 2.3, 4.0 };
private ScatterChart<Number, Number> chart;
@ -39,17 +39,17 @@ public class MapOptimizer extends Application {
double[][][] globe = MapAnalyzer.globe(0.02);
final Series<Number, Number> oldMaps = analyzeAll(globe, EXISTING_PROJECTIONS);
final Series<Number, Number> hyperMaps = optimizeHyperelliptical(globe);
final Series<Number, Number> roundMaps = optimizeElliptihypercosine(globe);
//final Series<Number, Number> hyperMaps = optimizeHyperelliptical(globe);
final Series<Number, Number> roundMaps = optimizeElliptabola(globe);
chart = new ScatterChart<Number, Number>(new NumberAxis("Size distortion", 0, 1, 0.1),
new NumberAxis("Shape distortion", 0, 0.5, 0.1));
chart = new ScatterChart<Number, Number>(new NumberAxis("Size distortion", 0, 1, 0.2),
new NumberAxis("Shape distortion", 0, 1, 0.2));
chart.getData().add(oldMaps);
chart.getData().add(hyperMaps);
//chart.getData().add(hyperMaps);
chart.getData().add(roundMaps);
System.out.println("Total time elapsed: "+
(System.currentTimeMillis()-startTime)/1000.);
(System.currentTimeMillis()-startTime)/1000.+"s");
stage.setTitle("Map Projections");
stage.setScene(new Scene(chart));
@ -84,7 +84,7 @@ public class MapOptimizer extends Application {
if (params[i] < bounds[i][1]) {
for (int j = 0; j < i; j ++)
params[j] = bounds[j][0];
params[i] += (bounds[i][1]-bounds[i][0])/8;
params[i] += bounds[i][2];
break;
}
}
@ -96,7 +96,7 @@ public class MapOptimizer extends Application {
final double avgSizeD = Stat.stdDev(dist[0]);
final double avgShapeD = Stat.mean(dist[1]);
for (int k = 0; k < WEIGHTS.length; k ++) {
final double avgDist = avgSizeD + WEIGHTS[k]*avgShapeD;
final double avgDist = Math.pow(avgSizeD,1.5) + WEIGHTS[k]*Math.pow(avgShapeD,1.5);
if (avgDist < currentBest[k][0]) {
currentBest[k][0] = avgDist;
currentBest[k][1] = avgSizeD;
@ -122,12 +122,14 @@ public class MapOptimizer extends Application {
private static Series<Number, Number> optimizeHyperelliptical(double[][][] points) { //optimize and plot some hyperelliptical maps
return optimizeFamily(MapOptimizer::hyperelliptical, "Hyperelliptic", new double[][] {{2,5}, {.75,1.75}}, points);
return optimizeFamily(MapOptimizer::hyperelliptical, "Hyperelliptic",
new double[][] {{2.5,5,.125}, {0.5,1.75,.125}, {1.0,1.5,.0625}}, points);
}
private static Series<Number, Number> optimizeElliptihypercosine(double[][][] points) { //optimize and plot some elliptical-hypebolic-cosine maps
return new Series<Number, Number>();
private static Series<Number, Number> optimizeElliptabola(double[][][] points) { //optimize and plot some elliptical-hypebolic-cosine maps
return optimizeFamily(MapOptimizer::elliptabola, "Elliptabola",
new double[][] {{.25,1,.125}, {.75,1.75,.125}, {1.25,2,.125}, {.25,1.5,.125}}, points);
}
@ -141,11 +143,30 @@ public class MapOptimizer extends Application {
private static final double[] hyperelliptical(double[] coords, double[] params) { //a hyperelliptic map projection with hyperellipse order k and lattitudinal spacind described by x^n/n
final double lat = coords[0], lon = coords[1];
final double k = params[0], n = params[1];
final double k = params[0], n = params[1], a = params[2];
return new double[] {
Math.pow(1 - Math.pow(Math.abs(lat/(Math.PI/2)), k),1/k)*lon,
(1-Math.pow(1-Math.abs(lat/(Math.PI/2)), n))/Math.sqrt(n)*Math.signum(lat)*Math.PI/2};
(1-Math.pow(1-Math.abs(lat/(Math.PI/2)), n))/Math.sqrt(n)*Math.signum(lat)*Math.PI/2*a};
}
private static double[] elliptabola(double[] coords, double[] params) {
final double lat = coords[0], lon = coords[1];
final double k1 = params[0], k2 = params[1], a = params[2], p2 = params[3];
final double f1p2 = Math.pow(1-Math.pow(1-Math.abs(lat)/(Math.PI/2), k1),2);
final double f2l2 = Math.pow(1-Math.pow(1-Math.abs(lon)/Math.PI, k2),2);
final double ax = f1p2/(p2*p2);
final double bx = 2*f1p2/(p2)+1/(f2l2);
final double cx = f1p2 - 1;
final double ay = f2l2;
final double by = p2/Math.sqrt(f1p2);
final double cy = -p2 - f2l2;
return new double[] {
Math.PI*Math.sqrt((-bx+Math.sqrt(bx*bx-4*ax*cx))/(2*ax))*Math.signum(lon),
Math.PI*(-by+Math.sqrt(by*by-4*ay*cy))/(2*ay)*Math.signum(lat)/a
};
}
}

View File

@ -42,6 +42,7 @@ import javafx.scene.text.Text;
import javafx.stage.FileChooser;
import javafx.stage.Stage;
import mfc.field.Complex;
import util.Dixon;
import util.ProgressBarDialog;
import util.Robinson;
import util.WinkelTripel;
@ -61,28 +62,24 @@ public class MapProjections extends Application {
private static final KeyCombination ctrlS = new KeyCodeCombination(KeyCode.S, KeyCodeCombination.CONTROL_DOWN);
private static final String[] PROJ_ARR = { "Equirectangular", "Mercator",
"Gall Stereographic", "Hobo-Dyer", "Polar", "Stereographic",
"Azimuthal Equal-Area", "Orthographic", "Gnomonic",
"Equidistant Conic", "Conformal Conic", "Albers", "Van der Grinten", "Robinson", "Winkel Tripel",
"Mollweide", "Hammer", "Sinusoidal", "Lemons",
"Pierce Quincuncial", "Guyou", "AuthaGraph", "TetraGraph",
"Magnifier", "Experimental" };
private static final double[] DEFA = { 2, 1, 4/3.0, 1.977, 1, 1, 1, 1, 1, 2,
2, 2, 1, 1.9716, Math.PI/2, 2, 2, 2, 2, 1, 2, 4.0 / Math.sqrt(3),
Math.sqrt(3), 1, 1 };
private static final String[] DESC = { "An equidistant cylindrical map", "A conformal cylindrical map",
"A compromising cylindrical map", "An equal-area cylindrical map", "An equidistant azimuthal map",
"A conformal azimuthal map", "An equal-area azimuthal map",
private static final String[] PROJ_ARR = { "Mercator", "Equirectangular", "Hobo-Dyer", "Gall Stereographic",
"Stereographic", "Polar", "Azimuthal Equal-Area", "Orthographic", "Gnomonic", "Conformal Conic",
"Equidistant Conic", "Albers", "Lee", "TetraGraph", "AuthaGraph", "Mollweide", "Hammer", "Sinusoidal", "Van der Grinten", "Robinson",
"Winkel Tripel", "Lemons", "Pierce Quincuncial", "Guyou", "Magnifier", "Experimental" };
private static final double[] DEFA = { 1, 2, 1.977, 4/3., 1, 1, 1, 1, 1, 2,
2, 2, Math.sqrt(3), Math.sqrt(3), 4/Math.sqrt(3), 2, 2, 2, 1, 1.9716, Math.PI/2, 2, 1, 2, 1, 1 };
private static final String[] DESC = { "A conformal cylindrical map", "An equidistant cylindrical map",
"An equal-area cylindrical map", "A compromising cylindrical map", "A conformal azimuthal map",
"An equidistant azimuthal map", "An equal-area azimuthal map",
"Represents earth viewed from an infinite distance",
"Every straight line on the map is a straight line on the sphere", "An equidistant conic map",
"A conformal conic map", "An equal-area conic map", "A circular compromise map",
"A circular compromise map", "A visually pleasing piecewise compromise map",
"The compromise map used by National Geographic", "An equal-area map shaped like an ellipse",
"An equal-area map shaped like an ellipse", "An equal-area map shaped like a sinusoid",
"Every straight line on the map is a straight line on the sphere", "A conformal conic map",
"An equidistant conic map", "An equal-area conic map", "A conformal tetrahedral map that really deserves more attention",
"An equidistant tetrahedral map that I invented", "An almost-equal-area tetrahedral map", "An equal-area map shaped like an ellipse",
"An equal-area map shaped like an ellipse", "An equal-area map shaped like a sinusoid","A circular compromise map",
"A visually pleasing piecewise compromise map",
"The compromise map used by National Geographic",
"BURN LIFE'S HOUSE DOWN!", "A conformal square map that uses complex math",
"A reorganized version of Pierce Quincuncial and actually the best map ever",
"An almost-equal-area map based on a tetrahedron.", "A compromising knockoff of the AuthaGraph projection",
"A rearranged Pierce Quincuncial map",
"A novelty map that swells the center to disproportionate scale",
"What happens when you apply a complex differentiable function to a stereographic projection?" };
@ -162,7 +159,7 @@ public class MapProjections extends Application {
}
}
});
projectionChooser.setValue(PROJ_ARR[1]);
projectionChooser.setValue("Mercator");
layout.getChildren().add(new HBox(3, lbl, projectionChooser));
projectionDesc = new Text(DESC[1]);
@ -465,6 +462,8 @@ public class MapProjections extends Application {
return albers(X, Y);
else if (p.equals("Robinson"))
return robinson(X, Y);
else if (p.equals("Lee"))
return lee(X, Y);
else
throw new IllegalArgumentException(p);
}
@ -526,7 +525,7 @@ public class MapProjections extends Application {
Complex k = new Complex(Math.sqrt(0.5)); // the rest comes from some fancy complex calculus
Complex ans = Jacobi.cn(u, k);
double p = 2 * Math.atan(ans.abs());
double theta = Math.atan2(ans.getIm(), ans.getRe()) - Math.PI/2;
double theta = ans.arg() - Math.PI/2;
double lambda = Math.PI/2 - p;
return new double[] {lambda, theta};
}
@ -535,7 +534,7 @@ public class MapProjections extends Application {
Complex z = new Complex(x*3, y*3);
Complex ans = z.sin();
double p = 2 * Math.atan(ans.abs());
double theta = Math.atan2(ans.getIm(), ans.getRe()) + Math.PI/2;
double theta = ans.arg();
double lambda = Math.PI/2 - p;
return new double[] {lambda, theta};
}
@ -619,11 +618,45 @@ public class MapProjections extends Application {
Complex k = new Complex(Math.sqrt(0.5)); // the rest comes from some fancy complex calculus
Complex ans = Jacobi.cn(u, k);
double p = 2 * Math.atan(ans.abs());
double theta = Math.atan2(ans.getIm(), ans.getRe());
double theta = ans.arg();
double lambda = Math.PI/2 - p;
return new double[] {lambda, theta};
}
private static double[] lee(double x, double y) { //a tessalatable rectangle map
final double x1, y1;
if (y > x+1) {
x1 = -x - 4/3.;
y1 = y - 1;
}
else if (y < -x-1) {
x1 = x + 2/3.;
y1 = -y - 1;
}
else if (y > 1-x) {
x1 = -4/3. + x;
y1 = 1 - y;
}
else if (y < x-1) {
x1 = 2/3. - x;
y1 = 1 + y;
}
else if (x < 0) {
x1 = x + 2/3.;
y1 = -y - 1;
}
else {
x1 = 2/3. - x;
y1 = y + 1;
}
Complex w = new Complex(x1, y1/Math.sqrt(3)).times(2.655);//2.65 2.66
Complex ans = Dixon.leeFunc(w).times(Math.pow(2, -5/6.));
double p = 2 * Math.atan(ans.abs());
double theta = ans.arg();
double lambda = p - Math.PI/2;
return new double[] {lambda, theta};
}
private static double[] mollweide(double x, double y) {
double tht = Math.asin(y);
return new double[] {
@ -632,8 +665,6 @@ public class MapProjections extends Application {
}
private static double[] winkel_tripel(double x, double y) {
final double tolerance = 10e-2;
double X = x * (1 + Math.PI/2);

33
src/util/Dixon.java Normal file
View File

@ -0,0 +1,33 @@
/**
*
*/
package util;
import mfc.field.Complex;
/**
* A class with a few handy Dixon elliptic functions as they pertain to the Lee conformal projection.
* All the algorithms here came directly from L.P. Lee's paper,
* "Some Conformal Projections Based on Elliptic Functions"
*
* Lee, L. P. Some Conformal Projections Based on Elliptic Functions.
* Geographical Review, vol. 55, no. 4, 1965, pp. 563580. JSTOR,
* www.jstor.org/stable/212415.
*
* @author jkunimune
*/
public class Dixon {
public static Complex leeFunc(Complex w) { //the 28th order McLaurin polynomial for 2sm(w/2)cm(w/2)
return w.minus(w.pow(4).times(.625e-1))
.plus(w.pow(7).times(.223214e-2))
.minus(w.pow(10).times(.0609754e-3))
.plus(w.pow(13).times(.020121e-4))
.minus(w.pow(16).times(.005539e-5))
.plus(w.pow(19).times(.001477e-6))
.minus(w.pow(22).times(.000385e-7))
.plus(w.pow(25).times(.000099e-8))
.minus(w.pow(28).times(.000025e-9));
}
}

View File

@ -165,7 +165,7 @@ public class MapProjections extends Application {
break;
}
}
updateMap();
updateMap(); //TODO: how cool would it be to animate all of the transitions?
}
});
projectionChooser.setPrefWidth(210);
@ -549,7 +549,7 @@ public class MapProjections extends Application {
else if (p.equals("Robinson"))
return robinson(lat, lon);
else if (p.equals("Test"))
return hyperelliptic(lat, lon);
return elliptabola(lat, lon);
else
throw new IllegalArgumentException(p);
}
@ -832,23 +832,34 @@ public class MapProjections extends Application {
}
private static double[] hyperelliptic(double lat, double lon) {
final double k = 3.5;
final double n = 1.3125;
private static final double[] hyperelliptical(double lat, double lon) { //a hyperelliptic map projection with hyperellipse order k and lattitudinal spacind described by x^n/n
final double k = 3.25, n = 1.25, a = 1.125;
return new double[] {
Math.pow(1 - Math.pow(Math.abs(lat/(Math.PI/2)), k),1/k)*lon,
(1-Math.pow(1-Math.abs(lat/(Math.PI/2)), n))/Math.sqrt(n)*Math.signum(lat)*Math.PI/2};
(1-Math.pow(1-Math.abs(lat/(Math.PI/2)), n))/Math.sqrt(n)*Math.signum(lat)*Math.PI/2*a};
}
/*private static double[] ellipticosh(double lat, double lon) {
final double a = 16/9.;
final double n = 1.5;
private static double[] elliptabola(double lat, double lon) {
final double k1 = 1.5;
final double k2 = .95;
final double a = 2;
final double p2 = 4;
final double f1p2 = Math.pow(1-Math.pow(1-Math.abs(lat)/(Math.PI/2), k1),2);
final double f2l2 = Math.pow(1-Math.pow(1-Math.abs(lon)/Math.PI, k2),2);
final double ax = f1p2/(p2*p2);
final double bx = 2*f1p2/(p2)+1/(f2l2);
final double cx = f1p2 - 1;
final double ay = f2l2;
final double by = p2/Math.sqrt(f1p2);
final double cy = -p2 - f2l2;
return new double[] {
0,
0
Math.PI*Math.sqrt((-bx+Math.sqrt(bx*bx-4*ax*cx))/(2*ax))*Math.signum(lon),
Math.PI*(-by+Math.sqrt(by*by-4*ay*cy))/(2*ay)*Math.signum(lat)/a
};
}*/
}