Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add EllipsoidRhumb class similar to EllipsoidGeodesic #7484

Merged
merged 12 commits into from
Jan 18, 2019
Merged
154 changes: 123 additions & 31 deletions Source/Core/EllipsoidRhumbLine.js
Original file line number Diff line number Diff line change
Expand Up @@ -98,17 +98,13 @@ define([
return Math.log(Math.tan(0.5 * (CesiumMath.PI_OVER_TWO + latitude))) - (ellipticity / 2.0 * Math.log((1 + eSinL) / (1 - eSinL)));
}

function calculateHeading(ellipsoidRhumbLine, major, minor, firstLongitude, firstLatitude, secondLongitude, secondLatitude) {
function calculateHeading(ellipsoidRhumbLine, firstLongitude, firstLatitude, secondLongitude, secondLatitude) {
var sigma1 = calculateSigma(ellipsoidRhumbLine._ellipticity, firstLatitude);
var sigma2 = calculateSigma(ellipsoidRhumbLine._ellipticity, secondLatitude);
return Math.atan2(CesiumMath.negativePiToPi(secondLongitude - firstLongitude), sigma2 - sigma1);
}

function calculateArcLength(ellipsoidRhumbLine, major, minor, firstLongitude, firstLatitude, secondLongitude, secondLatitude) {
//>>includeStart('debug', pragmas.debug);
Check.defined('heading', ellipsoidRhumbLine._heading);
//>>includeEnd('debug');

var heading = ellipsoidRhumbLine._heading;
var deltaLongitude = secondLongitude - firstLongitude;

Expand All @@ -135,7 +131,7 @@ define([
var scratchCart1 = new Cartesian3();
var scratchCart2 = new Cartesian3();

function initiailize(ellipsoidRhumbLine, start, end, ellipsoid) {
function initialize(ellipsoidRhumbLine, start, end, ellipsoid) {
var major = ellipsoid.maximumRadius;
var minor = ellipsoid.minimumRadius;
var majorSquared = major * major;
Expand All @@ -161,10 +157,9 @@ define([
Check.typeOf.number.greaterThanOrEquals('value', Math.abs(Math.abs(Cartesian3.angleBetween(firstCartesian, lastCartesian)) - Math.PI), 0.0125);
//>>includeEnd('debug');

initiailize(ellipsoidRhumbLine, start, end, ellipsoid);
initialize(ellipsoidRhumbLine, start, end, ellipsoid);

ellipsoidRhumbLine._heading = calculateHeading(ellipsoidRhumbLine, ellipsoid.maximumRadius, ellipsoid.minimumRadius,
start.longitude, start.latitude, end.longitude, end.latitude);
ellipsoidRhumbLine._heading = calculateHeading(ellipsoidRhumbLine, start.longitude, start.latitude, end.longitude, end.latitude);
ellipsoidRhumbLine._distance = calculateArcLength(ellipsoidRhumbLine, ellipsoid.maximumRadius, ellipsoid.minimumRadius,
start.longitude, start.latitude, end.longitude, end.latitude);
}
Expand Down Expand Up @@ -271,7 +266,7 @@ define([
* @param {Cartographic} start The initial planetodetic point on the path.
* @param {Cartographic} end The final planetodetic point on the path.
* @param {Ellipsoid} [ellipsoid=Ellipsoid.WGS84] The ellipsoid on which the rhumb line lies.
* @param {EllipsoidRhumbLine} [result] The object in whice to store the result.
* @param {EllipsoidRhumbLine} [result] The object in which to store the result.
* @returns {EllipsoidRhumbLine} The EllipsoidRhumbLine object.
*/
EllipsoidRhumbLine.fromStartAndEnd = function(start, end, ellipsoid, result) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is this function needed, it just matches the constructor, doesn't it??

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added it because it allows using a scratch result object.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Seems like overkill, can't you just call setEndPoints on an existing object? I doubt someone would need to change the ellipsoid on an existing instance. If we have a valid use case, then maybe it's worth it, but we try not to add code "just because it might be useful" since that introduces more code to maintain.

Expand Down Expand Up @@ -305,27 +300,21 @@ define([
Check.defined('start', start);
Check.defined('heading', heading);
Check.defined('distance', distance);
Check.typeOf.number.greaterThan('distance', distance, 0.0);
//>>includeEnd('debug');

var e = defaultValue(ellipsoid, Ellipsoid.WGS84);

if (defined(result)) {
result._ellipsoid = e;
initiailize(result, start, undefined, e);
result._heading = CesiumMath.negativePiToPi(heading);
result._distance = distance;

result._end = result.interpolateUsingSurfaceDistance(distance, new Cartographic());
return result;
if (!defined(result)) {
result = new EllipsoidRhumbLine(undefined, undefined, e);
}

var rhumbLine = new EllipsoidRhumbLine(undefined, undefined, e);
initiailize(rhumbLine, start, undefined, e);
rhumbLine._heading = CesiumMath.negativePiToPi(heading);
rhumbLine._distance = distance;
rhumbLine._end = rhumbLine.interpolateUsingSurfaceDistance(distance, new Cartographic());
initialize(result, start, undefined, e);
result._heading = CesiumMath.negativePiToPi(heading);
result._distance = distance;

return rhumbLine;
result._end = result.interpolateUsingSurfaceDistance(distance, new Cartographic());
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This whole function is kind of an odd implementation to me. Why not break out interpolateUsingSurfaceDistance into a standalone local private function that can be called with start, heading, distance, ellipsoid, us that to compute the end point then just create the fully contructed result. Using a half-initialized instance like this seems like something that could bite us.

Obviously we still expose the prototype version of interpolateUsingSurfaceDistance that just called the local one.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, this could work better. I'll make the change.

return result;
};

/**
Expand All @@ -347,29 +336,26 @@ define([
* Provides the location of a point at the indicated portion along the rhumb line.
*
* @param {Number} fraction The portion of the distance between the initial and final points.
* @param {Cartographic} result The object in which to store the result.
* @param {Cartographic} [result] The object in which to store the result.
* @returns {Cartographic} The location of the point along the rhumb line.
*/
EllipsoidRhumbLine.prototype.interpolateUsingFraction = function(fraction, result) {
//>>includeStart('debug', pragmas.debug);
Check.typeOf.number('distance', this._distance);
//>>includeEnd('debug');

return this.interpolateUsingSurfaceDistance(fraction * this._distance, result);
};

/**
* Provides the location of a point at the indicated distance along the rhumb line.
*
* @param {Number} distance The distance from the inital point to the point of interest along the rhumbLine.
* @param {Cartographic} result The object in which to store the result.
* @param {Cartographic} [result] The object in which to store the result.
* @returns {Cartographic} The location of the point along the rhumb line.
*
* @exception {DeveloperError} start and end must be set before calling function interpolateUsingSurfaceDistance
*/
EllipsoidRhumbLine.prototype.interpolateUsingSurfaceDistance = function(distance, result) {
//>>includeStart('debug', pragmas.debug);
Check.defined('distance', this._distance);
Check.defined('distance', distance);
Check.typeOf.number.greaterThan('this._distance', this._distance, CesiumMath.EPSILON12);
shehzan10 marked this conversation as resolved.
Show resolved Hide resolved
//>>includeEnd('debug');

var ellipsoid = this._ellipsoid;
Expand Down Expand Up @@ -427,5 +413,111 @@ define([
return new Cartographic(longitude, latitude, 0);
};

/**
* Provides the location of a point at the indicated longitude along the rhumb line.
*
* @param {Number} intersectionLongitude The longitude, in radians, at which to find the intersection point from the starting point using the heading.
* @param {Cartographic} [result] The object in which to store the result.
* @returns {Cartographic} The location of the intersection point along the rhumb line.
shehzan10 marked this conversation as resolved.
Show resolved Hide resolved
*
* @exception {DeveloperError} start and end must be set before calling function findIntersectionWithLongitude.
*/
EllipsoidRhumbLine.prototype.findIntersectionWithLongitude = function(intersectionLongitude, result) {
//>>includeStart('debug', pragmas.debug);
Check.defined('intersectionLongitude', intersectionLongitude);
shehzan10 marked this conversation as resolved.
Show resolved Hide resolved
Check.defined('distance', this._distance);
shehzan10 marked this conversation as resolved.
Show resolved Hide resolved
//>>includeEnd('debug');

var ellipticity = this._ellipticity;
var heading = this._heading;
var absHeading = Math.abs(heading);
var start = this._start;

intersectionLongitude = CesiumMath.negativePiToPi(intersectionLongitude);

if (!defined(result)) {
result = new Cartographic();
}

// If heading is -PI/2 or PI/2, this is an E-W rhumb line
// If heading is 0 or PI, this is an N-S rhumb line
if (Math.abs(CesiumMath.PI_OVER_TWO - absHeading) <= CesiumMath.EPSILON8) {
result.longitude = intersectionLongitude;
result.latitude = start.latitude;
result.height = 0;
return result;
} else if (CesiumMath.equalsEpsilon(Math.abs(CesiumMath.PI_OVER_TWO - absHeading), CesiumMath.PI_OVER_TWO, CesiumMath.EPSILON8)) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I assume these epsilons were from the ported code? (if not, how did we choose them?) Same for the rest of the PR.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, it was chosen as the epsilon at which Math.cos() returns 1.0 similar to http://help.agi.com/AGIComponents/html/F_AGI_Foundation_Trig_AngleEpsilon.htm.

if (CesiumMath.equalsEpsilon(intersectionLongitude, start.longitude, CesiumMath.EPSILON12)) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this be !CesiumMath.equalsEpsilon(...?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, this is correct. It's else if is checking if the line is N-S, and the nested if is checking if the longitude is the same, in which case there are infinite intersections.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmmm ok. So this case is "infinite intersections" and returns undefined, and the path below returns the intersection at the pole when considering the rhumbline to be an infinite spiral defined by a start and a heading? If so, that might be worth noting in the doc.

return undefined;
}

result.longitude = intersectionLongitude;
result.latitude = CesiumMath.PI_OVER_TWO * Math.sign(heading);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should the intersection latitude be limited by the rhumbline's start/end latitude?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I initially thought of doing that. I couldn't come up with a good reason to limit that functionality because technically rhumb lines are defined only by start and heading and are "infinite spirals" in both directions.

result.height = 0;
return result;
}

// Use iterative solver from Equation 9 from http://edwilliams.org/ellipsoid/ellipsoid.pdf
var phi1 = start.latitude;
var eSinPhi1 = ellipticity * Math.sin(phi1);
var leftComponent = Math.tan(0.5 * (CesiumMath.PI_OVER_TWO + phi1)) * Math.exp((intersectionLongitude - start.longitude) / Math.tan(heading));
var denominator = (1 + eSinPhi1) / (1 - eSinPhi1);

var newPhi = start.latitude;
var phi;
do {
phi = newPhi;
var eSinPhi = ellipticity * Math.sin(phi);
var numerator = (1 + eSinPhi) / (1 - eSinPhi);
newPhi = 2 * Math.atan(leftComponent * Math.pow(numerator / denominator, ellipticity / 2)) - CesiumMath.PI_OVER_TWO;
} while (!CesiumMath.equalsEpsilon(newPhi, phi, CesiumMath.EPSILON12));
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you know if there's a situation where this wouldn't converge, and we would want this to break after some number of iterations?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It will converge and it's converging pretty quickly in all the tests I did. In the test cases I ran, it took less than 20 iterations.


result.longitude = intersectionLongitude;
result.latitude = phi;
result.height = 0;
return result;
};

/**
* Provides the location of a point at the indicated latitude along the rhumb line.
*
* @param {Number} intersectionLatitude The latitude, in radians, at which to find the intersection point from the starting point using the heading.
* @param {Cartographic} [result] The object in which to store the result.
* @returns {Cartographic} The location of the intersection point along the rhumb line, undefined if there is no intersection, infinite intersections or if intersection point is not between start and end.
*
* @exception {DeveloperError} start and end must be set before calling function findIntersectionWithLongitude.
*/
EllipsoidRhumbLine.prototype.findIntersectionWithLatitude = function(intersectionLatitude, result) {
//>>includeStart('debug', pragmas.debug);
Check.defined('intersectionLatitude', intersectionLatitude);
Check.defined('distance', this._distance);
shehzan10 marked this conversation as resolved.
Show resolved Hide resolved
//>>includeEnd('debug');

var ellipticity = this._ellipticity;
var heading = this._heading;
var start = this._start;

// If start and end have same latitude, return undefined since it's either no intersection or infinite intersections
if (CesiumMath.equalsEpsilon(Math.abs(heading), CesiumMath.PI_OVER_TWO, CesiumMath.EPSILON8)) {
return;
}

// Can be solved using the same equations from interpolateUsingSurfaceDistance
var sigma1 = calculateSigma(ellipticity, start.latitude);
var sigma2 = calculateSigma(ellipticity, intersectionLatitude);
var deltaLongitude = Math.tan(heading) * (sigma2 - sigma1);
var longitude = CesiumMath.negativePiToPi(start.longitude + deltaLongitude);

if (defined(result)) {
result.longitude = longitude;
result.latitude = intersectionLatitude;
result.height = 0;

return result;
}

return new Cartographic(longitude, intersectionLatitude, 0);
};

return EllipsoidRhumbLine;
});
Loading