GPU中实现反距离加权插值(IDW)_cesium如何实现反距离加权平均插值法(idw)

利用GPU实现IDW(反距离加权插值)

IDW的实现比较简单，已知插值点位比较少的情况下，可以直接遍历所有插值点，来获取临近的几个点，进行插值运算。插值点较多时，需要可以使用kd-tree来加速临近点的查找。本次仅对小数据量情况进行讨论，将IDW算法在shader中实现，可以方便的使用webgl的裁剪面和强大的并行计算。具体计算流程如下：

插值点位提交到GPU

数据提交，我使用了uniform数组变量来存储，在渲染过程中，直接提交插值点位的数据，数据类型为vec3，webgl我使用了twgl第三方库（为了偷懒）, twgl可以减少编写原生webgl代码量。测试用的数据直接使用Math.random()生成。

//points为math.random生成的测试数据
const uniforms = {
    points: points
};
twgl.setUniforms(programInfo, uniforms);
1
2
3
4
5

需要注意，shader中声明的数据类型为vec3，使用x,y,z依次来存储位置，观测值。生成测试数据时，如果需要模拟300个点位，则需要生成的数组长度为300*3。即提交到uniform中的变量为：定点数的三倍。

Shader中实现IDW算法

IDW算法的生成需要用到，屏幕的像素坐标，所以应选用Fragment Shader，实现思路比较简单，插值点位的数据传入后，通过计算最近的两个点位的距离，来计算，该像素位置的观测值。webgl中内置的gl_FragCoord，可以方便的获取屏幕像素坐标。shader中的代码如下：

#version 300 es
precision mediump float;
uniform vec3 points[300];

//存储最近的点的索引，和距离
struct interpoPoint {
  int idx;
  float dist;
};

interpoPoint queues[2];

out vec4 outcolor;

//更新队列，
void updatequeue(int idx,float dist){
    if(dist<queues[0].dist){
        queues[1]=queues[0];
        queues[0].idx=idx;
        queues[0].dist=dist;
    }else if(dist<queues[1].dist){
        queues[1].idx=idx;
        queues[1].dist=dist;
    }
}

float getCoorddist(vec2 coord,vec2 point){
    float dertx=coord.x-point.x;
    float derty=coord.y-point.y;
    return dertx*dertx+derty*derty;
}

float getCoordColor(vec2 coord){
    //前两个点位的距离更新到queues
    float dist0=getCoorddist(coord,points[0].xy);
    float dist1=getCoorddist(coord,points[1].xy);
    if(dist0<dist1){
        queues[0].idx=0;
        queues[0].dist=dist0;
        queues[1].idx=1;
        queues[1].dist=dist1;
    }else{
        queues[0].idx=1;
        queues[0].dist=dist1;
        queues[1].idx=0;
        queues[1].dist=dist0;
    }

    for(int i=2;i<300;i++){
        float dist=getCoorddist(coord,points[i].xy);
        updatequeue(i,dist);
    }

    //idw
    float denominator=(1.0/queues[0].dist)+(1.0/queues[1].dist);
    float numerator1=points[queues[0].idx].z*(1.0/queues[0].dist)/denominator;
    float numerator2=points[queues[1].idx].z*(1.0/queues[1].dist)/denominator;

    return numerator1+numerator2;
}

void main(){
    float color=getCoordColor(gl_FragCoord.xy);

    outcolor=vec4(vec3(color),1.0);
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67

效果图

效果图有些奇怪，边界特别明显，我用的指数为2，最近的点位个数为2，印象里应该是比较平滑才对，这里理论知识需要在复习一下。

现有问题

• 最大的问题是，插值点位过多时，受限于uniform数据的个数限制，大量数据无法存放到shader的数组中

  • 可以使用生成texture的方法，来存放插值点的数据，同样的存在大量点位中查找最近的k个临近点的效率问题，暂时考虑使用kdtree来实现快速查找，但是需要在shader中实现对kdtree树的查找，这个断断续续看了两个月了，还没有实现，有大神一起做的话，欢迎讨论

• idw的指数，固定为2，想要直接使用的小伙伴，可以自行修改

• idw中，临近点的个数，shader中固定为2，这里需要实现一个有序的队列，来实现大于2的点位查找。

PS: 代码已提交，github地址，有在做相关内容的小伙伴，欢迎一起讨论