diff --git a/client/src/lib/Map.svelte b/client/src/lib/Map.svelte index 0334f4a..a485edc 100644 --- a/client/src/lib/Map.svelte +++ b/client/src/lib/Map.svelte @@ -231,12 +231,8 @@ gl.useProgram(this.program); - console.log(drd.radials.length); gl.uniform1i(gl.getUniformLocation(this.program, 'azimuthCount'), drd.radials.length); - console.log('done'); - console.log( - gl.getUniform(this.program, gl.getUniformLocation(this.program, 'azimuthCount')) - ); + gl.uniform1fv( gl.getUniformLocation(this.program, 'azimuthAngles'), new Float32Array( @@ -272,49 +268,31 @@ } } } - const dataTexture = gl.createTexture(); - gl.activeTexture(gl.TEXTURE0); - gl.bindTexture(gl.TEXTURE_2D, dataTexture); - const level = 0; - const internalFormat = gl.R32F; - const width = 720; - const height = 1832; - const border = 0; - const format = gl.RED; - const type = gl.FLOAT; - // interleave the data - const data: number[] = []; - for (let y = 0; y < height; y++) { - for (let x = 0; x < width; x++) { - const radial = drd.radials[x]; - if (radial.product && radial.product.data) { - data.push(scaleMomentData(radial, radial.product.data.data[y])); - } + const rdata: number[] = []; + for (const radial of drd.radials) { + for (const tdata of radial.product.data.data) { + rdata.push(scaleMomentData(radial, tdata)); } } - const typedData = new Float32Array(data); - gl.texImage2D( - gl.TEXTURE_2D, - level, - internalFormat, - width, - height, - border, - format, - type, - data - ); - gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.NEAREST); - gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.NEAREST); + const data = new Float32Array(rdata); + + gl.uniform1i(gl.getUniformLocation(this.program, 'scaledData'), 4); + + gl.activeTexture(gl.TEXTURE0 + 4); + this.texture = gl.createTexture(); + gl.bindTexture(gl.TEXTURE_2D, this.texture); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_S, gl.CLAMP_TO_EDGE); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_T, gl.CLAMP_TO_EDGE); - const alignment = 1; - gl.pixelStorei(gl.UNPACK_ALIGNMENT, alignment); - - gl.uniform1i(gl.getUniformLocation(this.program, 'scaledData'), 0); + gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.NEAREST); + gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.NEAREST); + gl.texImage2D(gl.TEXTURE_2D, 0, gl.R32F, 720, 1832, 0, gl.RED, gl.FLOAT, data); }, render(gl, args) { gl.useProgram(this.program); + + gl.activeTexture(gl.TEXTURE0 + 4); + gl.bindTexture(gl.TEXTURE_2D, this.texture); + gl.uniformMatrix4fv( gl.getUniformLocation(this.program, 'u_matrix'), false, diff --git a/client/src/lib/map/fragment.glsl b/client/src/lib/map/fragment.glsl index 0fecb24..64b64c4 100644 --- a/client/src/lib/map/fragment.glsl +++ b/client/src/lib/map/fragment.glsl @@ -19,211 +19,6 @@ void xyToLngLat(in float x, in float y, out float lat, out float lng) { in vec4 raw_pos; -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; - - return Ellipsoid(semiMajorAxisMeters, b, f, inverseFlattening); -} - -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)); -} - -struct GlobalCoordinates { - Angle latitude; - Angle longitude; -}; -void canonicalizeGlobalCoordinates(inout GlobalCoordinates coords) { - float latitudeRadians = coords.latitude.radians; - float longitudeRadians = coords.longitude.radians; - - latitudeRadians = mod((latitudeRadians + PI), TWO_PI); - if (latitudeRadians < 0.0) { - latitudeRadians += TWO_PI; - } - latitudeRadians -= PI; - - if (latitudeRadians > PI_OVER_TWO) { - latitudeRadians = PI - latitudeRadians; - longitudeRadians += PI; - } else if (latitudeRadians < NEGATIVE_PI_OVER_TWO) { - latitudeRadians = -PI - latitudeRadians; - longitudeRadians += PI; - } - - longitudeRadians = mod((longitudeRadians + PI), TWO_PI); - if (longitudeRadians <= 0.0) { - longitudeRadians += TWO_PI; - } - longitudeRadians -= PI; - - coords.latitude = radianAngle(latitudeRadians); - coords.longitude = radianAngle(longitudeRadians); -} - -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); - } -} - uniform int azimuthCount; uniform float[720] azimuthAngles; uniform float azimuthSpacing; @@ -266,7 +61,6 @@ LocateRadialResult locateRadial(float forAzimuth) { } void main() { - Ellipsoid WGS84 = fromAAndInverseF(6378137.0, 298.257223563); float lat; float lng; xyToLngLat(raw_pos.x, raw_pos.y, lat, lng); @@ -276,24 +70,24 @@ void main() { return; } - GlobalCoordinates radar = GlobalCoordinates(degreeAngle(radarLat), degreeAngle(radarLng)); - canonicalizeGlobalCoordinates(radar); - GlobalCoordinates - samplePoint = GlobalCoordinates(degreeAngle(lat), degreeAngle(lng)); - canonicalizeGlobalCoordinates(samplePoint); + float R = 6371.0 * pow(10.0, 3.0); // meters + float phi1 = radians(radarLat); + float phi2 = radians(lat); + float lambda1 = radians(radarLng); + float lambda2 = radians(lng); - Geodedic vincentyResult = vincenty(WGS84, radar, samplePoint); + float deltaPhi = radians(lat - radarLat); + float deltaLambda = radians(lng - radarLng); + float a = sin(deltaPhi / 2.0) * sin(deltaPhi / 2.0) + cos(phi1) * cos(phi2) * sin(deltaLambda / 2.0) * sin(deltaLambda / 2.0); + float c = 2.0 * atan(sqrt(a), sqrt(1.0 - a)); + float d = R * c; // meters - if (!vincentyResult.didGetGoodEstimate) { - fragColor = vec4(0.0, 0.0, 0.0, 0.0); - return; - } + float d_m = d; - float d_m = vincentyResult.s; - float azimuth = degrees(vincentyResult.a1.radians); - if (azimuth < 0.0) { - azimuth += 360.0; - } + float y = sin(lambda2 - lambda1) * cos(phi2); + float x = cos(phi1) * sin(phi2) - sin(phi1) * cos(phi2) * cos(lambda2 - lambda1); + float theta = atan(y, x); + float azimuth = (theta * 180.0 / PI + 360.0); // degrees LocateRadialResult maybeRadial = locateRadial(azimuth); if (!maybeRadial.didFindRadial) { @@ -306,7 +100,7 @@ void main() { float gate_spacing_km = sampleInterval / 1000.0; if (distance_km < first_gate_distance_km) { - fragColor = vec4(1.0, 0.0, 0.0, 1.0); + fragColor = vec4(0.0, 0.0, 0.0, 0.0); return; } @@ -316,8 +110,27 @@ void main() { return; } - vec2 coords = vec2(float(maybeRadial.radialIndex), float(gate)); - float rawValue = texture(scaledData, coords).r; + if (!maybeRadial.didFindRadial) { + fragColor = vec4(0.0, 0.0, 0.0, 0.0); + return; + } + + ivec2 coords = ivec2(maybeRadial.radialIndex, gate); + float rawValue = texelFetch(scaledData, coords, 0).r; + + const float BELOW_THRESHOLD = -9999.0; + const float RANGE_FOLDED = -9998.0; + + if (rawValue == BELOW_THRESHOLD) { + fragColor = vec4(0.0, 0.0, 0.0, 0.0); + return; + } + if (rawValue == RANGE_FOLDED) { + fragColor = vec4(1.0, 0.0, 1.0, 1.0); + return; + } + + //float rawValue = (float(value.r) - float(value.g)) / float(value.b); if (rawValue > 80.0) { fragColor = vec4(128.0 / 255.0, 128.0 / 255.0, 128.0 / 255.0, 1.0);