feat(vincenty): full convergence
This commit is contained in:
parent
7841a9ef2e
commit
2a5052f47f
4 changed files with 267 additions and 65 deletions
30
.idea/shelf/Uncommitted_changes_before_Update_at_3_24_25,_12_40 PM_[Changes]/shelved.patch
generated
Normal file
30
.idea/shelf/Uncommitted_changes_before_Update_at_3_24_25,_12_40 PM_[Changes]/shelved.patch
generated
Normal file
File diff suppressed because one or more lines are too long
7
.idea/shelf/Uncommitted_changes_before_Update_at_3_24_25__12_40PM__Changes_.xml
generated
Normal file
7
.idea/shelf/Uncommitted_changes_before_Update_at_3_24_25__12_40PM__Changes_.xml
generated
Normal file
|
|
@ -0,0 +1,7 @@
|
|||
<changelist name="Uncommitted_changes_before_Update_at_3_24_25,_12_40 PM_[Changes]" date="1742834400705" recycled="false" toDelete="true">
|
||||
<option name="PATH" value="$PROJECT_DIR$/.idea/shelf/Uncommitted_changes_before_Update_at_3_24_25,_12_40 PM_[Changes]/shelved.patch" />
|
||||
<option name="DESCRIPTION" value="Uncommitted changes before Update at 3/24/25, 12:40 PM [Changes]" />
|
||||
<binary>
|
||||
<option name="BEFORE_PATH" value="crates/client/dist/wxbox-client-505f3eac107f11ba_bg.wasm" />
|
||||
</binary>
|
||||
</changelist>
|
||||
41
.idea/workspace.xml
generated
41
.idea/workspace.xml
generated
|
|
@ -8,13 +8,8 @@
|
|||
</component>
|
||||
<component name="ChangeListManager">
|
||||
<list default="true" id="2d855648-9644-469a-afa2-59beb52bb1d6" name="Changes" comment="feat: client work">
|
||||
<change afterPath="$PROJECT_DIR$/client/src/lib/map/fragment.glsl" afterDir="false" />
|
||||
<change afterPath="$PROJECT_DIR$/client/src/lib/map/vertex.glsl" afterDir="false" />
|
||||
<change beforePath="$PROJECT_DIR$/.idea/.gitignore" beforeDir="false" />
|
||||
<change beforePath="$PROJECT_DIR$/.idea/vcs.xml" beforeDir="false" afterPath="$PROJECT_DIR$/.idea/vcs.xml" afterDir="false" />
|
||||
<change beforePath="$PROJECT_DIR$/.idea/wxbox.iml" beforeDir="false" afterPath="$PROJECT_DIR$/.idea/wxbox.iml" afterDir="false" />
|
||||
<change beforePath="$PROJECT_DIR$/client/src/lib/Map.svelte" beforeDir="false" afterPath="$PROJECT_DIR$/client/src/lib/Map.svelte" afterDir="false" />
|
||||
<change beforePath="$PROJECT_DIR$/crates/tiler/config.toml" beforeDir="false" afterPath="$PROJECT_DIR$/crates/tiler/config.toml" afterDir="false" />
|
||||
<change beforePath="$PROJECT_DIR$/.idea/workspace.xml" beforeDir="false" afterPath="$PROJECT_DIR$/.idea/workspace.xml" afterDir="false" />
|
||||
<change beforePath="$PROJECT_DIR$/client/src/lib/map/fragment.glsl" beforeDir="false" afterPath="$PROJECT_DIR$/client/src/lib/map/fragment.glsl" afterDir="false" />
|
||||
</list>
|
||||
<option name="SHOW_DIALOG" value="false" />
|
||||
<option name="HIGHLIGHT_CONFLICTS" value="true" />
|
||||
|
|
@ -39,22 +34,23 @@
|
|||
<option name="hideEmptyMiddlePackages" value="true" />
|
||||
<option name="showLibraryContents" value="true" />
|
||||
</component>
|
||||
<component name="PropertiesComponent"><![CDATA[{
|
||||
"keyToString": {
|
||||
"RunOnceActivity.ShowReadmeOnStart": "true",
|
||||
"RunOnceActivity.rust.reset.selective.auto.import": "true",
|
||||
"git-widget-placeholder": "core/client-rewrite-again",
|
||||
"node.js.detected.package.eslint": "true",
|
||||
"node.js.selected.package.eslint": "(autodetect)",
|
||||
"node.js.selected.package.tslint": "(autodetect)",
|
||||
"nodejs_package_manager_path": "npm",
|
||||
"org.rust.cargo.project.model.PROJECT_DISCOVERY": "true",
|
||||
"org.rust.first.attach.projects": "true",
|
||||
"settings.editor.selected.configurable": "preferences.pluginManager",
|
||||
"ts.external.directory.path": "/home/core/prj/personal/wxbox/client/node_modules/typescript/lib",
|
||||
"vue.rearranger.settings.migration": "true"
|
||||
<component name="PropertiesComponent">{
|
||||
"keyToString": {
|
||||
"RunOnceActivity.ShowReadmeOnStart": "true",
|
||||
"RunOnceActivity.git.unshallow": "true",
|
||||
"RunOnceActivity.rust.reset.selective.auto.import": "true",
|
||||
"git-widget-placeholder": "core/client-rewrite-again",
|
||||
"node.js.detected.package.eslint": "true",
|
||||
"node.js.selected.package.eslint": "(autodetect)",
|
||||
"node.js.selected.package.tslint": "(autodetect)",
|
||||
"nodejs_package_manager_path": "npm",
|
||||
"org.rust.cargo.project.model.PROJECT_DISCOVERY": "true",
|
||||
"org.rust.first.attach.projects": "true",
|
||||
"settings.editor.selected.configurable": "preferences.pluginManager",
|
||||
"ts.external.directory.path": "/home/core/prj/personal/wxbox/client/node_modules/typescript/lib",
|
||||
"vue.rearranger.settings.migration": "true"
|
||||
}
|
||||
}]]></component>
|
||||
}</component>
|
||||
<component name="RunManager">
|
||||
<configuration name="Run wxbox-ar2" type="CargoCommandRunConfiguration" factoryName="Cargo Command">
|
||||
<option name="command" value="run --package wxbox-ar2 --bin wxbox-ar2" />
|
||||
|
|
@ -130,6 +126,7 @@
|
|||
<workItem from="1748031993681" duration="8000" />
|
||||
<workItem from="1748032062973" duration="74000" />
|
||||
<workItem from="1748032143420" duration="1307000" />
|
||||
<workItem from="1748132890039" duration="1214000" />
|
||||
</task>
|
||||
<servers />
|
||||
</component>
|
||||
|
|
|
|||
|
|
@ -2,7 +2,11 @@
|
|||
|
||||
precision highp float;
|
||||
|
||||
#define PI 3.141592653589793238462643
|
||||
const float PI = 3.141592653589793238462643;
|
||||
const float PI_OVER_TWO = PI / 2.0;
|
||||
const float TWO_PI = PI + PI;
|
||||
const float NEGATIVE_PI_OVER_TWO = -PI_OVER_TWO;
|
||||
const float NEGATIVE_TWO_PI = -TWO_PI;
|
||||
|
||||
out highp vec4 fragColor;
|
||||
uniform mat4 u_matrix;
|
||||
|
|
@ -15,57 +19,213 @@ void xyToLngLat(in float x, in float y, out float lat, out float lng) {
|
|||
|
||||
in vec4 raw_pos;
|
||||
|
||||
#define a 6378137.0
|
||||
#define f 1.0 / 298.257223563
|
||||
#define b 6356752.314245
|
||||
struct Ellipsoid {
|
||||
float semiMajorAxisMeters;
|
||||
float semiMinorAxisMeters;
|
||||
float flattening;
|
||||
float inverseFlattening;
|
||||
};
|
||||
Ellipsoid fromAAndInverseF(float semiMajorAxisMeters, float inverseFlattening) {
|
||||
float f = 1.0 / inverseFlattening;
|
||||
float b = (1.0 - f) * semiMajorAxisMeters;
|
||||
|
||||
bool vincenty(in float phi1, in float phi2, in float l1, in float l2, out float a1, out float a2, out float s) {
|
||||
float u1 = atan((1.0 - f) * tan(phi1));
|
||||
float u2 = atan((1.0 - f) * tan(phi2));
|
||||
float L = l2 - l1;
|
||||
return Ellipsoid(semiMajorAxisMeters, b, f, inverseFlattening);
|
||||
}
|
||||
|
||||
float lambda = L;
|
||||
struct Angle {
|
||||
float radians;
|
||||
};
|
||||
float angleAsDegrees(Angle angle) {
|
||||
return degrees(angle.radians);
|
||||
}
|
||||
Angle radianAngle(float radians) {
|
||||
return Angle(radians);
|
||||
}
|
||||
Angle degreeAngle(float degrees) {
|
||||
return Angle(radians(degrees));
|
||||
}
|
||||
|
||||
float sinSigma;
|
||||
float cosSigma;
|
||||
float sigma;
|
||||
float sinAlpha;
|
||||
float cosSquaredAlpha;
|
||||
float cos2OmegaM;
|
||||
float C;
|
||||
struct GlobalCoordinates {
|
||||
Angle latitude;
|
||||
Angle longitude;
|
||||
};
|
||||
void canonicalizeGlobalCoordinates(inout GlobalCoordinates coords) {
|
||||
float latitudeRadians = coords.latitude.radians;
|
||||
float longitudeRadians = coords.longitude.radians;
|
||||
|
||||
// sometimes fails to converge. set maximum iteration count
|
||||
int iterCount = 0;
|
||||
latitudeRadians = mod((latitudeRadians + PI), TWO_PI);
|
||||
if (latitudeRadians < 0.0) {
|
||||
latitudeRadians += TWO_PI;
|
||||
}
|
||||
latitudeRadians -= PI;
|
||||
|
||||
while (lambda > pow(10.0, -12.0)) {
|
||||
iterCount++;
|
||||
if (iterCount > 10) {
|
||||
return false;
|
||||
}
|
||||
sinSigma = sqrt(pow(cos(u2)*sin(lambda), 2.0) + pow(cos(u1)*sin(u2)-sin(u1)*cos(u2)*cos(lambda),2.0));
|
||||
cosSigma = sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(lambda);
|
||||
sigma = atan(sinSigma,cosSigma);
|
||||
|
||||
sinAlpha = (cos(u1)*cos(u2)*sin(lambda))/(sinSigma);
|
||||
cosSquaredAlpha = 1.0 - pow(sinAlpha, 2.0);
|
||||
cos2OmegaM = cosSigma - (2.0*sin(u1)*sin(u2))/(cosSquaredAlpha);
|
||||
C = f/16.0 * cosSquaredAlpha * (4.0 + f * (4.0 - 3.0 * cosSquaredAlpha));
|
||||
lambda = L + (1.0 - C) * f * sinAlpha * (sigma + C * sinSigma * (cos2OmegaM + C*cosSigma * (-1.0 + 2.0 * pow(cos2OmegaM, 2.0))));
|
||||
if (latitudeRadians > PI_OVER_TWO) {
|
||||
latitudeRadians = PI - latitudeRadians;
|
||||
longitudeRadians += PI;
|
||||
} else if (latitudeRadians < NEGATIVE_PI_OVER_TWO) {
|
||||
latitudeRadians = -PI - latitudeRadians;
|
||||
longitudeRadians += PI;
|
||||
}
|
||||
|
||||
float uSquared = cosSquaredAlpha * ((pow(a, 2.0) - pow(b, 2.0))/pow(b, 2.0));
|
||||
float A = 1.0 + uSquared/16384.0 * (4096.0 + uSquared * (-768.0 + uSquared * (320.0 - 175.0*uSquared)));
|
||||
float B = uSquared / 1024.0 * (256.0 + uSquared * (-128.0 + uSquared * (74.0 - 47.0*uSquared)));
|
||||
longitudeRadians = mod((longitudeRadians + PI), TWO_PI);
|
||||
if (longitudeRadians <= 0.0) {
|
||||
longitudeRadians += TWO_PI;
|
||||
}
|
||||
longitudeRadians -= PI;
|
||||
|
||||
float deltaSigma = B * sinSigma * (cos2OmegaM + 1.0/4.0*B * (cosSigma * (-1.0 + 2.0*pow(cos2OmegaM,2.0)) - 1.0/6.0*B*cos2OmegaM*(-3.0+4.0*pow(sinSigma,2.0))*(-3.0+4.0*pow(cos2OmegaM,2.0))));
|
||||
s = b*A*(sigma - deltaSigma);
|
||||
a1 = atan(cos(u2)*sin(lambda), cos(u1)*sin(u2)-sin(u1)*cos(u2)*cos(lambda));
|
||||
a2 = atan(cos(u1)*sin(lambda), -sin(u1)*cos(u2)+cos(u1)*sin(u2)*cos(lambda));
|
||||
coords.latitude = radianAngle(latitudeRadians);
|
||||
coords.longitude = radianAngle(longitudeRadians);
|
||||
}
|
||||
|
||||
return true;
|
||||
struct Geodedic {
|
||||
float s;
|
||||
Angle a1;
|
||||
Angle a2;
|
||||
bool didGetGoodEstimate;
|
||||
};
|
||||
|
||||
const float tolerance = pow(10.0, -13.0);
|
||||
|
||||
Geodedic vincenty(Ellipsoid ellipsoid, GlobalCoordinates start, GlobalCoordinates end) {
|
||||
float a = ellipsoid.semiMajorAxisMeters;
|
||||
float b = ellipsoid.semiMinorAxisMeters;
|
||||
float f = ellipsoid.flattening;
|
||||
|
||||
float phi1 = start.latitude.radians;
|
||||
float lambda1 = start.longitude.radians;
|
||||
float phi2 = end.latitude.radians;
|
||||
float lambda2 = end.longitude.radians;
|
||||
|
||||
float a2 = a * a;
|
||||
float b2 = b * b;
|
||||
float a2b2b2 = (a2 - b2) / b2;
|
||||
|
||||
float omega = lambda2 - lambda1;
|
||||
|
||||
float tanphi1 = tan(phi1);
|
||||
float tanU1 = (1.0 - f) * tanphi1;
|
||||
float U1 = atan(tanU1);
|
||||
float sinU1 = sin(U1);
|
||||
float cosU1 = cos(U1);
|
||||
|
||||
float tanphi2 = tan(phi2);
|
||||
float tanU2 = (1.0 - f) * tanphi2;
|
||||
float U2 = atan(tanU2);
|
||||
float sinU2 = sin(U2);
|
||||
float cosU2 = cos(U2);
|
||||
|
||||
float sinU1sinU2 = sinU1 * sinU2;
|
||||
float cosU1sinU2 = cosU1 * sinU2;
|
||||
float sinU1cosU2 = sinU1 * cosU2;
|
||||
float cosU1cosU2 = cosU1 * cosU2;
|
||||
|
||||
// equ. 13
|
||||
float lambda = omega;
|
||||
|
||||
// intermediates to compute 's'
|
||||
float A = 0.0;
|
||||
float B = 0.0;
|
||||
float sigma = 0.0;
|
||||
float deltasigma = 0.0;
|
||||
float lambda0;
|
||||
bool converged = false;
|
||||
|
||||
for (int i = 0; i < 20; i++) {
|
||||
lambda0 = lambda;
|
||||
float sinlambda = sin(lambda);
|
||||
float coslambda = cos(lambda);
|
||||
// equ. 14
|
||||
float cosU1sinU2_sinU2cosU2coslambda = cosU1sinU2 - sinU1cosU2 * coslambda;
|
||||
float sin2sigma = (cosU2 * sinlambda * cosU2 * sinlambda) + (cosU1sinU2_sinU2cosU2coslambda * cosU1sinU2_sinU2cosU2coslambda);
|
||||
float sinsigma = sqrt(sin2sigma);
|
||||
|
||||
// equ. 15
|
||||
float cossigma = sinU1sinU2 + (cosU1cosU2 * coslambda);
|
||||
|
||||
// equ. 16
|
||||
sigma = atan(sinsigma, cossigma);
|
||||
|
||||
// equ. 17 - careful, sin2sigma might be almost 0
|
||||
float sinalpha = (sin2sigma == 0.0) ? 0.0 : cosU1cosU2 * sinlambda / sinsigma;
|
||||
float alpha = asin(sinalpha);
|
||||
float cosalpha = cos(alpha);
|
||||
float cos2alpha = cosalpha * cosalpha;
|
||||
|
||||
// equ. 18 - careful, cos2alpha might be almost 0
|
||||
float cos2sigmam = (cos2alpha == 0.0) ? 0.0 : cossigma - 2.0 * sinU1sinU2 / cos2alpha;
|
||||
float u2 = cos2alpha * a2b2b2;
|
||||
float cos2sigmam2 = cos2sigmam * cos2sigmam;
|
||||
|
||||
// equ. 3
|
||||
A = 1.0 + u2 / 16384.0 * (4096.0 + u2 * (-768.0 + u2 * (320.0 - 175.0 * u2)));
|
||||
// equ. 4
|
||||
B = u2 / 1024.0 * (256.0 + u2 * (-128.0 + u2 * (74.0 - 47.0 * u2)));
|
||||
// equ. 6
|
||||
deltasigma = B * sinsigma * (cos2sigmam + B / 4.0 * (cossigma * (-1.0 + 2.0 * cos2sigmam2) - B / 6.0 * cos2sigmam * (-3.0 + 4.0 * sin2sigma) * (-3.0 + 4.0 * cos2sigmam2)));
|
||||
// equ. 10
|
||||
float C = f / 16.0 * cos2alpha * (4.0 + f * (4.0 - 3.0 * cos2alpha));
|
||||
// equ. 11 (modified)
|
||||
lambda = omega + (1.0 - C) * f * sinalpha * (sigma + C * sinsigma * (cos2sigmam + C * cossigma * (-1.0 + 2.0 * cos2sigmam2)));
|
||||
|
||||
if (i < 2) {
|
||||
continue;
|
||||
}
|
||||
|
||||
float change = abs((lambda - lambda0) / lambda);
|
||||
if (change < tolerance) {
|
||||
converged = true;
|
||||
break;
|
||||
}
|
||||
|
||||
// equ. 19
|
||||
float s = b * A * (sigma - deltasigma);
|
||||
Angle alpha1;
|
||||
Angle alpha2;
|
||||
|
||||
bool didGetGoodEstimate = true;
|
||||
|
||||
// didn't converge? must be N/S
|
||||
if (!converged) {
|
||||
if (phi1 > phi2) {
|
||||
alpha1 = degreeAngle(180.0);
|
||||
alpha2 = degreeAngle(0.0);
|
||||
} else if (phi1 < phi2) {
|
||||
alpha1 = degreeAngle(0.0);
|
||||
alpha2 = degreeAngle(180.0);
|
||||
} else {
|
||||
alpha1 = Angle(0.0 / 0.0); // NaN
|
||||
alpha2 = Angle(0.0 / 0.0); // NaN
|
||||
didGetGoodEstimate = false;
|
||||
}
|
||||
} else {
|
||||
float radians;
|
||||
// equ. 20
|
||||
radians = atan(cosU2 * sin(lambda), (cosU1sinU2 - sinU1cosU2 * cos(lambda)));
|
||||
if (radians < 0.0) {
|
||||
radians += TWO_PI;
|
||||
}
|
||||
alpha1 = radianAngle(radians);
|
||||
|
||||
radians = atan(cosU1 * sin(lambda), (-sinU1cosU2 + cosU1sinU2 * cos(lambda))) + PI;
|
||||
if (radians < 0.0) {
|
||||
radians += TWO_PI;
|
||||
}
|
||||
alpha2 = radianAngle(radians);
|
||||
}
|
||||
|
||||
if (alpha1.radians >= TWO_PI) {
|
||||
alpha1 = radianAngle(alpha1.radians - TWO_PI);
|
||||
}
|
||||
if (alpha2.radians >= TWO_PI) {
|
||||
alpha2 = radianAngle(alpha2.radians - TWO_PI);
|
||||
}
|
||||
|
||||
return Geodedic(s, alpha1, alpha2, didGetGoodEstimate);
|
||||
}
|
||||
}
|
||||
|
||||
void main() {
|
||||
Ellipsoid WGS84 = fromAAndInverseF(6378137.0, 298.257223563);
|
||||
float lat;
|
||||
float lng;
|
||||
xyToLngLat(raw_pos.x, raw_pos.y, lat, lng);
|
||||
|
|
@ -76,17 +236,25 @@ void main() {
|
|||
float a1;
|
||||
float a2;
|
||||
float s;
|
||||
;
|
||||
|
||||
// should be 43.58777778
|
||||
// 88 > x > 89
|
||||
// >40.0
|
||||
float delta = 0.01;
|
||||
|
||||
GlobalCoordinates radar = GlobalCoordinates(degreeAngle(radarLat), degreeAngle(radarLng));
|
||||
canonicalizeGlobalCoordinates(radar);
|
||||
GlobalCoordinates
|
||||
samplePoint = GlobalCoordinates(degreeAngle(lat), degreeAngle(lng));
|
||||
canonicalizeGlobalCoordinates(samplePoint);
|
||||
|
||||
Geodedic vincentyResult = vincenty(WGS84, radar, samplePoint);
|
||||
|
||||
//
|
||||
// abs(lat-43.5877) <delta && abs(lng+96.72944444)<delta
|
||||
if (vincenty(radians(radarLat), radians(lat), radians(radarLng), radians(lng), a1, a2, s)) {
|
||||
if (vincentyResult.didGetGoodEstimate) {
|
||||
fragColor = vec4(0.0, 1.0, 0.0, 0.5);
|
||||
} else {
|
||||
fragColor = vec4(1.0, 0.0, 0.0, 0.5);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue