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 direct formulas #489

Merged
merged 17 commits into from Jun 20, 2018

Conversation

Projects
None yet
2 participants
@vissarion
Contributor

vissarion commented Jun 12, 2018

This PR is a split from #431 and contains 3 direct formulas, namely for the spherical case, the thomas first order formula for geographic CS (the direct equivalent for andoyer first order inverse formula) and a direct formula for the meridian geographic case.

@vissarion vissarion requested a review from awulkiew Jun 12, 2018

@vissarion vissarion referenced this pull request Jun 12, 2018

Open

Progress update #1

{
CT const f = formula::flattening<CT>(spheroid);
CT n = f / (CT(2) - f);
CT mp = 10001965.729;

This comment has been minimized.

@awulkiew

awulkiew Jun 12, 2018

Member

Quarter meridian should be calculated using spheroid, here WGS84 is always used.

// https://en.wikipedia.org/wiki/Meridian_arc#The_inverse_meridian_problem_for_the_ellipsoid
// latitudes are assumed to be in radians and in [-pi/2,pi/2]
template <typename T, typename Spheroid>
static CT apply(T m, Spheroid const& spheroid)

This comment has been minimized.

@awulkiew

awulkiew Jun 12, 2018

Member

What is the input and result of this formula? Why does it have different interface than the other direct formulas? AFAIU direct formula should take latitude, direction (north or south) and distance. The result should be the coordinates of the resulting point since they may be on the other side of the globe. What am I missing here?

This comment has been minimized.

@vissarion

vissarion Jun 18, 2018

Contributor

That was just the formula for computation. I added the interface.

@@ -217,6 +219,47 @@ inline int azimuth_side_value(T const& azi_a1_p, T const& azi_a1_a2)
: 1; // left
}
template <typename CT, typename Sphere>
inline result_direct<CT> spherical_direct(CT const& lon1,

This comment has been minimized.

@awulkiew

awulkiew Jun 12, 2018

Member

An alternative formula would convert to 3D cartesian ECEF and calculate the result using vector algebra. This is how it's done in the spherical intersection and side strategies. So should be careful to not mix things. But AFAIU this formula is a building block for other geographic formulas and not used in spherical for now. Is my understanding correct?

This comment has been minimized.

@vissarion

vissarion Jun 18, 2018

Contributor

Yes this is not used anywhere for now, I added it for completeness. I could be used in geographic pt-seg distance formula for optimizations.

result.lon2 = lon1 + lon2 - omg1;
result.lat2 = lat2;
result.reverse_azimuth = alp2;

This comment has been minimized.

@awulkiew

awulkiew Jun 12, 2018

Member

Are all of these results needed in all cases? In other formulas they are only calculated if compile-time flag is set. The same approach could be used here.

@@ -59,7 +59,8 @@ class thomas_direct
T const& la1,
Dist const& distance,
Azi const& azimuth12,
Spheroid const& spheroid)
Spheroid const& spheroid,
bool SecondOrder = true)

This comment has been minimized.

@awulkiew

awulkiew Jun 12, 2018

Member

Shouldn't this be template parameter?

Dist const& distance,
Azi const& azimuth12,
Spheroid const& spheroid)
{ return thomas_direct<

This comment has been minimized.

@awulkiew

awulkiew Jun 12, 2018

Member

I propose to implement thomas as one struct taking int Order or bool FirstOrder as the first or second template parameter, e.g. simply as thomas_direct. And then derive both thomas 1st and 2nd order from that if needed. You can implement all of them in one file. If the order flag was struct template parameter you could simply derive from it without the need to implement the interface again.

This comment has been minimized.

@vissarion

vissarion Jun 18, 2018

Contributor

I followed a common interface solution since some formulas rely on it e.g. https://github.com/boostorg/geometry/blob/develop/include/boost/geometry/formulas/gnomonic_intersection.hpp#L39

Now changed it to (the more convenient, I agree) interface with SecondOrder as template parameter and use strategies (in particular formula policies) for gnomonic_intersection in tests https://github.com/boostorg/geometry/pull/489/files#diff-faacdab3bb307b84e30e9bde5ca88591R74 However, this should be fixed in the future. Maybe by just moving formulas policies into formulas?

@@ -22,7 +22,7 @@
#include <boost/geometry/core/radius.hpp>
#include <boost/geometry/formulas/andoyer_inverse.hpp>
#include <boost/geometry/formulas/elliptic_arc_length.hpp>
#include <boost/geometry/formulas/elliptic_meridian_arc_inverse.hpp>

This comment has been minimized.

@awulkiew

awulkiew Jun 12, 2018

Member

Btw, for consistency, shouldn't it be called "spheroidal meridian arc"? Are we ok with using "spheroid" and "ellipsoid" interchangeably? It would probably only be a problem if we wanted to support 3-axial ellipsoids.

This comment has been minimized.

@vissarion

vissarion Jun 18, 2018

Contributor

ok for now I simply use the meridian_direct and meridian_inverse

IteratorType it = end;
do
{
result = result * x + *--it;

This comment has been minimized.

@awulkiew

awulkiew Jun 13, 2018

Member

What if begin == end? Either return 0 or check assertion.

ResultType const ab4 = (get_radius<0>(geometry)
+ get_radius<2>(geometry)) / c4;
return geometry::math::pi<ResultType>() * ab4 *
horner_evaluate(n*n, C, C+8) / C[0];

This comment has been minimized.

@awulkiew

awulkiew Jun 13, 2018

Member

So this function will calculate 8 nested levels of products of products or products, etc. My guess is that the error propagation will be big. Is this a problem? Would representing this series in a "typical" way improve the accuracy of the result, i.e. as sum of ratios multiplied by x^i?

This comment has been minimized.

@vissarion

vissarion Jun 14, 2018

Contributor

First, using horner rule (i.e. nested evaluation) should be faster than the standard evaluation, the first requires 2n operations while the second 3n-1 (if powers are computed iteratively). Second, horner seems to me more accurate since it uses less operations. Third, if this function is placed somewhere in geometry::math or in utilities then it could be used wherever we have to evaluate some polynomial. This will make the code shorter and more elegant than implementing manually polynomial evaluations in standard form when it is needed.

@vissarion

This comment has been minimized.

Contributor

vissarion commented Jun 18, 2018

@awulkiew Thanks for the comments. I addressed them all by commits or comments.

@awulkiew

This comment has been minimized.

Member

awulkiew commented Jun 18, 2018

Thanks! I'm ok with merging.

@@ -329,7 +329,7 @@ class area_formulas
std::size_t begin(0), end(0);
for(std::size_t i = 0; i <= SeriesOrder; i++){
end = begin + SeriesOrder + 1 - i;
coeffs2[i] = ((i==0) ? CT(1) : pow(var2,int(i)))
coeffs2[i] = ((i==0) ? CT(1) : pow(var2, CT(i)))

This comment has been minimized.

@awulkiew

awulkiew Jun 18, 2018

Member
test/algorithms/area/area.cpp:216:22:   required from here
../../boost/geometry/formulas/area_formulas.hpp:332:59: error: no matching function for call to ‘pow(double&, double)’
             coeffs2[i] = ((i==0) ? CT(1) : pow(var2, CT(i)))
                                                           ^

This comment has been minimized.

@awulkiew

awulkiew Jun 18, 2018

Member

In various places in area_formulas.hpp variables are passed by value instead of const&. May this be the cause?

This comment has been minimized.

@vissarion

vissarion Jun 19, 2018

Contributor

I use pass by reference but I get the same error in circleci. Note that this PR does not touch area_formulas.hpp i.e. the error also appear in 12f7a22 It seems that the error is in calling boost pow function from boost/rational.hpp although I cannot find pow function there and I cannot reproduce the error in my pc with gcc-4.9 and debug mode.

This comment has been minimized.

@vissarion

vissarion Jun 20, 2018

Contributor

@awulkiew tried with the std pow function in area and now pow functions in projections yield errors see: https://circleci.com/gh/vissarion/geometry/249#tests/containers/1

Did you have similar problems in other branches?

This comment has been minimized.

@vissarion

vissarion Jun 20, 2018

Contributor

Ok there is a general issue #491 with pow function independent from this PR.

This comment has been minimized.

@awulkiew

awulkiew Jun 20, 2018

Member

@vissarion right, so this should be fixed library-wise at our side. AFAIK the correct way of fixing it is bringing std namespace into the scope and then calling pow() so:

using namespace std;
pow(1, 2);

This way the compiler should choose the correct one. And actually we could have our own function doing exactly that (in util/math.hpp):

namespace geometry { namespace math {

template <typename T>
T pow(T const& v1, T const& v2)
{
    using namespace std;
    return pow(v1, v2);
}

}}

and then call:

math::pow(v1, v2);

This comment has been minimized.

@awulkiew

awulkiew Jun 20, 2018

Member

If it didn't work the function could be brought into the scope. This is how it's done in Boost.Math (https://github.com/boostorg/math/blob/develop/include/boost/math/tools/config.hpp#L308):

using std::pow;
pow(1, 2);

This is probably better since we're not enabling the whole namespace.

@awulkiew

This comment has been minimized.

Member

awulkiew commented Jun 20, 2018

I'm ok with merging it now and fixing pow() afterwards. Since the problem is mainly in projections I'll fix it as described above, by implementing geometry::math::pow().

@awulkiew awulkiew merged commit 53ff12a into boostorg:develop Jun 20, 2018

1 check failed

ci/circleci Your tests failed on CircleCI
Details

@vissarion vissarion deleted the vissarion:feature_direct_formulas branch Jun 20, 2018

@awulkiew awulkiew added this to the 1.68 milestone Jun 20, 2018

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment