{"action":"created","issue":{"url":"https://api.github.com/repos/elki-project/elki/issues/29","repository_url":"https://api.github.com/repos/elki-project/elki","labels_url":"https://api.github.com/repos/elki-project/elki/issues/29/labels{/name}","comments_url":"https://api.github.com/repos/elki-project/elki/issues/29/comments","events_url":"https://api.github.com/repos/elki-project/elki/issues/29/events","html_url":"https://github.com/elki-project/elki/issues/29","id":220482284,"number":29,"title":"Geodetic mindist","user":{"login":"brano543","id":6738774,"avatar_url":"https://avatars0.githubusercontent.com/u/6738774?v=3","gravatar_id":"","url":"https://api.github.com/users/brano543","html_url":"https://github.com/brano543","followers_url":"https://api.github.com/users/brano543/followers","following_url":"https://api.github.com/users/brano543/following{/other_user}","gists_url":"https://api.github.com/users/brano543/gists{/gist_id}","starred_url":"https://api.github.com/users/brano543/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/brano543/subscriptions","organizations_url":"https://api.github.com/users/brano543/orgs","repos_url":"https://api.github.com/users/brano543/repos","events_url":"https://api.github.com/users/brano543/events{/privacy}","received_events_url":"https://api.github.com/users/brano543/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"assignees":[],"milestone":null,"comments":4,"created_at":"2017-04-09T16:35:54Z","updated_at":"2017-04-10T13:19:53Z","closed_at":null,"body":"Hello,\r\n\r\ni was reading your paper \"Geodetic Distance Queries on R-Trees for Indexing Geographic Data\" and i wasn´t sure what did you mean by all those subscript 360 functions, so i checked your code for that defined here https://github.com/elki-project/elki/blob/e8f3c6fdf54e1e0aa8444b94ad5374ad518dfc0c/elki-core-math/src/main/java/de/lmu/ifi/dbs/elki/math/geodesy/SphereUtil.java, it doesn´t follow the pseudocode. For instance lines 806-813\r\n```\r\n// Determine whether going east or west is shorter.\r\n    double lngE = rminlng - plng;\r\n    if(lngE < 0) {\r\n      lngE += MathUtil.TWOPI;\r\n    }\r\n    double lngW = rmaxlng - plng; // we keep this negative!\r\n    if(lngW > 0) {\r\n      lngW -= MathUtil.TWOPI;\r\n}\r\n```\r\nwhere i think you meant to do\r\n```\r\n//mod360(rminlng-plng) <= mod360(plng-rmaxlng)\r\n// Determine whether going east or west is shorter.\r\n        double lngW = fmod(rminlng - plng, 2 * M_PI);\r\n        if (lngW < 0)\r\n        {\r\n            lngW += 2 * M_PI;\r\n        }\r\n        double lngE = fmod(plng - rmaxlng, 2 * M_PI);\r\n        if (lngE < 0)\r\n        {\r\n            lngE += 2 * M_PI;\r\n}\r\n```\r\nAlso my theory is that 360 subscript means we are working in radians, not degrees. In fact when i followed your pseudocode in paper exactly it is working for me in C++, but i am still not sure what you meant by sentence \"In order to distinguish the other cases, we first need to test whether we are on the left or on the right side by rotating the mean longitude of the rectangle by 180◦ – the meridian opposite of the rectangle.\" . Do i understand correctly, that you wanted to do modulo 360(2 PI) to normalize the angle to range 0-360 [0..2 PI]?\r\n\r\nAlso i have found out,that the haversine formula should use absolute difference of latitudes and longitudes. Correct me if i am wrong,but if you dont use absolute values there,you may get negative distances. I attach my version of the algorithm including the Haversine formula\r\n\r\n```\r\nstatic const double wgs84_radius_m = 6378137;\r\nstatic const double wgs84_flattening = 1.0 / 298.257223563;\r\nstatic const double earth_radius_m = wgs84_radius_m * (1 - wgs84_flattening);\r\nstatic const double earth_radius_km = earth_radius_m / 1000.0;\r\n struct WGS84Mindist\r\n{\r\n    //input is in radians\r\n    //returns distance in kilometers\r\n    static double haversineFormulaRad(double lat1, double lon1, double lat2, double lon2)\r\n    {\r\n        double d_lat = abs(lat1 - lat2);\r\n        double d_lon = abs(lon1 - lon2);\r\n\r\n        double a = pow(sin(d_lat / 2), 2) + cos(lat1) * cos(lat2) * pow(sin(d_lon / 2), 2);\r\n\r\n        //double d_sigma = 2 * atan2(sqrt(a), sqrt(1 - a));\r\n        double d_sigma = 2 * asin(sqrt(a));\r\n\r\n        return earth_radius_km * d_sigma;\r\n    }\r\n\r\n    //input is in radians\r\n    //returns angular distance on unit sphere\r\n    static double haversineFormulaRadAngular(double lat1, double lon1, double lat2, double lon2)\r\n    {\r\n        double d_lat = abs(lat1 - lat2);\r\n        double d_lon = abs(lon1 - lon2);\r\n\r\n        double a = pow(sin(d_lat / 2), 2) + cos(lat1) * cos(lat2) * pow(sin(d_lon / 2), 2);\r\n\r\n        //double d_sigma = 2 * atan2(sqrt(a), sqrt(1 - a));\r\n        double d_sigma = 2 * asin(sqrt(a));\r\n\r\n        return d_sigma;\r\n    }\r\n\r\n    //input is in radians\r\n    //returns angle in radians\r\n    static double getBearingRad(double lat1, double lon1, double lat2, double lon2)\r\n    {\r\n        double dLon = lon2 - lon1;\r\n        double y = sin(dLon) * cos(lat2);\r\n        double x = cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(dLon);\r\n        double radiansBearing = atan2(y, x);\r\n\r\n        return radiansBearing;\r\n    }\r\n\r\n    //start, end, point - coordinates in radians\r\n    static double getCrossTrackDistanceRad(double lat1, double lon1, double lat2, double lon2, double lat3, double lon3)\r\n    {\r\n        double angDist1Q = haversineFormulaRadAngular(lat1, lon1, lat3, lon3);\r\n\r\n        double cos_lat1 = cos(lat1);\r\n        double sin_lat1 = sin(lat1);\r\n        double cos_lat3 = cos(lat3);\r\n        double cos_lat2 = cos(lat2);\r\n\r\n        //double theta1Q = getBearing(lat1, lon1, lat3, lon3);\r\n        // double dLon = lon3 - lon1;\r\n        // double y = sin(dLon) * cos(lat3);\r\n        // double x = cos(lat1) * sin(lat3) - sin(lat1) * cos(lat3) * cos(dLon);\r\n        // double radiansBearing = atan2(y, x);\r\n        double dLon1 = lon3 - lon1;\r\n        double y1 = sin(dLon1) * cos_lat3;\r\n        double x1 = cos_lat1 * sin(lat3) - sin_lat1 * cos_lat3 * cos(dLon1);\r\n        double theta1Q = atan2(y1, x1);\r\n\r\n        //double theta12 = getBearing(lat1, lon1, lat2, lon2);\r\n        // double dLon = lon2 - lon1;\r\n        // double y = sin(dLon) * cos(lat2);\r\n        // double x = cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(dLon);\r\n        // double radiansBearing = atan2(y, x);\r\n        double dLon2 = lon2 - lon1;\r\n        double y2 = sin(dLon2) * cos_lat2;\r\n        double x2 = cos_lat1 * sin(lat2) - sin_lat1 * cos_lat2 * cos(dLon2);\r\n        double theta12 = atan2(y2, x2);\r\n\r\n        return asin(sin(angDist1Q) * sin(theta1Q - theta12)) * earth_radius_km;\r\n    }\r\n\r\n    static double latlngMinDistDeg(double &plat, double &plng, double &rminlat, double &rminlng, double &rmaxlat, double &rmaxlng)\r\n    {\r\n        return latlngMinDistRad(GeoDistance::deg2rad(plat), GeoDistance::deg2rad(plng),\r\n                                GeoDistance::deg2rad(rminlat), GeoDistance::deg2rad(rminlng),\r\n                                GeoDistance::deg2rad(rmaxlat), GeoDistance::deg2rad(rmaxlng));\r\n    }\r\n\r\n    //returns distance in kilometers\r\n    static double latlngMinDistRad(double plat, double plng, double rminlat, double rminlng, double rmaxlat, double rmaxlng)\r\n    {\r\n        // The simplest case is when the query point is in the same \"slice\":\r\n        if (rminlng <= plng && plng <= rmaxlng)\r\n        {\r\n            if (plat < rminlat)\r\n            {\r\n                return (rminlat - plat) * earth_radius_km; //Sout of MBR\r\n            }\r\n            else if (plat > rmaxlat)\r\n            {\r\n                return (plat - rmaxlat) * earth_radius_km; //North of MBR\r\n            }\r\n            return 0.0; // INSIDE MBR\r\n        }\r\n\r\n        // Determine whether going east or west is shorter.\r\n        double lngW = fmod(rminlng - plng, 2 * M_PI);\r\n        if (lngW < 0)\r\n        {\r\n            lngW += 2 * M_PI;\r\n        }\r\n        double lngE = fmod(plng - rmaxlng, 2 * M_PI);\r\n        if (lngE < 0)\r\n        {\r\n            lngE += 2 * M_PI;\r\n        }\r\n\r\n        if (lngW <= lngE)\r\n        {\r\n            // West of MBR\r\n            double tau = tan(plat);\r\n\r\n            if (lngW >= 0.5 * M_PI) // Large delta of longitude\r\n            {\r\n                if (tau <= tan((rmaxlat + rminlat) * 0.5) * cos(rminlng - plng))\r\n                {\r\n                    return haversineFormulaRad(plat, plng, rmaxlat, rminlng); //North-West\r\n                }\r\n                else\r\n                {\r\n                    return haversineFormulaRad(plat, plng, rminlat, rminlng); //South-West\r\n                }\r\n            }\r\n\r\n            if (tau >= tan(rmaxlat) * cos(rminlng - plng))\r\n            {\r\n                return haversineFormulaRad(plat, plng, rmaxlat, rminlng); //North-West\r\n            }\r\n\r\n            if (tau <= tan(rminlat) * cos(rminlng - plng))\r\n            {\r\n                return haversineFormulaRad(plat, plng, rminlat, rminlng); //South-West\r\n            }\r\n\r\n            return abs(getCrossTrackDistanceRad(rminlat, rminlng, rmaxlat, rminlng, plat, plng)); // West\r\n        }     \r\n        else\r\n        {\r\n             // East of MBR\r\n             double tau = tan(plat);\r\n\r\n             if (lngE >= 0.5 * M_PI) // Large delta of longitude\r\n             {\r\n                 if (tau <= tan((rmaxlat + rminlat) * 0.5) * cos(rmaxlng - plng))\r\n                 {\r\n                    return haversineFormulaRad(plat, plng, rmaxlat, rmaxlng); //North-East\r\n                 }\r\n                 else\r\n                 {\r\n                    return haversineFormulaRad(plat, plng, rminlat, rmaxlng); //Sout-East\r\n                 }\r\n             }\r\n\r\n             if (tau >= tan(rmaxlat) * cos(rmaxlng - plng))\r\n             {\r\n                 return haversineFormulaRad(plat, plng, rmaxlat, rmaxlng); //North-East\r\n             }\r\n\r\n             if (tau <= tan(rminlat) * cos(rmaxlng - plng))\r\n             {\r\n                 return haversineFormulaRad(plat, plng, rminlat, rmaxlng); //Sout-East\r\n             }\r\n\r\n             return abs(getCrossTrackDistanceRad(rmaxlat, rmaxlng, rminlat, rmaxlng, plat, plng)); // East\r\n        }\r\n    }\r\n};\r\n```\r\nThe last thing i dont understand why in pseudocode you are returning c * (rminlat-plat)/360 for the N/S case when i think it should be c * (rminlat-plat) , because you say in text we are calculating length of meridian arc, where c is radius of earth.\r\n\r\nThank you for your answer."},"comment":{"url":"https://api.github.com/repos/elki-project/elki/issues/comments/292946895","html_url":"https://github.com/elki-project/elki/issues/29#issuecomment-292946895","issue_url":"https://api.github.com/repos/elki-project/elki/issues/29","id":292946895,"user":{"login":"brano543","id":6738774,"avatar_url":"https://avatars0.githubusercontent.com/u/6738774?v=3","gravatar_id":"","url":"https://api.github.com/users/brano543","html_url":"https://github.com/brano543","followers_url":"https://api.github.com/users/brano543/followers","following_url":"https://api.github.com/users/brano543/following{/other_user}","gists_url":"https://api.github.com/users/brano543/gists{/gist_id}","starred_url":"https://api.github.com/users/brano543/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/brano543/subscriptions","organizations_url":"https://api.github.com/users/brano543/orgs","repos_url":"https://api.github.com/users/brano543/repos","events_url":"https://api.github.com/users/brano543/events{/privacy}","received_events_url":"https://api.github.com/users/brano543/received_events","type":"User","site_admin":false},"created_at":"2017-04-10T13:19:53Z","updated_at":"2017-04-10T13:19:53Z","body":"So okay, you say the angles after mod360 should be in range -pi,+pi, what do you think about the solution presented here http://stackoverflow.com/questions/11498169/dealing-with-angle-wrap-in-\r\n```\r\nNormalize to [-180,180):\r\n\r\ndouble constrainAngle(double x){\r\n    x = fmod(x + 180,360);\r\n    if (x < 0)\r\n        x += 360;\r\n    return x - 180;\r\n}\r\n```"}}
{"actor":{"display_login":"brano543"}}