Skip to content

Commit

Permalink
Bug#30461595 IN SOME CASES QUERIES WITH ST_CONTAINS DO NOT RETURN ANY…
Browse files Browse the repository at this point in the history 
… RESULTS

This was caused by a bug in boost::geometry, causing too large Minimum
Bounding Rectangles to be calculated. Fixed by a patch from our boost
geometry developers.

Change-Id: Ia0b6d7f993e31c75e382ce0818d650ed94ddd037
  • Loading branch information
@torje
torje committed Mar 6, 2020
1 parent 48a00ad commit d2395fe
Show file tree
Hide file tree
Showing 3 changed files with 382 additions and 0 deletions.
293 changes: 293 additions & 0 deletions include/boost_1_70_0/patches/boost/geometry/formulas/andoyer_inverse.hpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,293 @@
// Boost.Geometry

// Copyright (c) 2018 Adam Wulkiewicz, Lodz, Poland.

// Copyright (c) 2015-2020 Oracle and/or its affiliates.

// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle

// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_GEOMETRY_FORMULAS_ANDOYER_INVERSE_HPP
#define BOOST_GEOMETRY_FORMULAS_ANDOYER_INVERSE_HPP


#include <boost/math/constants/constants.hpp>

#include <boost/geometry/core/radius.hpp>

#include <boost/geometry/util/condition.hpp>
#include <boost/geometry/util/math.hpp>

#include <boost/geometry/formulas/differential_quantities.hpp>
#include <boost/geometry/formulas/flattening.hpp>
#include <boost/geometry/formulas/result_inverse.hpp>


namespace boost { namespace geometry { namespace formula
{

/*!
\brief The solution of the inverse problem of geodesics on latlong coordinates,
Forsyth-Andoyer-Lambert type approximation with first order terms.
\author See
- Technical Report: PAUL D. THOMAS, MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS, 1965
http://www.dtic.mil/docs/citations/AD0627893
- Technical Report: PAUL D. THOMAS, SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY, 1970
http://www.dtic.mil/docs/citations/AD703541
*/
template <
typename CT,
bool EnableDistance,
bool EnableAzimuth,
bool EnableReverseAzimuth = false,
bool EnableReducedLength = false,
bool EnableGeodesicScale = false
>
class andoyer_inverse
{
static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale;
static const bool CalcAzimuths = EnableAzimuth || EnableReverseAzimuth || CalcQuantities;
static const bool CalcFwdAzimuth = EnableAzimuth || CalcQuantities;
static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities;

public:
typedef result_inverse<CT> result_type;

template <typename T1, typename T2, typename Spheroid>
static inline result_type apply(T1 const& lon1,
T1 const& lat1,
T2 const& lon2,
T2 const& lat2,
Spheroid const& spheroid)
{
result_type result;

// coordinates in radians

if ( math::equals(lon1, lon2) && math::equals(lat1, lat2) )
{
return result;
}

CT const c0 = CT(0);
CT const c1 = CT(1);
CT const pi = math::pi<CT>();
CT const f = formula::flattening<CT>(spheroid);

CT const dlon = lon2 - lon1;
CT const sin_dlon = sin(dlon);
CT const cos_dlon = cos(dlon);
CT const sin_lat1 = sin(lat1);
CT const cos_lat1 = cos(lat1);
CT const sin_lat2 = sin(lat2);
CT const cos_lat2 = cos(lat2);

// H,G,T = infinity if cos_d = 1 or cos_d = -1
// lat1 == +-90 && lat2 == +-90
// lat1 == lat2 && lon1 == lon2
CT cos_d = sin_lat1*sin_lat2 + cos_lat1*cos_lat2*cos_dlon;
// on some platforms cos_d may be outside valid range
if (cos_d < -c1)
cos_d = -c1;
else if (cos_d > c1)
cos_d = c1;

CT const d = acos(cos_d); // [0, pi]
CT const sin_d = sin(d); // [-1, 1]

if ( BOOST_GEOMETRY_CONDITION(EnableDistance) )
{
CT const K = math::sqr(sin_lat1-sin_lat2);
CT const L = math::sqr(sin_lat1+sin_lat2);
CT const three_sin_d = CT(3) * sin_d;

CT const one_minus_cos_d = c1 - cos_d;
CT const one_plus_cos_d = c1 + cos_d;
// cos_d = 1 means that the points are very close
// cos_d = -1 means that the points are antipodal

CT const H = math::equals(one_minus_cos_d, c0) ?
c0 :
(d + three_sin_d) / one_minus_cos_d;
CT const G = math::equals(one_plus_cos_d, c0) ?
c0 :
(d - three_sin_d) / one_plus_cos_d;

CT const dd = -(f/CT(4))*(H*K+G*L);

CT const a = CT(get_radius<0>(spheroid));

result.distance = a * (d + dd);
}

if ( BOOST_GEOMETRY_CONDITION(CalcAzimuths) )
{
// sin_d = 0 <=> antipodal points (incl. poles) or very close
if (math::equals(sin_d, c0))
{
// T = inf
// dA = inf
// azimuth = -inf

// TODO: The following azimuths are inconsistent with distance
// i.e. according to azimuths below a segment with antipodal endpoints
// travels through the north pole, however the distance returned above
// is the length of a segment traveling along the equator.
// Furthermore, this special case handling is only done in andoyer
// formula.
// The most correct way of fixing it is to handle antipodal regions
// correctly and consistently across all formulas.

// points very close
if (cos_d >= c0)
{
result.azimuth = c0;
result.reverse_azimuth = c0;
}
// antipodal points
else
{
// Set azimuth to 0 unless the first endpoint is the north pole
if (! math::equals(sin_lat1, c1))
{
result.azimuth = c0;
result.reverse_azimuth = pi;
}
else
{
result.azimuth = pi;
result.reverse_azimuth = 0;
}
}
}
else
{
CT const c2 = CT(2);

CT A = c0;
CT U = c0;
if (math::equals(cos_lat2, c0))
{
if (sin_lat2 < c0)
{
A = pi;
}
}
else
{
CT const tan_lat2 = sin_lat2/cos_lat2;
CT const M = cos_lat1*tan_lat2-sin_lat1*cos_dlon;
A = atan2(sin_dlon, M);
CT const sin_2A = sin(c2*A);
U = (f/ c2)*math::sqr(cos_lat1)*sin_2A;
}

CT B = c0;
CT V = c0;
if (math::equals(cos_lat1, c0))
{
if (sin_lat1 < c0)
{
B = pi;
}
}
else
{
CT const tan_lat1 = sin_lat1/cos_lat1;
CT const N = cos_lat2*tan_lat1-sin_lat2*cos_dlon;
B = atan2(sin_dlon, N);
CT const sin_2B = sin(c2*B);
V = (f/ c2)*math::sqr(cos_lat2)*sin_2B;
}

CT const T = d / sin_d;

// even with sin_d == 0 checked above if the second point
// is somewhere in the antipodal area T may still be great
// therefore dA and dB may be great and the resulting azimuths
// may be some more or less arbitrary angles

if (BOOST_GEOMETRY_CONDITION(CalcFwdAzimuth))
{
CT const dA = V*T - U;
result.azimuth = A - dA;
normalize_azimuth(result.azimuth, A, dA);
}

if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth))
{
CT const dB = -U*T + V;
if (B >= 0)
result.reverse_azimuth = pi - B - dB;
else
result.reverse_azimuth = -pi - B - dB;
normalize_azimuth(result.reverse_azimuth, B, dB);
}
}
}

if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
{
CT const b = CT(get_radius<2>(spheroid));

typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 1> quantities;
quantities::apply(dlon, sin_lat1, cos_lat1, sin_lat2, cos_lat2,
result.azimuth, result.reverse_azimuth,
b, f,
result.reduced_length, result.geodesic_scale);
}

return result;
}

private:
static inline void normalize_azimuth(CT & azimuth, CT const& A, CT const& dA)
{
CT const c0 = 0;

if (A >= c0) // A indicates Eastern hemisphere
{
if (dA >= c0) // A altered towards 0
{
if (azimuth < c0)
{
azimuth = c0;
}
}
else // dA < 0, A altered towards pi
{
CT const pi = math::pi<CT>();
if (azimuth > pi)
{
azimuth = pi;
}
}
}
else // A indicates Western hemisphere
{
if (dA <= c0) // A altered towards 0
{
if (azimuth > c0)
{
azimuth = c0;
}
}
else // dA > 0, A altered towards -pi
{
CT const minus_pi = -math::pi<CT>();
if (azimuth < minus_pi)
{
azimuth = minus_pi;
}
}
}
}
};

}}} // namespace boost::geometry::formula


#endif // BOOST_GEOMETRY_FORMULAS_ANDOYER_INVERSE_HPP
1 change: 1 addition & 0 deletions unittest/gunit/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -206,6 +206,7 @@ SET(SERVER_TESTS
gis_is_simple
gis_isvalid
gis_relops
gis_rtree_support
gis_srs
gis_wkb_parser
gis_wkb_writer
Expand Down
88 changes: 88 additions & 0 deletions unittest/gunit/gis_rtree_support-t.cc
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,88 @@

/*
Copyright (c) 2020, Oracle and/or its affiliates. All rights reserved.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License, version 2.0,
as published by the Free Software Foundation.
This program is also distributed with certain software (including
but not limited to OpenSSL) that is licensed under separate terms,
as designated in a particular file or component or in included license
documentation. The authors of MySQL hereby grant you an additional
permission to link the program and your derivative works with the
separately licensed software that they have included with MySQL.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License, version 2.0, for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
#include <gtest/gtest.h>
#include "sql/dd/dd.h"
#include "sql/dd/impl/types/spatial_reference_system_impl.h"
#include "sql/dd/types/spatial_reference_system.h"
#include "sql/gis/geometries.h"
#include "sql/gis/geometries_cs.h"
#include "sql/gis/mbr_utils.h"

namespace mbr_utils_unittest {
gis::Geographic_linearring linearringFromVector(std::vector<double> data) {
if (data.size() % 2 != 0) {
throw std::exception(); /* purecov: dead code */
}
gis::Geographic_linearring lr;
for (size_t i = 0; i + 1 < data.size(); i += 2) {
lr.push_back(gis::Geographic_point(data[i], data[i + 1]));
}
return lr;
}

gis::Geographic_polygon polygonFromVector(std::vector<double> data) {
gis::Geographic_polygon py;
py.push_back(linearringFromVector(data));
return py;
}

std::unique_ptr<dd::Spatial_reference_system_impl> get_srs() {
// EPSG 4326, but with long-lat axes (E-N).
std::unique_ptr<dd::Spatial_reference_system_impl> m_srs(
dynamic_cast<dd::Spatial_reference_system_impl *>(
dd::create_object<dd::Spatial_reference_system>()));
m_srs->set_id(4326);
m_srs->set_name("WGS 84");
m_srs->set_created(0UL);
m_srs->set_last_altered(0UL);
m_srs->set_organization("EPSG");
m_srs->set_organization_coordsys_id(4326);
m_srs->set_definition(
"GEOGCS[\"WGS 84\",DATUM[\"World Geodetic System "
"1984\",SPHEROID[\"WGS "
"84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],"
"AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY["
"\"EPSG\",\"8901\"]],UNIT[\"degree\",0.017453292519943278,"
"AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Lon\",EAST],AXIS[\"Lat\","
"NORTH],AUTHORITY[\"EPSG\",\"4326\"]]");
m_srs->set_description("");
m_srs->parse_definition();
return m_srs;
}

TEST(MBR_Utils_test, test) {
auto polygon = polygonFromVector({0.26377745335935482, 0.84009247535954357,
0.26387320519564444, 0.84008767903244352,
0.26387320679785603, 0.84008767099694759,
0.26377745335935482, 0.84009247535954357});
gis::Geographic_box box;
std::unique_ptr<dd::Spatial_reference_system> srs = get_srs();
gis::box_envelope(&polygon, srs.get(), &box);
EXPECT_GT(box.max_corner().x(), 0.26);
EXPECT_GT(box.max_corner().y(), 0.84);
EXPECT_LT(box.max_corner().x(), 0.27);
EXPECT_LT(box.max_corner().y(), 0.85);
}
} // namespace mbr_utils_unittest

0 comments on commit d2395fe

Please sign in to comment.