Geodesic intersections More...

#include <GeographicLib/Intersect.hpp>

Public Types
typedef std::pair< Math::real, Math::real >	Point

Public Member Functions
Constructor
	Intersect (const Geodesic &geod)

Finding intersections
Point	Closest (Math::real latX, Math::real lonX, Math::real aziX, Math::real latY, Math::real lonY, Math::real aziY, const Point &p0=Point(0, 0), int *c=nullptr) const

Point	Closest (const GeodesicLine &lineX, const GeodesicLine &lineY, const Point &p0=Point(0, 0), int *c=nullptr) const

Point	Segment (Math::real latX1, Math::real lonX1, Math::real latX2, Math::real lonX2, Math::real latY1, Math::real lonY1, Math::real latY2, Math::real lonY2, int &segmode, int *c=nullptr) const

Point	Segment (const GeodesicLine &lineX, const GeodesicLine &lineY, int &segmode, int *c=nullptr) const

Point	Next (Math::real latX, Math::real lonX, Math::real aziX, Math::real aziY, int *c=nullptr) const

Point	Next (const GeodesicLine &lineX, const GeodesicLine &lineY, int *c=nullptr) const

Finding all intersections
std::vector< Point >	All (Math::real latX, Math::real lonX, Math::real aziX, Math::real latY, Math::real lonY, Math::real aziY, Math::real maxdist, std::vector< int > &c, const Point &p0=Point(0, 0)) const

std::vector< Point >	All (Math::real latX, Math::real lonX, Math::real aziX, Math::real latY, Math::real lonY, Math::real aziY, Math::real maxdist, const Point &p0=Point(0, 0)) const

std::vector< Point >	All (const GeodesicLine &lineX, const GeodesicLine &lineY, Math::real maxdist, std::vector< int > &c, const Point &p0=Point(0, 0)) const

std::vector< Point >	All (const GeodesicLine &lineX, const GeodesicLine &lineY, Math::real maxdist, const Point &p0=Point(0, 0)) const

Diagnostic counters
long long	NumInverse () const

long long	NumBasic () const

long long	NumChange () const

long long	NumCorner () const

long long	NumOverride () const

Insepctor function
const Geodesic &	GeodesicObject () const

Static Public Member Functions
static Math::real	Dist (const Point &p, const Point &p0=Point(0, 0))

Static Public Attributes
static const unsigned	LineCaps

Detailed Description

Geodesic intersections

Find the intersections of two geodesics X and Y. Four calling sequences are supported.

The geodesics are defined by a position (latitude and longitude) and an azimuth. In this case the closest intersection is found.
The geodesics are defined by two endpoints. The intersection of the two segments is found. If they don't intersect, then the closest intersection is returned.
The geodesics are defined as an intersection point, a single position and two azimuths. In this case, the next closest intersection is found.
The geodesics are defined as in the first case and all intersection within a specified distance are returned.

In all cases the position of the intersection is given by the signed displacements x and y along the geodesics from the starting point (the first point in the case of a geodesic segment). The closest itersection is defined as the one that minimizes the L1 distance, Intersect::Dist([x, y) = |x| + |y|.

The routines also optionally return a coincidence indicator c. This is typically 0. However if the geodesics lie on top of one another at the point of intersection, then c is set to +1, if they are parallel, and −1, if they are antiparallel.

Example of use:

// Example of using the GeographicLib::Intersect class
 
#include <iostream>
#include <exception>
#include <GeographicLib/Intersect.hpp>
#include <GeographicLib/Constants.hpp>
 
using namespace std;
using namespace GeographicLib;
 
int main() {
  try {
    Geodesic geod(Constants::WGS84_a(), Constants::WGS84_f());
    Intersect intersect(geod);
    // Define the two geodesics
    GeodesicLine lineX(geod, 0, 0, 45, Intersect::LineCaps);
    GeodesicLine lineY(geod, 45, 10, 135, Intersect::LineCaps);
    // Find displacement to closest intersection
    Intersect::Point point = intersect.Closest(lineX, lineY);
    // Check position at intersection
    double latx, lonx, laty, lony;
    lineX.Position(point.first, latx, lonx);
    lineY.Position(point.second, laty, lony);
    cout << "X intersection displacement + position "
         << point.first << " " << latx << " " << lonx << "\n";
    cout << "Y intersection displacement + position "
         << point.second << " " << laty << " " << lony << "\n";
  }
  catch (const exception& e) {
    cerr << "Caught exception: " << e.what() << "\n";
    return 1;
  }
}

IntersectTool is a command-line utility providing access to the functionality of this class.

This solution for intersections is described in

C. F. F. Karney,
Geodesic intersections, J. Surveying Eng. 150(3), 04024005:1–9 (2024); preprint arxiv:2308.00495.

It is based on the work of

S. Baseldga and J. C. Martinez-Llario, Intersection and point-to-line solutions for geodesics on the ellipsoid, Stud. Geophys. Geod. 62, 353–363 (2018); DOI: 10.1007/s11200-017-1020-z.

Definition at line 71 of file Intersect.hpp.

Member Typedef Documentation

◆ Point

typedef std::pair<Math::real, Math::real> GeographicLib::Intersect::Point

The type used to hold the two displacement along the geodesics. This is just a std::pair with x = first and y = second.

Definition at line 79 of file Intersect.hpp.

Constructor & Destructor Documentation

◆ Intersect()

GeographicLib::Intersect::Intersect ( const Geodesic & geod )

Constructor for an ellipsoid with

Parameters

[in] geod a Geodesic object. This sets the parameters a and f for the ellipsoid.

Exceptions

GeographicErr if the eccentricity of the elliposdoid is too large.

Note: This class has been validated for -1/4 ≤ f ≤ 1/5. It may give satisfactory results slightly outside this range; however sufficient far outside the range, some internal checks will fail and an exception thrown.; If |f| > 1/50, then the Geodesic object should be constructed with exact = true.

Definition at line 20 of file Intersect.cpp.

References GeographicLib::Geodesic::Inverse(), and GeographicLib::Math::pi().

Member Function Documentation

◆ Closest() [1/2]

Intersect::Point GeographicLib::Intersect::Closest	(	Math::real	latX,
		Math::real	lonX,
		Math::real	aziX,
		Math::real	latY,
		Math::real	lonY,
		Math::real	aziY,
		const Point &	p0 = `Point(0, 0)`,
		int *	c = `nullptr`
	)		const

Find the closest intersection point, with each geodesic specified by position and azimuth.

Parameters

[in]	latX	latitude of starting point for geodesic X (degrees).
[in]	lonX	longitude of starting point for geodesic X (degrees).
[in]	aziX	azimuth at starting point for geodesic X (degrees).
[in]	latY	latitude of starting point for geodesic Y (degrees).
[in]	lonY	longitude of starting point for geodesic Y (degrees).
[in]	aziY	azimuth at starting point for geodesic Y (degrees).
[in]	p0	an optional offset for the starting points (meters), default = [0,0].
[out]	c	optional pointer to an integer coincidence indicator.

Returns: p the intersection point closest to p0.

The returned intersection minimizes Intersect::Dist(p, p0).

Definition at line 55 of file Intersect.cpp.

References Closest(), GeographicLib::Geodesic::Line(), and LineCaps.

Referenced by Closest(), and main().

◆ Closest() [2/2]

Intersect::Point GeographicLib::Intersect::Closest	(	const GeodesicLine &	lineX,
		const GeodesicLine &	lineY,
		const Point &	p0 = `Point(0, 0)`,
		int *	c = `nullptr`
	)		const

Find the closest intersection point, with each geodesic given as a GeodesicLine.

Parameters

[in]	lineX	geodesic X.
[in]	lineY	geodesic Y.
[in]	p0	an optional offset for the starting points (meters), default = [0,0].
[out]	c	optional pointer to an integer coincidence indicator.

Returns: p the intersection point closest to p0.

The returned intersection minimizes Intersect::Dist(p, p0).

Note: lineX and lineY should be created with minimum capabilities Intersect::LineCaps. The methods for creating a GeodesicLine include all these capabilities by default.

Definition at line 64 of file Intersect.cpp.

◆ Segment() [1/2]

Intersect::Point GeographicLib::Intersect::Segment	(	Math::real	latX1,
		Math::real	lonX1,
		Math::real	latX2,
		Math::real	lonX2,
		Math::real	latY1,
		Math::real	lonY1,
		Math::real	latY2,
		Math::real	lonY2,
		int &	segmode,
		int *	c = `nullptr`
	)		const

Find the intersection of two geodesic segments defined by their endpoints.

Parameters

[in]	latX1	latitude of first point for segment X (degrees).
[in]	lonX1	longitude of first point for segment X (degrees).
[in]	latX2	latitude of second point for segment X (degrees).
[in]	lonX2	longitude of second point for segment X (degrees).
[in]	latY1	latitude of first point for segment Y (degrees).
[in]	lonY1	longitude of first point for segment Y (degrees).
[in]	latY2	latitude of second point for segment Y (degrees).
[in]	lonY2	longitude of second point for segment Y (degrees).
[out]	segmode	an indicator equal to zero if the segments intersect (see below).
[out]	c	optional pointer to an integer coincidence indicator.

Returns: p the intersection point if the segments intersect, otherwise the intersection point closest to the midpoints of the two segments.

Warning: The results are only well defined if there's a unique shortest geodesic between the endpoints of the two segments.

segmode codes up information about the closest intersection in the case where the segments intersect. Let x₁₂ be the length of the segment X and x = p.first, the position of the intersection on segment X. Define

k_x = −1, if x < 0,
k_x = 0, if 0 ≤ x ≤ x₁₂,
k_x = 1, if x₁₂ < x.

and similarly for segment Y. Then segmode = 3 k_x + k_y.

Definition at line 72 of file Intersect.cpp.

References GeographicLib::Geodesic::InverseLine(), LineCaps, and Segment().

Referenced by main(), and Segment().

◆ Segment() [2/2]

Intersect::Point GeographicLib::Intersect::Segment	(	const GeodesicLine &	lineX,
		const GeodesicLine &	lineY,
		int &	segmode,
		int *	c = `nullptr`
	)		const

Find the intersection of two geodesic segments each defined by a GeodesicLine.

Parameters

[in]	lineX	segment X.
[in]	lineY	segment Y.
[out]	segmode	an indicator equal to zero if the segments intersect (see below).
[out]	c	optional pointer to an integer coincidence indicator.

Returns: p the intersection point if the segments intersect, otherwise the intersection point closest to the midpoints of the two segments.

Warning: lineX and lineY must represent shortest geodesics, e.g., they can be created by Geodesic::InverseLine. The results are only well defined if there's a unique shortest geodesic between the endpoints of the two segments.

Note: lineX and lineY should be created with minimum capabilities Intersect::LineCaps. The methods for creating a GeodesicLine include all these capabilities by default.

See previous definition of Intersect::Segment for more information on segmode.

Definition at line 83 of file Intersect.cpp.

◆ Next() [1/2]

Intersect::Point GeographicLib::Intersect::Next	(	Math::real	latX,
		Math::real	lonX,
		Math::real	aziX,
		Math::real	aziY,
		int *	c = `nullptr`
	)		const

Find the next closest intersection point to a given intersection, specified by position and two azimuths.

Parameters

[in]	latX	latitude of starting points for geodesics X and Y (degrees).
[in]	lonX	longitude of starting points for geodesics X and Y (degrees).
[in]	aziX	azimuth at starting point for geodesic X (degrees).
[in]	aziY	azimuth at starting point for geodesic Y (degrees).
[out]	c	optional pointer to an integer coincidence indicator.

Returns: p the next closest intersection point.

The returned intersection minimizes Intersect::Dist(p) (excluding p = [0,0]).

Note: Equidistant closest intersections are surprisingly common. If this may be a problem, use Intersect::All with a sufficiently large maxdist to capture close intersections.

Definition at line 91 of file Intersect.cpp.

References GeographicLib::Geodesic::Line(), LineCaps, and Next().

Referenced by main(), and Next().

◆ Next() [2/2]

Intersect::Point GeographicLib::Intersect::Next	(	const GeodesicLine &	lineX,
		const GeodesicLine &	lineY,
		int *	c = `nullptr`
	)		const

Find the next closest intersection point to a given intersection, with each geodesic specified a GeodesicLine.

Parameters

[in]	lineX	geodesic X.
[in]	lineY	geodesic Y.
[out]	c	optional pointer to an integer coincidence indicator.

Returns: p the next closest intersection point.

Warning: lineX and lineY must both have the same starting point, i.e., the distance between [lineX.Latitude(), lineX.Longitude()] and [lineY.Latitude(), lineY.Longitude()] must be zero.

Note: lineX and lineY should be created with minimum capabilities Intersect::LineCaps. The methods for creating a GeodesicLine include all these capabilities by default.; Equidistant closest intersections are surprisingly common. If this may be a problem, use Intersect::All with a sufficiently large maxdist to capture close intersections.

Definition at line 98 of file Intersect.cpp.

◆ All() [1/4]

std::vector< Intersect::Point > GeographicLib::Intersect::All	(	Math::real	latX,
		Math::real	lonX,
		Math::real	aziX,
		Math::real	latY,
		Math::real	lonY,
		Math::real	aziY,
		Math::real	maxdist,
		std::vector< int > &	c,
		const Point &	p0 = `Point(0, 0)`
	)		const

Find all intersections within a certain distance, with each geodesic specified by position and azimuth.

Parameters

[in]	latX	latitude of starting point for geodesic X (degrees).
[in]	lonX	longitude of starting point for geodesic X (degrees).
[in]	aziX	azimuth at starting point for geodesic X (degrees).
[in]	latY	latitude of starting point for geodesic Y (degrees).
[in]	lonY	longitude of starting point for geodesic Y (degrees).
[in]	aziY	azimuth at starting point for geodesic Y (degrees).
[in]	maxdist	the maximum distance for the returned intersections (meters).
[out]	c	vector of coincidences.
[in]	p0	an optional offset for the starting points (meters), default = [0,0].

Returns: plist a vector for the intersections closest to p0.

Each intersection point satisfies Intersect::Dist(p, p0) ≤ maxdist. The vector of returned intersections is sorted on the distance from p0.

Definition at line 115 of file Intersect.cpp.

References All(), GeographicLib::Geodesic::Line(), and LineCaps.

Referenced by All(), All(), and main().

◆ All() [2/4]

std::vector< Intersect::Point > GeographicLib::Intersect::All	(	Math::real	latX,
		Math::real	lonX,
		Math::real	aziX,
		Math::real	latY,
		Math::real	lonY,
		Math::real	aziY,
		Math::real	maxdist,
		const Point &	p0 = `Point(0, 0)`
	)		const

Find all intersections within a certain distance, with each geodesic specified by position and azimuth. Don't return vector of coincidences.

Parameters

[in]	latX	latitude of starting point for geodesic X (degrees).
[in]	lonX	longitude of starting point for geodesic X (degrees).
[in]	aziX	azimuth at starting point for geodesic X (degrees).
[in]	latY	latitude of starting point for geodesic Y (degrees).
[in]	lonY	longitude of starting point for geodesic Y (degrees).
[in]	aziY	azimuth at starting point for geodesic Y (degrees).
[in]	maxdist	the maximum distance for the returned intersections (meters).
[in]	p0	an optional offset for the starting points (meters), default = [0,0].

Returns: plist a vector for the intersections closest to p0.

Each intersection point satisfies Intersect::Dist(p, p0) ≤ maxdist. The vector of returned intersections is sorted on the distance from p0.

Definition at line 106 of file Intersect.cpp.

References All(), GeographicLib::Geodesic::Line(), and LineCaps.

◆ All() [3/4]

std::vector< Intersect::Point > GeographicLib::Intersect::All	(	const GeodesicLine &	lineX,
		const GeodesicLine &	lineY,
		Math::real	maxdist,
		std::vector< int > &	c,
		const Point &	p0 = `Point(0, 0)`
	)		const

Find all intersections within a certain distance, with each geodesic specified by a GeodesicLine.

Parameters

[in]	lineX	geodesic X.
[in]	lineY	geodesic Y.
[in]	maxdist	the maximum distance for the returned intersections (meters).
[out]	c	vector of coincidences.
[in]	p0	an optional offset for the starting points (meters), default = [0,0].

Returns: plist a vector for the intersections closest to p0.

Each intersection point satisfies Intersect::Dist(p, p0) ≤ maxdist. The vector of returned intersections is sorted on the distance from p0.

Note: lineX and lineY should be created with minimum capabilities Intersect::LineCaps. The methods for creating a GeodesicLine include all these capabilities by default.

Definition at line 132 of file Intersect.cpp.

◆ All() [4/4]

std::vector< Intersect::Point > GeographicLib::Intersect::All	(	const GeodesicLine &	lineX,
		const GeodesicLine &	lineY,
		Math::real	maxdist,
		const Point &	p0 = `Point(0, 0)`
	)		const

Find all intersections within a certain distance, with each geodesic specified by a GeodesicLine. Don't return vector or coincidences.

Parameters

[in]	lineX	geodesic X.
[in]	lineY	geodesic Y.
[in]	maxdist	the maximum distance for the returned intersections (meters).
[in]	p0	an optional offset for the starting points (meters), default = [0,0].

Returns: plist a vector for the intersections closest to p0.

Each intersection point satisfies Intersect::Dist(p, p0) ≤ maxdist. The vector of returned intersections is sorted on the distance from p0.

Note: lineX and lineY should be created with minimum capabilities Intersect::LineCaps. The methods for creating a GeodesicLine include all these capabilities by default.

Definition at line 125 of file Intersect.cpp.

◆ NumInverse()

long long GeographicLib::Intersect::NumInverse ( ) const

inline

Returns: the cumulative number of invocations of h.

This is a count of the number of times the spherical triangle needs to be solved. Each involves a call to Geodesic::Inverse and this is a good metric for the overall cost. This counter is set to zero by the constructor.

Warning: The counter is a mutable variable and so is not thread safe.

Definition at line 504 of file Intersect.hpp.

◆ NumBasic()

long long GeographicLib::Intersect::NumBasic ( ) const

inline

Returns: the cumulative number of invocations of b.

This is a count of the number of invocations of the basic algorithm, which is used by all the intersection methods. This counter is set to zero by the constructor.

Warning: The counter is a mutable variable and so is not thread safe.

Definition at line 514 of file Intersect.hpp.

◆ NumChange()

long long GeographicLib::Intersect::NumChange ( ) const

inline

Returns: the number of times intersection point was changed in Intersect::Closest and Intersect::Next.

If this counter is incremented by just 1 in Intersect::Closest, then the initial result of the basic algorithm was eventually accepted. This counter is set to zero by the constructor.

Note: This counter is also incremented by Intersect::Segment, which calls Intersect::Closest.

Warning: The counter is a mutable variable and so is not thread safe.

Definition at line 528 of file Intersect.hpp.

◆ NumCorner()

long long GeographicLib::Intersect::NumCorner ( ) const

inline

Returns: the number of times a corner point is checked in Intersect::Segment.

This counter is set to zero by the constructor.

Warning: The counter is a mutable variable and so is not thread safe.

Definition at line 537 of file Intersect.hpp.

◆ NumOverride()

long long GeographicLib::Intersect::NumOverride ( ) const

inline

Returns: the number of times a corner point is returned by Intersect::Segment.

This counter is set to zero by the constructor.

Note: A conjecture is that a corner point never results in an intersection that overrides the intersection closest to the midpoints of the segments; i.e., NumCorner() always returns 0.

Warning: The counter is a mutable variable and so is not thread safe.

Definition at line 550 of file Intersect.hpp.

◆ GeodesicObject()

const Geodesic & GeographicLib::Intersect::GeodesicObject ( ) const

inline

Returns: geod the Geodesic object used in the constructor.

This can be used to query Geodesic::EquatorialRadius(), Geodesic::Flattening(), Geodesic::Exact(), and Geodesic::EllipsoidArea().

Definition at line 563 of file Intersect.hpp.

◆ Dist()

static Math::real GeographicLib::Intersect::Dist	(	const Point &	p,
		const Point &	p0 = `Point(0, 0)`
	)

inlinestatic

The L1 distance.

Parameters

[in]	p	the position along geodesics X and Y.
[in]	p0	[optional] the reference position, default = [0, 0].

Returns: the L1 distance of p from p0, i.e., |p_x − p0_x| + |p_y − p0_y|.

Definition at line 576 of file Intersect.hpp.

Referenced by main().

Member Data Documentation

◆ LineCaps

const unsigned GeographicLib::Intersect::LineCaps

static

Initial value:

= Geodesic::LATITUDE | Geodesic::LONGITUDE |
      Geodesic::AZIMUTH | Geodesic::REDUCEDLENGTH | Geodesic::GEODESICSCALE |
      Geodesic::DISTANCE_IN

The minimum capabilities for GeodesicLine objects which are passed to this class.

Definition at line 84 of file Intersect.hpp.

Referenced by All(), All(), Closest(), main(), Next(), and Segment().

The documentation for this class was generated from the following files:

Public Types

Public Member Functions

Static Public Member Functions

Static Public Attributes

Detailed Description

Member Typedef Documentation

◆ Point

Constructor & Destructor Documentation

◆ Intersect()

Member Function Documentation

◆ Closest() [1/2]

◆ Closest() [2/2]

◆ Segment() [1/2]

◆ Segment() [2/2]

◆ Next() [1/2]

◆ Next() [2/2]

◆ All() [1/4]

◆ All() [2/4]

◆ All() [3/4]

◆ All() [4/4]

◆ NumInverse()

◆ NumBasic()

◆ NumChange()

◆ NumCorner()

◆ NumOverride()

◆ GeodesicObject()

◆ Dist()

Member Data Documentation

◆ LineCaps