以下內(nèi)容轉(zhuǎn)載自木的樹的文章《WebGL著色器32位浮點(diǎn)數(shù)精度損失問(wèn)題》

作者：木的樹

鏈接：https://www.cnblogs.com/dojo-lzz/p/11250327.html

來(lái)源：博客園

著作權(quán)歸作者所有白群。商業(yè)轉(zhuǎn)載請(qǐng)聯(lián)系作者獲得授權(quán)贾陷，非商業(yè)轉(zhuǎn)載請(qǐng)注明出處彼宠。

前言

Javascript API GL是基于WebGL技術(shù)打造的3D版地圖API堂竟，3D化的視野更為自由，交互更加流暢缴挖。
提供豐富的功能接口亲桦，包括點(diǎn)殉摔、線炮车、面繪制舵变，自定義圖層、個(gè)性化樣式及繪圖瘦穆、測(cè)距工具等棋傍，使開發(fā)者更加容易的實(shí)現(xiàn)產(chǎn)品構(gòu)思。
充分發(fā)揮GPU的并行計(jì)算能力难审，同時(shí)結(jié)合WebWorker多線程技術(shù)，大幅度提升了大數(shù)據(jù)量的渲染性能亿絮。最高支持百萬(wàn)級(jí)點(diǎn)告喊、線、面繪制派昧，同時(shí)可以保持高幀率運(yùn)行黔姜。

同步推出基于Javascript API GL的 位置數(shù)據(jù)可視化API庫(kù)，歡迎體驗(yàn)蒂萎。

問(wèn)題

WebGL浮點(diǎn)數(shù)精度最大的問(wèn)題是就是因?yàn)閖s是64位精度的秆吵，js往著色器里面穿的時(shí)候只能是32位浮點(diǎn)數(shù)，有效數(shù)是8位五慈，精度丟失比較嚴(yán)重纳寂。

分析

在基礎(chǔ)底圖中，所有的要素拿到的都是瓦片里面的相對(duì)坐標(biāo)泻拦，坐標(biāo)范圍在0-256之間毙芜。在每次渲染時(shí)都會(huì)重新實(shí)時(shí)計(jì)算瓦片相對(duì)中心點(diǎn)的一個(gè)偏移來(lái)計(jì)算瓦片自己的矩陣，這種情況下精度損失比較小争拐，而且每個(gè)zoom級(jí)別都會(huì)加載新的瓦片腋粥，不會(huì)出現(xiàn)精度損失過(guò)大問(wèn)題。

但是對(duì)于一些覆蓋物架曹，比如marker隘冲、polyline、label使用的都是經(jīng)緯度绑雄，經(jīng)緯度小數(shù)點(diǎn)后位數(shù)比較多展辞，從js的數(shù)字傳入到gl中使用的gl.FLOAT是32位浮點(diǎn)數(shù)，小數(shù)點(diǎn)只能保證到后4位或者5位绳慎。在18級(jí)會(huì)出現(xiàn)嚴(yán)重的抖動(dòng)問(wèn)題纵竖。

文章中提到了幾種解決方案漠烧，像mapbox使用的是第二種方案，將覆蓋物比如marker靡砌、polyline已脓、polygon都按照瓦片切分，經(jīng)緯都轉(zhuǎn)換成瓦片網(wǎng)格里面的0-256數(shù)字通殃。這種方法每次zoom變換都要按照新的網(wǎng)格來(lái)重新切分度液。尤其到了18級(jí)往后，比如室內(nèi)圖22級(jí)画舌，網(wǎng)格非常小堕担，導(dǎo)致切分時(shí)間特別長(zhǎng)。
繼續(xù)嘗試發(fā)現(xiàn)mapbox中也有類似問(wèn)題：https://github.com/mapbox/mapbox-gl-js/issues/7268

mapbox這里也是使用了轉(zhuǎn)換到視空間曲聂。但這種方式并不適合我們霹购。

繼續(xù)思考，實(shí)際這個(gè)問(wèn)題原因是32位浮點(diǎn)數(shù)有效位不夠朋腋，我們要找一個(gè)相對(duì)坐標(biāo)為基準(zhǔn)齐疙，其他的覆蓋物坐標(biāo)都是以這個(gè)點(diǎn)為基準(zhǔn)，這個(gè)相對(duì)原點(diǎn)的坐標(biāo)保留大部分?jǐn)?shù)字旭咽，剩下的相對(duì)坐標(biāo)數(shù)字盡量小贞奋，這樣有效位盡量留給更多的小數(shù)位。然后把這個(gè)相對(duì)坐標(biāo)分為兩部分Math.fround（lat）穷绵，lat - Math.fround(lat)轿塔；然后兩部分分別在著色器重進(jìn)行計(jì)算結(jié)果在相加。

6.17號(hào)第一次按照這個(gè)邏輯執(zhí)行了仲墨，搞到凌晨四點(diǎn)多勾缭，發(fā)現(xiàn)并不能解決浮點(diǎn)數(shù)精度問(wèn)題。18號(hào)跟安哥討論了下宗收，首先這個(gè)高位和低位不能直接在著色器里相加后進(jìn)行計(jì)算漫拭。盡管設(shè)置了highp類型的float還是不行，這里面可能是因?yàn)楹竺嬗凶隽艘恍┐髷?shù)的乘法計(jì)算導(dǎo)致精度被消磨掉了混稽。而后有做了高位的低位分別計(jì)算最后在相加采驻，結(jié)果也不行，猜測(cè)是因?yàn)槔锩孀隽送咂鴺?biāo)轉(zhuǎn)換匈勋，有一部分256 x 2^n這種計(jì)算礼旅，導(dǎo)致精度損失。也有可能是在某些機(jī)型上即使設(shè)置了highp實(shí)際使用的浮點(diǎn)數(shù)也是32位的洽洁，按照這篇文章說(shuō)法https://blog.csdn.net/abcdu1/article/details/75095781來(lái)看痘系，下面這個(gè)確實(shí)是得到32位浮點(diǎn)數(shù)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中找到了一種解決方案，也是將傳入的數(shù)據(jù)拆分成一個(gè)高位和低位饿自。

project_uCoordinateOrigin使用的是地圖中心點(diǎn)的經(jīng)緯度坐標(biāo)

其中著色器中的一部分關(guān)鍵是project_uCommonUnitsPerWorldUnit和project_uCommonUnitsPerWorldUnit2這兩個(gè)uniform量汰翠。跟蹤代碼后發(fā)現(xiàn)在這里有計(jì)算：

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

對(duì)于project_uCommonUnitsPerWorldUnit來(lái)說(shuō)就是計(jì)算在精度和緯度上龄坪，一度代表的瓦片像素?cái)?shù)目。對(duì)于project_uCommonUnitsPerWorldUnit2來(lái)說(shuō)這里面用了一個(gè)泰勒級(jí)數(shù)的二階展開（咨詢了下管戈复唤，泰勒級(jí)數(shù)展開項(xiàng)越多代表模擬值誤差越小健田，這里用到了第二級(jí)）主要是在著色器中在project_uCommonUnitsPerWorldUnit + project_uCommonUnitsPerWorldUnit2 * dy這里做精度補(bǔ)償

這里也有一些疑點(diǎn)，這里數(shù)字也不小佛纫，有效位的保留也不多妓局，難道是uniform這種能夠保留的有效位多一些？（也可能是轉(zhuǎn)化成了瓦片像素坐標(biāo)不需要那么高的精度吧呈宇。只需要整數(shù)的瓦片位好爬，個(gè)人猜測(cè)可能不對(duì)）

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

整體來(lái)說(shuō)使用這種方案解決精度損失引起的抖動(dòng)問(wèn)題，為后續(xù)的點(diǎn)甥啄、線存炮、面、seiya都做了精度基礎(chǔ)蜈漓。

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坐標(biāo)系下的坐標(biāo), 采用高精度模式
    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));

    }

最終效果：