Introduce formula for Karney's direct geodesic method #486

adl1995 · 2018-05-31T15:10:30Z

This pull request introduces the direct geodesic method as proposed by Karney (2011). The following are the highlights of the changes introduced:

Implementation of the direct geodesic method.
Series expansions of integrals given in the paper. All of these were calculated using Maxima with the geod.mac script (associated with GeographicLib), except for co-efficients of C3. I noticed there was a variation in the results produced by GeographicLib and the function generated through Maxima, which caused lon2 to be inaccurate. Replacing the function with Geodesic::C3coeff() from GeographicLib removed this inaccuracy, so I have used that instead.
Utility functions in math.hpp, such as normalization, Clenshaw's summation, and Horner's evaluation. Some of these overlap with formulas/area_formulas.hpp. I will merge them in a single file as part of a separate PR.
Test cases with existing dataset and antipodal points dataset, collected from GeodTest, associated with GeographicLib. These were converted into C++ array format using this Python script. Geodesic scale (M12) is missing from the GeodTest dataset, so I generated it manually using this C++ script.

However, the tests fail when comparing geodesic scale (M12) with Boost implementation of the direct method. Also, I noticed that using the C++ header GeographicLib/Geodesic.hpp and GeodSolve CLI tool both provided different results. As an example, using these points:

lat1: 31.863784754278001, lon1: 0, azi1: 29.572313410210999, s12: 19993324.8682601

with GeographicLib, GeodSolve CLI, and Boost gave the following results:

GeographicLib version:      -0.99598315084972932
GeodSolve tool:             -0.99588800900114727
Boost version:              -0.99625199787331664

I'm not sure why this happens, because GeodSolve uses the same C++ code at the back-end. The difference might not seem much, but since only a tolerance within 1e-7 is acceptable, these tests fail.

…unction

The paper can be found at: https://arxiv.org/pdf/1109.4448.pdf This commit also introduces the evaluate_series_A1 function for evaluating the series expantion, which was generated using Maxima: http://maxima.sourceforge.net

- Add SED script for converting x to CT(x) - Improve code documentation

- Update function headers

- Fix conversion from degree to radian in sin_cos_degrees function

…ethod - Link to paper: https://arxiv.org/pdf/1109.4448.pdf

…ies expansion

…nal difference

…pansion

…'s direct method

- Add separate function for evaluating C3 from C3x coefficient

- Add minus sign for B12 evaluation

Dataset is collected from: https://zenodo.org/record/32156 It is then parsed using a Python script.

The geodesic length is calculated manually using GeographicLib/Geodesic.hpp in C++. However, this value differs when calculated using the CLI tool GeodSolve.

vissarion · 2018-06-04T11:20:07Z

include/boost/geometry/formulas/karney_direct.hpp

+*/
+template <
+    typename CT,
+    size_t SeriesOrder = 8,


maybe place it at the end of parameters' list to be more consistent with similar strategies

Thanks! Just updated this.

vissarion · 2018-06-04T11:25:20Z

include/boost/geometry/formulas/karney_direct.hpp

+        CT const lat1 = la1;
+        Azi const azi12 = math::normalize_angle<CT>(azimuth12);
+
+        if (math::equals(distance, Dist(0)) || distance < Dist(0))


why not creating a const for 0?

Thanks! Just updated this.

vissarion · 2018-06-04T11:47:33Z

include/boost/geometry/formulas/karney_direct.hpp

+                                                math::normalize_angle(lon12));
+        }
+
+        if (BOOST_GEOMETRY_CONDITION(CalcQuantities))


this part of code duplicates code from https://github.com/boostorg/geometry/blob/develop/include/boost/geometry/formulas/differential_quantities.hpp is it possible to integrate ?

I actually tried using that code. However, almost all the test cases fail if I do. I find that odd, since it seems to be working well with Vincenty and Thomas. I will try and dig a little deeper into the code and see if I'm missing out on a step.

vissarion · 2018-06-04T12:09:24Z

include/boost/geometry/util/series_expansion.hpp

+
+namespace boost { namespace geometry { namespace series_expansion {
+
+    /*


there is a lot of duplication in the maxima scripts inside comments, could you try to decrease it?

I just removed the duplicated code in this commit. Please let me know if I should further reduce the comments.

It could be more useful to have all maxima script in one place. Because for example if someone want to reproduce series for I2 then the definition of ataylor is missing.

Also the following script (taken from area formula, right?) is repeated 7 times!

To replace each number x by CT(x) the following script can be used: sed -e 's/[0-9]\+/CT(&)/g; s/\[CT/\[/g; s/)\]/\]/g; s/case\sCT(/case /g; s/):/:/g; s/epsCT(2)/eps2/g;'

Finally, is the script taken from http://geographiclib.sourceforge.net/html/geod.mac without modifications? If so you can just add the link in the comments.

Also the following script (taken from area formula, right?) is repeated 7 times!
...

Yes, it matches the one from area formula, except the last part s/epsCT(2)/eps2/g;. Before that the script was converting eps2 -> epsCT(2). So, I will just keep a single copy of this script and remove the duplicated ones.

Finally, is the script taken from http://geographiclib.sourceforge.net/html/geod.mac without modifications? If so you can just add the link in the comments.

The scripts are indeed taken from (http://geographiclib.sourceforge.net/html/geod.mac), however, I had to make minor adjustments e.g. change the function header and data types to keep it consistent with Boost Geometry guidelines. Is it okay to discard those changes? If so, I will only provide a link to geod.mac, otherwise, I will add them in a separate file.

The scripts are indeed taken from (http://geographiclib.sourceforge.net/html/geod.mac), however, I had to make minor adjustments e.g. change the function header and data types to keep it consistent with Boost Geometry guidelines. Is it okay to discard those changes? If so, I will only provide a link to geod.mac, otherwise, I will add them in a separate file.

Ok, since they contain changes, I think it is useful to keep them in Boost.Geometry.

I have moved the Maxima scripts to: doc/other/maxima/geod.mac.

vissarion · 2018-06-04T12:21:30Z

Thanks for this PR!

Series expansions of integrals given in the paper. All of these were calculated using Maxima with the geod.mac script (associated with GeographicLib), except for co-efficients of C3. I noticed there was a variation in the results produced by GeographicLib and the function generated through Maxima, which caused lon2 to be inaccurate. Replacing the function with Geodesic::C3coeff() from GeographicLib removed this inaccuracy, so I have used that instead.

lon2 is inaccurate with respect to the data set in https://zenodo.org/record/32156 ? How do you know that Geodesic::C3coeff() is more accurate than the function generated from Maxima?

Test cases with existing dataset and antipodal points dataset, collected from GeodTest, associated with GeographicLib. These were converted into C++ array format using this Python script. Geodesic scale (M12) is missing from the GeodTest dataset, so I generated it manually using this C++ script.

Could you please add this information in the code as a comment for future reference?

In some cases, this allows the caller to ignore the CT template argument, as it is deduced from the argument list.

Other changes include the update of series expansion function calls, as the template arguments are reversed.

Previously, the array size was passed in as a separate parameter.

…n function Modified functions are: - evaluate_coeffs_C3x - normalize_values

awulkiew · 2018-06-26T22:39:33Z

include/boost/geometry/formulas/karney_direct.hpp

+
+        math::normalize_values<CT>(sin_beta1, cos_beta1);
+
+        cos_beta1 = std::max(math::sqrt(c0), cos_beta1);


std::max won't compile with msvc. Take it in brackets. Furthermore sqrt(0) is 0. So:

(std::max)(c0, cos_beta1);

or am I missing something?

awulkiew · 2018-06-26T22:41:05Z

include/boost/geometry/formulas/karney_direct.hpp

+        CT const lat1 = la1;
+
+        Azi azi12 = azimuth12;
+        math::normalize_angle<degree, Azi>(azi12);


So if you called it normalize_azimuth() then it'd be clear that it has to be in range (-180, 180].

awulkiew · 2018-06-26T22:47:37Z

include/boost/geometry/util/math.hpp

+    // Convert to radians.
+    remainder *= d2r<T>();
+
+    T s = std::sin(remainder), c = std::cos(remainder);


Before you used math functions without std::. If you specify std:: then the compiler will not use functions from the namespace where a user-defined numeric type is defined. So please remove std::. If you like you can bring std:: math functions into scope:

using std::floor; using std::sin; using std::cos;

but they are not used like that in the whole library so I'd say this is not needed.

awulkiew · 2018-06-26T22:48:01Z

include/boost/geometry/util/math.hpp

+    // the argument to the range [-45, 45] before converting it to radians.
+    T remainder; int quotient;
+
+    remainder = std::fmod(x, T(360));


Use math::mod()

awulkiew · 2018-06-26T22:48:15Z

include/boost/geometry/util/math.hpp

+        return 0;
+    }
+
+    T y = std::abs(x);


math::abs()

awulkiew · 2018-06-26T22:51:52Z

include/boost/geometry/util/series_expansion.hpp

+    template <typename Coeffs, typename CT>
+    inline void evaluate_coeffs_C3x(Coeffs &c, size_t const& SeriesOrder, CT const& n) {
+        size_t const coeff_size = Coeffs::static_size;
+        BOOST_GEOMETRY_ASSERT(coeff_size == (SeriesOrder * (SeriesOrder - 1)) / 2);


If you took SeriesOrder as template parameter you could check this at compile time:

static size_t const coeff_size = Coeffs::static_size; static size_t const expected_size = (SeriesOrder * (SeriesOrder - 1)) / 2; BOOST_MPL_ASSERT((coeff_size == expected_size));

I'm having some trouble passing SeriesOrder as the template parameter.

I changed the evaluate_coeffs_C3x function header to:

template <typename Coeffs, typename CT, size_t SeriesOrder> inline void evaluate_coeffs_C3x(Coeffs &c, CT const& n) {

and called it as:

evaluate_coeffs_C3x(*this, n);

but I get this error:

error: ‘evaluate_coeffs_C3x’ was not declared in this scope

I also tried changing the order of parameters as:

template <size_t SeriesOrder, typename CT, typename Coeffs> inline void evaluate_coeffs_C3x(CT const& n, Coeffs &c) {

but calling the function as:

evaluate_coeffs_C3x(n, *this);

or:

evaluate_coeffs_C3x<SeriesOrder, CT>(n, *this);

doesn't work either.

It produces the error:

error: expected primary-expression before ‘>’ token

and:

couldn't deduce template parameter ‘SeriesOrder’

respectively.

If you're calling these functions inside structs' constructors then they should be defined before structs.

Either of these should work:

template <size_t SeriesOrder, typename CT, typename Coeffs> inline void evaluate_coeffs_C3x(CT const& n, Coeffs &c){} // this is more strict, it takes only the struct you like template <size_t SeriesOrder, typename CT> inline void evaluate_coeffs_C3x(CT const& n, coeffs_C3x<SeriesOrder, CT> &c){}

and either of these should work:

evaluate_coeffs_C3x<SeriesOrder>(n, *this); evaluate_coeffs_C3x<SeriesOrder, CT>(n, *this); evaluate_coeffs_C3x<SeriesOrder, CT, coeffs_C3x>(n, *this); // assuming it's called in coeffs_C3x's constructor

Thanks, this worked.

Out of curiosity, why was this working before for the previous function calls? Does this have something to do with passing explicit template parameters? As I only got the error with evaluate_coeffs_C3x<SeriesOrder>.

No idea. You'd have to know how exactly type deduction is done in your compiler. Are you using Visual Studio? Its compiler allows things it shouldn't sometimes.

No, I'm using g++ version 7.2.0 on Ubuntu.

awulkiew · 2018-06-26T22:54:22Z

include/boost/geometry/util/series_expansion.hpp

+            int m = Coeffs1::static_size - l - 1;
+            mult *= eps;
+
+            std::vector<CT> coeffs2_slice(coeffs2.begin(), coeffs2.begin() + offset);


You don't have to create a container here, one each iteration btw. This will be horribly slow. Pass range into math::polyval().

awulkiew · 2018-06-26T22:56:54Z

include/boost/geometry/util/math.hpp

+\brief Evaluate the polynomial.
+*/
+template<typename CT>
+inline CT polyval(std::vector<CT> coeff,


template <typename It> inline CT polyval(It first, It last, std::iterator_traits<It>::value_type const& eps)

or just

template <typename It, typename T> inline CT polyval(It first, It last, T const& eps)

awulkiew · 2018-06-26T22:57:27Z

include/boost/geometry/util/math.hpp

+inline CT polyval(std::vector<CT> coeff,
+                  const CT eps)
+{
+    int N = boost::size(coeff) - 1;


std::distance(first, last) - 1

awulkiew · 2018-06-26T22:58:16Z

include/boost/geometry/util/math.hpp

+    int N = boost::size(coeff) - 1;
+    int index = 0;
+
+    CT y = N < 0 ? 0 : coeff[index++];


*(first + (index++))

For normalizing longitudes, the normalize_longitude function is used instead.

…) function The coefficient container structs are moved to the bottom of the file.

This is done to avoid creating a separate container in each iteration.

Conflicts: include/boost/geometry/util/math.hpp test/formulas/direct.cpp The conflicting files have been updated.

adl1995 · 2018-06-29T09:57:25Z

Thank you for the detailed reviews. I think all the requested changes have been addressed. I have also merged develop into this branch, and resolved all conflicting files.

Please let me know if there are any further changes required.

awulkiew · 2018-07-02T13:22:53Z

Thanks! I'm ok with merging. I'd prefer to wait until after the 1.68 release. Or do you need this PR merged in order to work on the next step of the project?

adl1995 · 2018-07-02T14:51:55Z

That's fine. I can work on local branches for the next part i.e. the inverse method.

Btw, when is the next release scheduled for?

awulkiew · 2018-07-03T12:31:47Z

Here is the schedule: https://www.boost.org/development/index.html
The branches for final release are closed 1 Aug however after beta is released we only do bugfixing.

I think the whole month is too long to keep you wait so I'm ok with merging it now. If needed I'll create a branch for bugfixing based on previous develop state.

adl1995 · 2018-07-03T14:05:41Z

Okay, as you deem fit. If it is merged before August 14th, I can mention this in my final evaluation report.

adl1995 added 22 commits May 12, 2018 22:08

[utils] Add function sin_cos_degrees for evaluating sine and cosine f…

302e240

…unction

[formulas] Add function for evaluating coefficients for C1

c73ea92

- Add SED script for converting x to CT(x) - Improve code documentation

[utils] Add function for normalizing and rounding off an angle

320891c

[formulas] Compute sin cos series using Clenshaw summation

b6fcee8

- Update function headers

[util] Move series expansion functions to util/series_expansion.hpp

3700b4f

[util] Evaluate coefficients for C1p using series expansion

6ce1b09

- Fix conversion from degree to radian in sin_cos_degrees function

[formulas] Compute the latitude for second point following Karney's m…

881c73c

…ethod - Link to paper: https://arxiv.org/pdf/1109.4448.pdf

[util] Add functions to evaluate coefficients for A3 and C3 using ser…

9152191

…ies expansion

[formulas] Compute the longitude for second point using the longitudi…

9e2c74a

…nal difference

[util] Evaluate series and coefficients for A2 and C2 using series ex…

a0a55f3

…pansion

[formulas] Compute the reduced length and geodesic scale using Karney…

ac4e483

…'s direct method

[formulas] Make variable declarations constant

0c2b8cd

[formulas][util] Add missing import and function return type

6a2897d

[test] Add direct Karney's method to direct.cpp test cases

afe7bc6

[util] Modify function for evaluting C3x coefficient

1e31876

- Add separate function for evaluating C3 from C3x coefficient

[util] Add functions for normalizing and evaluating polynomial

3dd6bce

[formulas] Fix direct geodesic method by performing normalization

b8a225e

- Add minus sign for B12 evaluation

[test] Add nearly antipodal points dataset for direct geodesic problem

fdbb388

Dataset is collected from: https://zenodo.org/record/32156 It is then parsed using a Python script.

[test] Add geodesic length to antipodal points dataset

3c21b13

The geodesic length is calculated manually using GeographicLib/Geodesic.hpp in C++. However, this value differs when calculated using the CLI tool GeodSolve.

[test] Test Karney's method on antipodal points dataset

6ff97a6

[util] Update series expansion for C3x

4fa4a82

adl1995 mentioned this pull request Jun 1, 2018

Progress update BoostGSoC18/geometry#1

Open

vissarion reviewed Jun 4, 2018

View reviewed changes

adl1995 added 2 commits June 4, 2018 19:09

[formulas] Move SeriesOrder to the end of template parameter list

1ed5f10

[formulas] Use constant type variables for comparison

8318737

adl1995 added 5 commits June 25, 2018 18:39

[util] Reverse template argument order for series expansion functions

19f06fa

In some cases, this allows the caller to ignore the CT template argument, as it is deduced from the argument list.

[formulas] Use assignment operator on the same line for consistency

06eb057

Other changes include the update of series expansion function calls, as the template arguments are reversed.

[util] Avoid passing array size using std::vector

9c96bec

Previously, the array size was passed in as a separate parameter.

[util] Add BOOST_GEOMETRY_ASSERT in series expansion and normalizatio…

4d39eae

…n function Modified functions are: - evaluate_coeffs_C3x - normalize_values

[util] Change static inline to inline in series_expansion.hpp file

c26483b

awulkiew reviewed Jun 26, 2018

View reviewed changes

include/boost/geometry/util/math.hpp Outdated

return 0;

}

T y = std::abs(x);

Copy link

Member

awulkiew Jun 26, 2018

Choose a reason for hiding this comment

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

math::abs()

awulkiew reviewed Jun 26, 2018

View reviewed changes

adl1995 added 5 commits June 27, 2018 10:34

[util][formulas] Rename normalize_angle function to normalize_azimuth

9699b49

For normalizing longitudes, the normalize_longitude function is used instead.

[util] Use functions from math namespace instead of std

9ff4fbe

[util] Rename math::normalize_values to math::normalize_unit_vector

b86a93b

[util] Pass SeriesOrder as template parameter in evaluate_coeffs_C3x(…

1fe3b3b

…) function The coefficient container structs are moved to the bottom of the file.

[util] Pass range into math::polyval() instead of std::vector

6219503

This is done to avoid creating a separate container in each iteration.

adl1995 force-pushed the feature/geodesic_direct branch from ab45a02 to 6219503 Compare June 28, 2018 06:00

Merge branch 'develop' into feature/geodesic_direct

687df8e

Conflicts: include/boost/geometry/util/math.hpp test/formulas/direct.cpp The conflicting files have been updated.

awulkiew approved these changes Jul 2, 2018

View reviewed changes

awulkiew added this to the 1.69 milestone Jul 25, 2018

awulkiew merged commit 79ef70f into boostorg:develop Jul 25, 2018

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Introduce formula for Karney's direct geodesic method #486

Introduce formula for Karney's direct geodesic method #486


		namespace boost { namespace geometry { namespace series_expansion {

		/*


		math::normalize_values<CT>(sin_beta1, cos_beta1);

		cos_beta1 = std::max(math::sqrt(c0), cos_beta1);

Introduce formula for Karney's direct geodesic method #486

Introduce formula for Karney's direct geodesic method #486

Conversation

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment