preface

Javascript API GL It is a 3D map API based on WebGL technology, with more free vision and smoother interaction. It provides a variety of functional interfaces, including point, line, surface drawing, custom layer, personalized style and drawing, ranging tools, etc., making it easier for developers to realize product ideas. The parallel computing capability of GPU is fully utilized, and the rendering performance of large amount of data is greatly improved by combining with WebWorker multithreading technology. It supports drawing of points, lines and faces up to a million levels, and can maintain high frame rate operation.

Synchronous launch of Javascript API GL based Location data visualization API library , welcome to experience.

problem

The biggest problem of WebGL floating-point number precision is that js is 64 bit precision. js can only wear 32-bit floating-point number into the shader, the effective number is 8 bits, and the precision loss is serious.

analysis

In the basic base map, all the elements get the relative coordinates in the tile, and the coordinate range is 0-256. In each rendering, a real-time calculation of the offset of the tiles relative to the center point will be performed to calculate the tile's own matrix. In this case, the accuracy loss is relatively small, and each zoom level will load new tiles, and there will be no problem of excessive accuracy loss.

But for some coverings, such as marker, polyline and label, longitude and latitude are used. There are many decimal places after longitude and latitude. The number from js is transferred into GL gl.FLOAT It is a 32-bit floating-point number, and the decimal point can only be guaranteed to the last 4 or 5 digits. There will be serious jitter problems at level 18.

Several solutions are mentioned in the article. For example, mapbox uses the second solution, which cuts the covers such as marker, polyline and polygon according to tiles, and converts the longitude and latitude into 0-256 numbers in tile grid. In this method, each zoom transformation is re segmented according to the new mesh. Especially after level 18, such as level 22 of interior drawing, the grid is very small, resulting in a long segmentation time. Continue to try to find similar problems in mapbox: https://github.com/mapbox/mapbox-gl-js/issues/7268

mapbox also uses conversion to view space. But it's not for us.

Continue to think, the reason for this problem is that the 32-bit floating-point number effective bit is not enough. We need to find a relative coordinate as the benchmark, and other covering coordinates are based on this point. The coordinates of this relative origin retain most of the numbers, and the remaining relative coordinate numbers are as small as possible, so the effective bit is left as many decimal places as possible. And then divide the relative coordinates into two parts Math.fround (lat)，lat - Math.fround(lat); then the two parts are added in the colorator.

6.17 was executed according to this logic for the first time, and it was found that the problem of floating-point precision could not be solved when it was more than 4:00 in the morning. No. 18 discussed with angor. First, the high and low bits cannot be added directly in the shader for calculation. Although it is not possible to set the float of high type, it may be because some large number multiplication calculation is done later and the precision is eliminated. After that, the high-order low order is calculated separately and finally added, but the result is not good. It is speculated that the tile coordinate conversion is done inside, and part of the calculation is 256 x 2^n, resulting in the loss of accuracy. It is also possible that on some models, even if the floating-point number used by high P is set, it is 32-bit, according to this article https://blog.csdn.net/abcdu1/article/details/75095781 To see, this is to get 32-bit floating-point numbers https://developer.mozilla.org/en-US/docs/Web/API/WebGL_API/WebGL_best_practices

map.renderEngin.gl.getShaderPrecisionFormat( map.renderEngin.gl.VERTEX_SHADER, map.renderEngin.gl.HIGH_FLOAT )

solve

Finally from deck.gl A solution was found in to split the incoming data into a high and a low bit.

project_uCoordinateOrigin uses the latitude and longitude coordinates of the map center point

Part of the key in shaders is project_uCommonUnitsPerWorldUnit and project_uCommonUnitsPerWorldUnit2 are two uniform quantities. After tracking the code, it is found that there are calculations here:

getDistanceScales() {
        // {latitude, longitude, zoom, scale, highPrecision = false}

        let center = this.center;
        let latitude = center.lat;
        let longitude = center.lng;
        let scale = this.zoomScale(this.zoom);
        let highPrecision = true;
        // Calculate scale from zoom if not provided
        scale = scale !== undefined ? scale : this.zoomToScale(zoom);

        // assert(Number.isFinite(latitude) && Number.isFinite(longitude) && Number.isFinite(scale));
      
        const result = {};
        const worldSize = TILE_SIZE * scale;
        const latCosine = Math.cos(latitude * DEGREES_TO_RADIANS);
      
        /**
         * Number of pixels occupied by one degree longitude around current lat/lon:
           pixelsPerDegreeX = d(lngLatToWorld([lng, lat])[0])/d(lng)
             = scale * TILE_SIZE * DEGREES_TO_RADIANS / (2 * PI)
           pixelsPerDegreeY = d(lngLatToWorld([lng, lat])[1])/d(lat)
             = -scale * TILE_SIZE * DEGREES_TO_RADIANS / cos(lat * DEGREES_TO_RADIANS)  / (2 * PI)
         */
        const pixelsPerDegreeX = worldSize / 360;
        const pixelsPerDegreeY = pixelsPerDegreeX / latCosine;
      
        /**
         * Number of pixels occupied by one meter around current lat/lon:
         */
        const altPixelsPerMeter = worldSize / EARTH_CIRCUMFERENCE / latCosine;
      
        /**
         * LngLat: longitude -> east and latitude -> north (bottom left)
         * UTM meter offset: x -> east and y -> north (bottom left)
         * World space: x -> east and y -> south (top left)
         *
         * Y needs to be flipped when converting delta degree/meter to delta pixels
         */
        result.pixelsPerMeter = [altPixelsPerMeter, altPixelsPerMeter, altPixelsPerMeter];
        result.metersPerPixel = [1 / altPixelsPerMeter, 1 / altPixelsPerMeter, 1 / altPixelsPerMeter];
      
        result.pixelsPerDegree = [pixelsPerDegreeX, pixelsPerDegreeY, altPixelsPerMeter];
        result.degreesPerPixel = [1 / pixelsPerDegreeX, 1 / pixelsPerDegreeY, 1 / altPixelsPerMeter];
      
        /**
         * Taylor series 2nd order for 1/latCosine
           f'(a) * (x - a)
             = d(1/cos(lat * DEGREES_TO_RADIANS))/d(lat) * dLat
             = DEGREES_TO_RADIANS * tan(lat * DEGREES_TO_RADIANS) / cos(lat * DEGREES_TO_RADIANS) * dLat
         */
        if (highPrecision) {
          const latCosine2 = DEGREES_TO_RADIANS * Math.tan(latitude * DEGREES_TO_RADIANS) / latCosine;
          const pixelsPerDegreeY2 = pixelsPerDegreeX * latCosine2 / 2;
          const altPixelsPerDegree2 = worldSize / EARTH_CIRCUMFERENCE * latCosine2;
          const altPixelsPerMeter2 = altPixelsPerDegree2 / pixelsPerDegreeY * altPixelsPerMeter;
      
          result.pixelsPerDegree2 = [0, pixelsPerDegreeY2, altPixelsPerDegree2];
          result.pixelsPerMeter2 = [altPixelsPerMeter2, 0, altPixelsPerMeter2];
        }
      
        // Main results, used for converting meters to latlng deltas and scaling offsets
        return result;
    }

For project_uCommonUnitsPerWorldUnit is to calculate the number of tile pixels once represented in precision and latitude. For project_ For ucommonunitsperworldunit2, there is a second-order expansion of Taylor series (the more the expansion terms of Taylor series refer to the smaller the simulation value error, and the second level is used here) mainly in the shader_ uCommonUnitsPerWorldUnit + project_uCommonUnitsPerWorldUnit2 * dy do precision compensation here

There are also some doubts here. The number here is not small, and the number of significant bits is not much. Is it that there are more significant bits that can be retained by uniform? (it may also be that it does not need that high precision to convert the tile pixel coordinates. Only the tile bits of integer are needed, personal guess may not be correct)

gl.uniform3f(this.project_uCommonUnitsPerWorldUnit,distanceScles.pixelsPerDegree[0],distanceScles.pixelsPerDegree[1],distanceScles.pixelsPerDegree[2]);

As a whole, this scheme is used to solve the jitter problem caused by the loss of accuracy, which makes the accuracy basis for subsequent points, lines, surfaces and seiya.

vec2 project_offset(vec2 offset) {
      float dy = offset.y;
      // if (project_uCoordinateSystem == COORDINATE_SYSTEM_LNGLAT_AUTO_OFFSET) {
        dy = clamp(dy, -1., 1.);
      // }
      vec3 commonUnitsPerWorldUnit = project_uCommonUnitsPerWorldUnit + project_uCommonUnitsPerWorldUnit2 * dy;
      // return vec4(offset.xyz * commonUnitsPerWorldUnit, offset.w);
      return vec2(offset.xy * commonUnitsPerWorldUnit.xy);
    }

    // Return the coordinates in 3d coordinate system of v3 api, using high precision mode
    vec2 project_view_local_position3(vec2 latlngHigh, vec2 latlngLow) {
      vec2 centerCoordHigh = project_position(center.xy + center.zw, zoom);

      // Subtract high part of 64 bit value. Convert remainder to float32, preserving precision.
      float X = latlngHigh.x - center.x;
      float Y = latlngHigh.y - center.y;

      return project_offset(vec2(X + latlngLow.x, Y + latlngLow.y));

    }

Final effect: