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

木的樹發表於2019-07-29
Web

問題

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

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

 分析

在基礎底圖中，所有的要素拿到的都是瓦片裡面的相對座標，座標範圍在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));

    }

最終效果：

 

 

相關文章