风场可视化：绘制粒子

引子

了解风场数据之后，接着去看如何绘制粒子。

源库：webgl-wind
Origin
My GitHub

绘制地图粒子

查看源库，发现单独有一个 Canvas 绘制地图，获取的世界地图海岸线坐标，主要格式如下：

{"type": "FeatureCollection","features": [{"type": "Feature","properties": {"scalerank": 1,"featureclass": "Coastline"},"geometry": {"type": "LineString","coordinates": [[-163.7128956777287,-78.59566741324154],// 数据省略]}},// 数据省略]
}

这些坐标对应的点连起来就可以形成整体的轮廓，主要逻辑如下：

  // 省略for (let i = 0; i < len; i++) {const coordinates = data[i].geometry.coordinates || [];const coordinatesNum = coordinates.length;for (let j = 0; j < coordinatesNum; j++) {context[j ? "lineTo" : "moveTo"](((coordinates[j][0] + 180) * node.width) / 360,((-coordinates[j][1] + 90) * node.height) / 180);}// 省略

按照 Canvas 实际的宽高度，与生成的风场图片宽高按比例映射。

绘制地图的单独逻辑示例见这里。

绘制风粒子

查看源库，单独有一个 Canvas 绘制风粒子。看源码的时候，发现其中的逻辑涉及较多状态，计划先单独弄明白绘制静态粒子的逻辑。

静态风粒子效果见示例。

先理一下实现的主要思路：

风速映射到像素颜色编码的 R 和 G 分量，由此生成了图片 W 。
创建显示用的颜色数据，并存放到纹理 T1 中。
根据粒子数，创建存储粒子索引的数据并缓冲。还创建每个粒子相关信息的数据，并存放到纹理 T2 中。
加载图片 W 并将图片数据存放到纹理 T3 中。
顶点着色器处理的时候，会根据粒子索引从纹理 T2 中获取对应数据，进行转换会生成一个位置 P 传递给片元着色器。
片元着色器根据位置 P 从图片纹理 T3 中得到数据并进行线性混合得到一个值 N ，根据 N 在颜色纹理 T1 中得到对应的颜色。

下面就看看具体的实现。

颜色数据

生成颜色数据主要逻辑：

function getColorRamp(colors) {const canvas = document.createElement("canvas");const ctx = canvas.getContext("2d");canvas.width = 256;canvas.height = 1;// createLinearGradient 用法： https://developer.mozilla.org/en-US/docs/Web/API/CanvasRenderingContext2D/createLinearGradientconst gradient = ctx.createLinearGradient(0, 0, 256, 0);for (const stop in colors) {gradient.addColorStop(+stop, colors[stop]);}ctx.fillStyle = gradient;ctx.fillRect(0, 0, 256, 1);return new Uint8Array(ctx.getImageData(0, 0, 256, 1).data);
}

这里通过创建一个渐变的 Canvas 得到数据，由于跟颜色要对应，一个颜色分量存储为 8 位二进制，总共 256 种。

Canvas 里面的数据放到纹理中，需要足够的大小：16 * 16 = 256 。这里的宽高在后面的片元着色器会用到，需要这两个地方保持一致才能达到预期结果。

this.colorRampTexture = util.createTexture(this.gl,this.gl.LINEAR,getColorRamp(colors),16,16
);

顶点数据和状态数据

主要逻辑：

set numParticles(numParticles) {const gl = this.gl;const particleRes = (this.particleStateResolution = Math.ceil(Math.sqrt(numParticles)));// 总粒子数this._numParticles = particleRes * particleRes;// 所有粒子的颜色信息const particleState = new Uint8Array(this._numParticles * 4);for (let i = 0; i < particleState.length; i++) {// 生成随机颜色，颜色会对应到图片中的位置particleState[i] = Math.floor(Math.random() * 256);}// 创建存储所有粒子颜色信息的纹理this.particleStateTexture = util.createTexture(gl,gl.NEAREST,particleState,particleRes,particleRes);// 粒子索引const particleIndices = new Float32Array(this._numParticles);for (let i = 0; i < this._numParticles; i++) particleIndices[i] = i;this.particleIndexBuffer = util.createBuffer(gl, particleIndices);
}

粒子的颜色信息会存在纹理中，这里创建了宽高相等的纹理，每个粒子颜色 RGBA 4 个分量，每个分量 8 位。注意这里生成随机颜色分量的大小范围是 [0, 256) 。

从后面逻辑可知，这里顶点数据 particleIndexBuffer 是用来辅助计算最终位置，而实际位置跟纹理有关。更加详细见下面顶点着色器的具体实现。

顶点着色器

顶点着色器和对应绑定的变量：

const drawVert = `precision mediump float;attribute float a_index;uniform sampler2D u_particles;uniform float u_particles_res;varying vec2 v_particle_pos;void main(){vec4 color=texture2D(u_particles,vec2(fract(a_index/u_particles_res),floor(a_index/u_particles_res)/u_particles_res));// 从像素的 RGBA 值解码当前粒子位置v_particle_pos=vec2(color.r / 255.0 + color.b,color.g / 255.0 + color.a);gl_PointSize = 1.0;gl_Position = vec4(2.0 * v_particle_pos.x - 1.0, 1.0 - 2.0 * v_particle_pos.y, 0, 1);}
`;// 代码省略
util.bindAttribute(gl, this.particleIndexBuffer, program.a_index, 1);
// 代码省略
util.bindTexture(gl, this.particleStateTexture, 1);
// 代码省略
gl.uniform1i(program.u_particles, 1);
// 代码省略
gl.uniform1f(program.u_particles_res, this.particleStateResolution);

从这些分散的逻辑中，找到着色器中变量对应的实际值：

a_index ： particleIndices 里面的粒子索引数据。
u_particles ：所有粒子颜色信息的纹理 particleStateTexture 。
u_particles_res ： particleStateResolution 的值，与纹理 particleStateTexture 的宽高一致，也是总粒子数的平方根，也是粒子索引数据长度的平方根。

根据这些对应值，再来看主要的处理逻辑：

vec4 color=texture2D(u_particles,vec2(fract(a_index/u_particles_res),floor(a_index/u_particles_res)/u_particles_res));

先介绍两个函数信息：

floor(x) ：返回小于等于 x 的最大整数值。
fract(x) ：返回 x - floor(x) ，即返回 x 的小数部分。

假设总粒子数是 4 ，那么 particleIndices = [0,1,2,3]、u_particles_res = 2 ，那么二维坐标依次是 vec2(0,0)、 vec2(0.5,0)、 vec2(0,0.5)、 vec2(0.5,0.5) 。这里的计算方式确保了得到的坐标都在 0 到 1 之间，这样才能在纹理 particleStateTexture 中采集到颜色信息。

这里需要注意的是 texture2D 采集返回的值范围是 [0, 1] ，具体原理见这里。

v_particle_pos=vec2(color.r / 255.0 + color.b,color.g / 255.0 + color.a);

源码注释说“从像素的 RGBA 值解码当前粒子位置”，结合前面数据来看，这样的计算方式得到分量理论范围是 [0, 256/255] ，。变量 v_particle_pos 会在片元着色器中用到。

gl_Position = vec4(2.0 * v_particle_pos.x - 1.0, 1.0 - 2.0 * v_particle_pos.y, 0, 1);

gl_Position 变量是顶点转换到裁剪空间中的坐标值，裁减空间范围 [-1.0, +1.0] ，想要显示就必须要在这个范围内，这里的计算方式达到了这个目的。

片元着色器

片元着色器和对应绑定的变量：

const drawFrag = `precision mediump float;uniform sampler2D u_wind;uniform vec2 u_wind_min;uniform vec2 u_wind_max;uniform sampler2D u_color_ramp;varying vec2 v_particle_pos;void main() {vec2 velocity = mix(u_wind_min, u_wind_max, texture2D(u_wind, v_particle_pos).rg);float speed_t = length(velocity) / length(u_wind_max);vec2 ramp_pos = vec2(fract(16.0 * speed_t),floor(16.0 * speed_t) / 16.0);gl_FragColor = texture2D(u_color_ramp, ramp_pos);}
`;// 代码省略
util.bindTexture(gl, this.windTexture, 0);
// 代码省略
gl.uniform1i(program.u_wind, 0); // 风纹理数据
// 代码省略
util.bindTexture(gl, this.colorRampTexture, 2);
// 代码省略
gl.uniform1i(program.u_color_ramp, 2); // 颜色数据
// 代码省略
gl.uniform2f(program.u_wind_min, this.windData.uMin, this.windData.vMin);
gl.uniform2f(program.u_wind_max, this.windData.uMax, this.windData.vMax);

从这些分散的逻辑中，找到着色器中变量对应的实际值：

u_wind ：风场图片生成的纹理 windTexture 。
u_wind_min ：风场数据分量最小值。
u_wind_max ：风场数据分量最大值。
u_color_ramp ：创建的颜色纹理 colorRampTexture 。
v_particle_pos ：在顶点着色器里面生成的位置。

vec2 velocity = mix(u_wind_min, u_wind_max, texture2D(u_wind, v_particle_pos).rg);
float speed_t = length(velocity) / length(u_wind_max);

先介绍内置函数：

mix(x, y, a) ：会返回 x 和 y 的线性混合，计算方式等同于 x*(1-a) + y*a 。

velocity 的值确保在 u_wind_min 和 u_wind_max 之间，那么 speed_t 的结果一定是小于或等于 1 。根据 speed_t 按照一定规则得到位置 ramp_pos ，在颜色纹理 colorRampTexture 中得到输出到屏幕的颜色。

绘制

在以上逻辑准备好后，绘制按照正常的顺序执行即可。

虽然是绘制静态的粒子，但在单独抽离的过程中发现，不同数量的粒子，如果只执行一次绘制 wind.draw() ，可能无法完成绘制。

静态风粒子效果见示例。

小结

经过了上面代码逻辑分析后，再回头看看一开始的主要思路，换个方式表述一下：

根据需要显示的粒子数，随机初始化每一个粒子的颜色编码信息并存放到纹理 T2 中；创建最终显示粒子的颜色纹理 T1 ；加载风速生成的图片 W 并存放到纹理 T3 中。
最终的目的是从颜色纹理 T1 中获取到颜色并显示，这个过程的方式就是根据纹理 T2 从纹理 T3 中找到一个对应的风速映射点，然后根据这个点从 T1 找到对应的显示颜色。

感觉比一开始的主要思路好懂了一些，但还是有一些疑问。

为什么不直接将纹理 T3 与颜色纹理 T1 关联映射？

目前这里只是整个风场可视化逻辑的一部分重现，回头看看完整的实现效果：是动态的。那么为了跟踪每一个粒子的移动，增加一个相关记录变量的实现方式，个人感觉在逻辑上会更加清晰一些，纹理 T2 主要是用来记录粒子数及状态，后续会继续深入相关逻辑。

顶点着色器中用于纹理采样的二维向量计算依据是什么？

对应的就是为什么用下面这个逻辑：

vec2(fract(a_index/u_particles_res),floor(a_index/u_particles_res)/u_particles_res
)

在前面的具体解释中有说，这样的计算方式确保了得到的坐标都在 0 到 1 之间，但能生成这个范围内的方式应该不止这一种，为什么偏偏选这种，个人也不太清楚。后面片元着色器中计算最终位置 ramp_pos 时也用了这样类似的方式。

片元着色器本来就已经得到一个位置了，为什么还要计算 `velocity` 重新得到一个位置？

也就是为什么要有下面这段逻辑：

vec2 velocity = mix(u_wind_min, u_wind_max, texture2D(u_wind, v_particle_pos).rg);
float speed_t = length(velocity) / length(u_wind_max);

从顶点着色器中得到位置 v_particle_pos 是基于随机生成的颜色纹理 T2 得到的，前面有说分量值计算理论范围是 [0, 256/255] ，无法保证一定可以在风场图片中找到对应的点，那么通过 mix 函数就可以生成一种关联。

片元着色器中计算 `ramp_pos` 相乘的系数为什么是 16.0 ？

就是下面这段逻辑：

vec2 ramp_pos = vec2(fract(16.0 * speed_t),floor(16.0 * speed_t) / 16.0);

通过尝试发现这里的 16.0 是跟前面生成最终显示用的颜色纹理 T1 的宽高需要一致，猜测这样一致才能达到均匀的效果。

参考资料

How I built a wind map with WebGL