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.

Sign up for GitHub

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

Jump to bottom

DM-12924: SpherePoint.offset should work at the poles and for negative offsets #302

Merged
merged 1 commit into from
Jan 8, 2018
Merged

DM-12924: SpherePoint.offset should work at the poles and for negative offsets #302

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
73 changes: 34 additions & 39 deletions include/lsst/afw/geom/SpherePoint.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,8 @@ class SpherePoint final {
/**
* Construct a SpherePoint from a longitude and latitude.
*
* @param longitude The longitude of the point.
* @param longitude The longitude of the point. `+/-inf` is treated as `nan`
* (because longitude is wrapped to a standard range and infinity cannot be wrapped)
* @param latitude The latitude of the point. Must be in the
* interval [-π/2, π/2] radians.
*
Expand All @@ -76,7 +77,8 @@ class SpherePoint final {
/**
* Construct a SpherePoint from double longitude and latitude.
*
* @param longitude The longitude of the point.
* @param longitude The longitude of the point. `+/-inf` is treated as `nan`
* (because longitude is wrapped to a standard range and infinity cannot be wrapped)
* @param latitude The latitude of the point. Must be in the
* interval [-π/2, π/2] radians.
* @param units The units of longitude and latitude
Expand All @@ -99,6 +101,9 @@ class SpherePoint final {
* Thrown if `vector` is the zero vector.
*
* @exceptsafe Provides strong exception guarantee.
*
* @note If the SpherePoint is at a pole then longitude is set to 0.
* That provides predictable behavior for @ref bearingTo and @ref offset.
*/
explicit SpherePoint(Point3D const& vector);

Expand Down Expand Up @@ -148,11 +153,6 @@ class SpherePoint final {
/**
* The longitude of this point.
*
* If this point is at a coordinate pole, the longitude is undefined, and
* this method may return any value. If the SpherePoint implementation
* allows multiple values of longitude from a pole, they shall all be
* treated as valid representations of the same point.
*
* @returns the longitude, in the interval [0, 2π) radians.
*
* @exceptsafe Shall not throw exceptions.
Expand Down Expand Up @@ -202,10 +202,8 @@ class SpherePoint final {
* @exceptsafe Shall not throw exceptions.
*/
bool atPole() const noexcept {
// Unit tests indicate we don't need to worry about epsilon-errors
// Objects constructed from lat=90*degrees, <0,0,1>, etc. all have
// atPole() = true. More complex operations like bearingTo have also
// been tested near the poles with no ill effects
// Unit tests indicate we don't need to worry about epsilon-errors, in that
// Objects constructed from lat=90*degrees, <0,0,1>, etc. all have atPole() = true.
return fabs(_latitude) >= HALFPI;
}

Expand All @@ -225,17 +223,18 @@ class SpherePoint final {
*/

/**
* `true` if two points represent the same position.
*
* Points are considered equal if and only if they represent the same
* location, regardless of how they were constructed. In particular,
* representations of the north or south pole with different longitudes
* are considered equal.
* `true` if two points have the same longitude and latitude.
*
* @param other the point to test for equality
* @returns true if this point matches `other` exactly, false otherwise
* @returns true if this point has exactly the same values of getLongitude() and getLatitude()
* as `other`, false otherwise
*
* @exceptsafe Shall not throw exceptions.
*
* @note Two points at the same pole will compare unequal if they have different longitudes,
* despite representing the same point on the unit sphere.
* This is important because the behavior of @ref offset and @ref bearingTo
* depend on longitude, even at the pole.
*
* @warning Points whose @ref getLongitude(), @ref getLatitude(), or
* @ref getVector() methods return `NaN` shall not compare
Expand All @@ -254,20 +253,19 @@ class SpherePoint final {
bool operator!=(SpherePoint const& other) const noexcept;

/**
* Direction from one point to another.
* Orientation at this point of the great circle arc to another point.
*
* The orientation is measured in a plane tangent to the sphere at this point, with 0 degrees along
* the direction of increasing longitude and 90 degrees along the direction of increasing latitude.
* Thus for any point except the pole, the arc is due east at 0 degrees and due north at 90 degrees.
* To understand the behavior at the poles it may help to imagine the point has been shifted slightly
* by decreasing the absolute value of its latitude.
Copy link
Member

Choose a reason for hiding this comment

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

I think this makes it pretty clear how the pole is handled. 👍

*
* This method finds the shortest (great-circle) arc between two
* points, and characterizes its direction by the angle between
* it and a line of latitude passing through this point. 0 represents
* due east, &pi;/2 represents due north. If the points are on
* opposite sides of the sphere, the bearing may be any value.
* If the points are on opposite sides of the sphere then the returned bearing may be any value.
*
* @param other the point to which to measure the bearing
* @returns the direction, as defined above, in the interval [0, 2&pi;).
*
* @throws pex::exceptions::DomainError
* Thrown if `this.atPole()`.
*
Copy link
Member

Choose a reason for hiding this comment

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

Please add a paragraph documenting the behavior at the poles. The current description implies that the bearing must always be -π/2 at the north pole and π/2 at the south pole, which I suspect is no longer what happens.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

OK. I rewrote the description to one whose behavior is clearer at the poles.

* @exceptsafe Provides strong exception guarantee.
*
* @note For two points `A` and `B`, `A.bearingTo(B)` will in
Expand Down Expand Up @@ -296,7 +294,8 @@ class SpherePoint final {
* @param axis a point defining the north pole of the rotation axis.
* @param amount the amount of rotation, where positive values
* represent right-handed rotations around `axis`.
* @returns a new point created by rotating this point
* @returns a new point created by rotating this point. If the new point is at the pole
* then its longitude will be 0.
*
* @exceptsafe Shall not throw exceptions.
*/
Expand All @@ -305,20 +304,16 @@ class SpherePoint final {
/**
* Return a point offset from this one along a great circle.
*
* For any point `A` not at a coordinate pole, and any two angles `b`
* and `delta`, `A.bearingTo(A.offset(b, delta))` = `b` and
* `A.separationTo(A.offset(b, delta))` = `delta`.
* For any point `A`, any angle `bearing` and any non-negative angle `amount`,
* if `B` = A.offset(bearing, amount)` then `A.bearingTo(B) = amount` and `A.separationTo(B) = amount`.
* Negative values of `amount` are supported in the obvious manner:
* `A.offset(bearing, delta) = A.offset(bearing + 180*degrees, -delta)`
*
* @param bearing the direction in which to move this point, following
* the conventions described in @ref bearingTo.
* @param amount the distance by which to move along the great
* circle defined by `bearing`
* @returns a new point created by shifting this point
*
* @throws pex::exceptions::DomainError
* Thrown if `this.atPole()`.
* @throws pex::exceptions::InvalidParameterError
* Thrown if `amount` is negative.
* @param amount the distance by which to move along the great circle defined by `bearing`
* @returns a new point created by rotating this point. If the new point is at the pole
* then its longitude will be 0.
*
* @exceptsafe Provides strong exception guarantee.
*/
Expand Down
32 changes: 8 additions & 24 deletions src/geom/SpherePoint.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,10 +90,14 @@ SpherePoint::SpherePoint(Point3D const& vector) {
double const y = vector.getY() / norm;
double const z = vector.getZ() / norm;

// Need to convert to Angle, Angle::wrap, and convert back to radians
// to handle _longitude = -1e-16 without code duplication
_longitude = (atan2(y, x) * radians).wrap().asRadians();
_latitude = asin(z);
if (!atPole()) {
// Need to convert to Angle, Angle::wrap, and convert back to radians
// to handle _longitude = -1e-16 without code duplication
_longitude = (atan2(y, x) * radians).wrap().asRadians();
} else {
_longitude = 0;
}
}

SpherePoint::SpherePoint(SpherePoint const& other) noexcept = default;
Expand Down Expand Up @@ -124,22 +128,12 @@ bool SpherePoint::isFinite() const noexcept { return isfinite(_longitude) && isf
bool SpherePoint::operator==(SpherePoint const& other) const noexcept {
// Deliberate override of Style Guide 5-12
// Approximate FP comparison would make object equality intransitive
if (atPole()) {
return _latitude == other._latitude;
} else {
return _longitude == other._longitude && _latitude == other._latitude;
}
return _longitude == other._longitude && _latitude == other._latitude;
}

bool SpherePoint::operator!=(SpherePoint const& other) const noexcept { return !(*this == other); }

Angle SpherePoint::bearingTo(SpherePoint const& other) const {
if (atPole()) {
stringstream buffer;
buffer << "Cannot calculate offset from pole " << *this << ".";
throw pexExcept::DomainError(buffer.str());
}

Copy link
Member

Choose a reason for hiding this comment

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

As I understand it, bearingTo can give different results for the same point (e.g., the north pole) depending on how the longitude has been initialized. If this is the desired behavior, consider adding a warning to the documentation to make this explicit.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I think the documentation is sufficiently clear, now that I have reworked it.

Angle const deltaLon = other.getLongitude() - this->getLongitude();

double const sinDelta1 = sin(getLatitude().asRadians());
Expand Down Expand Up @@ -167,16 +161,6 @@ SpherePoint SpherePoint::rotated(SpherePoint const& axis, Angle const& amount) c

SpherePoint SpherePoint::offset(Angle const& bearing, Angle const& amount) const {
double const phi = bearing.asRadians();
if (atPole()) {
stringstream buffer;
buffer << "Cannot define offset direction from pole " << *this << ".";
throw pexExcept::DomainError(buffer.str());
}
if (amount < 0.0) {
stringstream buffer;
buffer << "Negative offset of " << amount.asDegrees() << " degrees is not allowed.";
throw pexExcept::InvalidParameterError(buffer.str());
}

// let v = vector in the direction bearing points (tangent to surface of sphere)
// To do the rotation, use rotate() method.
Expand Down