WebGL着色器32位浮點數精度損失問題

問題

WebGL浮點數精度最大的問題是就是因為js是64位精度的，js往着色器里面穿的時候只能是32位浮點數，有效數是8位，精度丟失比較嚴重。

這篇文章里講了一些處理方式，但是視坐標這種方式放在我們的場景里不適用

http://www.yanglajiao.com/article/baidumap2018/80270150

 

 

 分析

在基礎底圖中，所有的要素拿到的都是瓦片里面的相對坐標，坐標范圍在0-256之間。在每次渲染時都會重新實時計算瓦片相對中心點的一個偏移來計算瓦片自己的矩陣，這種情況下精度損失比較小，而且每個zoom級別都會加載新的瓦片，不會出現精度損失過大問題。

 

但是對於一些覆蓋物，比如marker、polyline、label使用的都是經緯度，經緯度小數點后位數比較多，從js的數字傳入到gl中使用的gl.FLOAT是32位浮點數，小數點只能保證到后4位或者5位。在18級會出現嚴重的抖動問題。

 

文章中提到了幾種解決方案，像mapbox使用的是第二種方案，將覆蓋物比如marker、polyline、polygon都按照瓦片切分，經緯都轉換成瓦片網格里面的0-256數字。這種方法每次zoom變換都要按照新的網格來重新切分。尤其到了18級往后，比如室內圖22級，網格非常小，導致切分時間特別長。

繼續嘗試發現mapbox中也有類似問題：https://github.com/mapbox/mapbox-gl-js/issues/7268

mapbox這里也是使用了轉換到視空間。但這種方式並不適合我們。

 

繼續思考，實際這個問題原因是32位浮點數有效位不夠，我們要找一個相對坐標為基准，其他的覆蓋物坐標都是以這個點為基准，這個相對原點的坐標保留大部分數字，剩下的相對坐標數字盡量小，這樣有效位盡量留給更多的小數位。然后把這個相對坐標分為兩部分Math.fround（lat），lat - Math.fround(lat)；然后兩部分分別在着色器重進行計算結果在相加。

 

6.17號第一次按照這個邏輯執行了，搞到凌晨四點多，發現並不能解決浮點數精度問題。18號跟安哥討論了下，首先這個高位和低位不能直接在着色器里相加后進行計算。盡管設置了highp類型的float還是不行，這里面可能是因為后面有做了一些大數的乘法計算導致精度被消磨掉了。而后有做了高位的低位分別計算最后在相加，結果也不行，猜測是因為里面做了瓦片坐標轉換，有一部分256 x 2^n這種計算，導致精度損失。也有可能是在某些機型上即使設置了highp實際使用的浮點數也是32位的，按照這篇文章說法（ https://blog.csdn.net/abcdu1/article/details/75095781）來看，下面這個確實是得到32位浮點數 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 )

解決

最終從deck.gl中找到了一種解決方案，也是將傳入的數據拆分成一個高位和低位。

 

project_uCoordinateOrigin使用的是地圖中心點的經緯度坐標

其中着色器中的一部分關鍵是project_uCommonUnitsPerWorldUnit和project_uCommonUnitsPerWorldUnit2這兩個uniform量。跟蹤代碼后發現在這里有計算：

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;
    }

對於project_uCommonUnitsPerWorldUnit來說就是計算在精度和緯度上，一度代表的瓦片像素數目。對於project_uCommonUnitsPerWorldUnit2來說這里面用了一個泰勒級數的二階展開（咨詢了下管戈，泰勒級數展開項越多代表模擬值誤差越小，這里用到了第二級）主要是在着色器中在`project_uCommonUnitsPerWorldUnit + project_uCommonUnitsPerWorldUnit2 * dy`這里做精度補償

 

這里也有一些疑點，這里數字也不小，有效位的保留也不多，難道是uniform這種能夠保留的有效位多一些？（也可能是轉化成了瓦片像素坐標不需要那么高的精度吧。只需要整數的瓦片位，個人猜測可能不對）

 

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

整體來說使用這種方案解決精度損失引起的抖動問題，為后續的點、線、面、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);
    }

    // 返回在v3 api中的3d坐標系下的坐標, 采用高精度模式
    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));

    }

最終效果：