fix: bilinear resample dont handle nodata

hongfaqiu · Sep 22, 2024 · 6b1d709 · 6b1d709
1 parent fb5ed2c
commit 6b1d709
Show file tree

Hide file tree

Showing 3 changed files with 40 additions and 22 deletions.
diff --git a/packages/TIFFImageryProvider/src/TIFFImageryProvider.ts b/packages/TIFFImageryProvider/src/TIFFImageryProvider.ts
@@ -607,6 +607,7 @@ export class TIFFImageryProvider {
         window: [x0, y0, x1, y1],
         method: this.renderOptions.resampleMethod,
         buffer,
+        nodata: this.noData,
       })));
 
       return {

diff --git a/packages/TIFFImageryProvider/src/helpers/resample.ts b/packages/TIFFImageryProvider/src/helpers/resample.ts
@@ -9,6 +9,8 @@ export type ResampleDataOptions = {
   window: [number, number, number, number];
   method?: 'bilinear' | 'nearest';
   buffer?: number;
+  /** No-data value */
+  nodata?: number; // 新增 nodata 属性
 }
 
 export function resampleNearest(data: TypedArray, options: ResampleDataOptions) {
@@ -33,7 +35,7 @@ export function resampleNearest(data: TypedArray, options: ResampleDataOptions)
 
 export function resampleData(data: TypedArray, options: ResampleDataOptions) {
   const { method = 'nearest' } = options;
-  
+
   switch (method) {
     case "nearest":
       return resampleNearest(data, options);
@@ -68,30 +70,43 @@ export function resampleBilinear(data: TypedArray, options: ResampleDataOptions)
   const effectiveSourceWidth = sourceWidth - 2 * buffer;
   const effectiveSourceHeight = sourceHeight - 2 * buffer;
 
+  const invTargetWidth = 1 / targetWidth;
+  const invTargetHeight = 1 / targetHeight;
+
   for (let y = 0; y < targetHeight; y++) {
-    const rawY = effectiveSourceHeight * (y0 + y / targetHeight * windowHeight) + buffer;
+    const yRatio = y * invTargetHeight;
+    const yMapped = y0 + yRatio * windowHeight;
+    const rawY = effectiveSourceHeight * yMapped + buffer;
     const yl = Math.floor(rawY);
     const yh = Math.min(Math.ceil(rawY), sourceHeight - buffer - 1);
+    const ty = rawY - yl; // Equivalent to rawY % 1 but more efficient
 
     for (let x = 0; x < targetWidth; x++) {
-      const rawX = effectiveSourceWidth * (x0 + x / targetWidth * windowWidth) + buffer;
-      const tx = rawX % 1;
-
+      const xRatio = x * invTargetWidth;
+      const xMapped = x0 + xRatio * windowWidth;
+      const rawX = effectiveSourceWidth * xMapped + buffer;
       const xl = Math.floor(rawX);
       const xh = Math.min(Math.ceil(rawX), sourceWidth - buffer - 1);
-
-      const ll = data[(yl * sourceWidth) + xl];
-      const hl = data[(yl * sourceWidth) + xh];
-      const lh = data[(yh * sourceWidth) + xl];
-      const hh = data[(yh * sourceWidth) + xh];
-
-      const value = lerp(
-        lerp(ll, hl, tx),
-        lerp(lh, hh, tx),
-        rawY % 1,
-      );
-      newArray[(y * targetWidth) + x] = value;
+      const tx = rawX - xl; // Equivalent to rawX % 1 but more efficient
+
+      const ll = data[yl * sourceWidth + xl];
+      const hl = data[yl * sourceWidth + xh];
+      const lh = data[yh * sourceWidth + xl];
+      const hh = data[yh * sourceWidth + xh];
+
+      // Check if any of the four neighboring pixels is nodata
+      if (ll === options.nodata || hl === options.nodata || lh === options.nodata || hh === options.nodata) {
+        newArray[y * targetWidth + x] = options.nodata;
+        continue;
+      }
+
+      // Perform bilinear interpolation
+      const v0 = ll * (1 - tx) + hl * tx;
+      const v1 = lh * (1 - tx) + hh * tx;
+      const value = v0 * (1 - ty) + v1 * ty;
+      newArray[y * targetWidth + x] = value;
     }
   }
+
   return newArray;
 }
diff --git a/vite-example/src/index.ts b/vite-example/src/index.ts
@@ -1,7 +1,7 @@
 import { ArcGisMapServerImageryProvider, ImageryLayer, Viewer } from 'cesium';
 import './index.css';
 import TIFFImageryProvider from 'tiff-imagery-provider';
-import proj4 from 'proj4-fully-loaded'; 
+import proj4 from 'proj4-fully-loaded';
 
 const viewer = new Viewer('cesiumContainer', {
   baseLayer: ImageryLayer.fromProviderAsync(ArcGisMapServerImageryProvider.fromUrl('https://services.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer', {
@@ -20,13 +20,13 @@ const viewer = new Viewer('cesiumContainer', {
   orderIndependentTranslucency: false,
 });
 
-TIFFImageryProvider.fromUrl('/cogtif.tif', {
+TIFFImageryProvider.fromUrl('/zeta.tif', {
   enablePickFeatures: true,
   projFunc: (code) => {
     if (![4326, 3857, 900913].includes(code)) {
       {
         try {
-          let prj = proj4("EPSG:4326", `EPSG:${code}`, )
+          let prj = proj4("EPSG:4326", `EPSG:${code}`,)
           if (prj) return {
             project: prj.forward,
             unproject: prj.inverse
@@ -40,16 +40,18 @@ TIFFImageryProvider.fromUrl('/cogtif.tif', {
   },
   renderOptions: {
     single: {
+      band: 1,
       colorScale: 'rainbow'
-    }
+    },
+    resampleMethod: 'bilinear'
   },
 }).then((provider) => {
   console.log(provider);
   const imageryLayer = viewer.imageryLayers.addImageryProvider(provider as any);
   const legend = document.getElementById("legend") as HTMLImageElement;
   const img = provider.plot?.colorScaleCanvas.toDataURL();
   legend.src = img;
-  
+
   viewer.flyTo(imageryLayer, {
     duration: 1,
   });