流动的点描：基于 WCSPH 流体动力学的图像风格化实践

最近我在做一个和点描（Stippling）相关的小项目：根据图片的灰度信息生成大量点的位置。实现这一目标的方法其实很多，比如误差扩散、Voronoi 质心剖分等等。每种方法都有它独特的美感和挑战。

在尝试这些算法的过程中，我突然想到：这些生成的点不就是一个个粒子吗？如果把它们看作粒子，那是不是可以借用平滑粒子流体动力学（SPH）的思想来探索新的可能性？或许能让点描的分布更自然、更具动态感。

虽然这只是一个实验性的想法，但正是这种跨领域的联想让项目变得更有趣。

平滑粒子自由流体动力学（SPH）

物理量与物理量的导数

在 SPH 方法中，空间中的每个粒子都“携带”特定的物理量。为了近似空间中任意位置 $\boldsymbol x$ 处的物理属性 $\boldsymbol A(\boldsymbol x)$，我们可以对邻域内的粒子进行采样，并利用核函数 $W$（Kernel）进行加权求和：

$$\boldsymbol A(\boldsymbol x) \approx \sum_i \boldsymbol A_i \frac{m_i}{\rho_i} W(\|\boldsymbol x - \boldsymbol x_i\|, h) \tag {1}$$

其中，$\boldsymbol A_i$ 是粒子 $i$ 所携带的物理属性，$m_i$ 与 $\rho_i$ 分别为其质量与密度，$h$ 为平滑长度（Smoothing Length），决定了核函数的影响范围。

在流体仿真中，我们通常更关注物理量对空间位置的导数。为了保证计算的数值稳定性，通常采用如下对称形式的经验公式来计算粒子 $i$ 处的梯度：

$$\nabla \boldsymbol A_i = \rho_i \sum_j m_j \left( \frac{\boldsymbol A_i}{\rho_i^2} + \frac{\boldsymbol A_j}{\rho_j^2} \right) \nabla_{\boldsymbol x_i} W(\|\boldsymbol x_i - \boldsymbol x_j\|, h) \tag {2}$$

需要注意的是： 这种形式并非数学意义上的严格导数，但在物理仿真中至关重要。它通过强制算子的对称性，确保了粒子间相互作用力的守恒（如动量守恒）。如果导数形式不对称，系统可能会产生非物理的数值驱动力，导致能量不守恒或模拟崩溃。

实现步骤

本项目的整体实现流程与流体动力学仿真基本一致，主要分为以下四个核心步骤：

1. 密度场计算：遍历每个粒子，根据其邻域粒子的分布计算局部密度 $\rho$。
2. 受力与加速度求解：基于密度场和状态方程（EOS）计算压力梯度，进而导出每个粒子的加速度 $\boldsymbol a$。
3. 速度更新：利用数值积分（如显式欧拉法）根据加速度迭代粒子的瞬时速度 $\boldsymbol v$。
4. 位置平移：根据更新后的速度调整粒子坐标，完成空间位置的迭代。

密度场计算

在本方案中，我采用了 弱可压缩 SPH (WCSPH) 的相关理论。首先，根据粒子间的相对位置，利用公式$(1)$近似计算连续的密度场：

$$\rho(\boldsymbol x) \approx \sum_i m_i W(\|\boldsymbol x - \boldsymbol x_i\|, h)$$

受力与加速度求解

随后，需要根据密度场推导速度对时间的导数（即加速度 $\boldsymbol a$）。在流体力学中，粒子的运动方程由下式给出：

$$\frac{\mathrm{D} \boldsymbol v}{\mathrm{D} t} = -\frac{1}{\rho} \nabla p + \boldsymbol g \tag {3}$$

其中，$p$ 为标量压力场，$\boldsymbol g$ 代表重力或其他外部驱动力，$\rho$ 即为前文求解的局部密度。为了闭合方程组，我们引入状态方程 (Equation of State, EOS) 来显式计算压力：

$$p = B \left( \left( \frac{\rho}{\rho_0} \right)^\gamma - 1 \right) \tag {4}$$

在该公式中，$B$ 为体积模量 (Bulk Modulus)，用于衡量流体抵抗压缩的能力，其作用类似于弹簧质点系统中的杨氏模量。$\rho_0$ 为参考密度 (Rest Density)。从物理直观上理解：当局部密度 $\rho$ 大于 $\rho_0$ 时，粒子感受到正向压力并趋于向外扩张；反之，当密度低于参考密度时，产生的负压（吸引力）会促使粒子向内收缩，从而维持系统的质量分布平衡。

为了使迭代后的粒子分布与目标图像的灰度分布一致，我将图像的灰度值映射为粒子的参考密度 $\rho_0$。具体实现上，通过对粒子位置处的灰度图进行双线性插值来确定该点的 $\rho_0$。此外，还可以对原始灰度值进行非线性变换，以灵活调整最终生成的粒子疏密程度。

注意到公式 $(3)$ 中涉及的是压力梯度项，因此我们利用公式 $(2)$ 给出压力梯度的离散化表达：

$$\nabla p_i = \rho_i \sum_j m_j \left( \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2} \right) \nabla_{\boldsymbol x_i} W(\|\boldsymbol x_i - \boldsymbol x_j\|, h)$$

通过该式，我们便得到了物质导数 $\frac{\mathrm{D} \boldsymbol v}{\mathrm{D} t}$ 中的压力贡献项（即内力项）。

对于第二项外力 $\boldsymbol g$，我们可以直接对目标灰度图求取梯度。引入该项的物理意义在于：通过灰度梯度引导粒子运动，从而有效减少粒子从高灰度（高目标密度）区域向低灰度区域的非预期扩散，强化粒子分布对图像特征的捕捉能力。当然这也只是一个可选项，在这个任务中 $\boldsymbol g$ 可以直接设置为0。

顺带一提，我在这里选用的核函数（Kernel Function）为经典的三次样条核（Cubic Spline Kernel），其数学表达式如下：

$$W(q, h) = \sigma_d \begin{cases} 6(q^3 - q^2) + 1 & 0 \leq q \leq \frac{1}{2} \\ 2(1 - q)^3 & \frac{1}{2} < q \leq 1 \\ 0 & \text{otherwise} \end{cases}$$

其中，$q = \frac{\|\boldsymbol r\|}{h}$ 代表粒子间距与平滑长度的比值。$\sigma_d$ 为保证积分单位化的归一化系数。由于本项目处理的是二维空间情况，对应的归一化系数取值为 $\sigma_d = \frac{40}{7\pi h^2}$。

该核函数具有良好的 紧支集（Compact Support） 特性，即当粒子间距超过 $h$ 时，相互作用力自动降为零，这极大地优化了近邻搜索的计算效率。

速度与位置更新

在时间积分方案上，本项目直接采用 显式欧拉法（Explicit Euler Method） 来更新粒子的运动状态。其离散迭代公式如下：

1. 速度更新：

   $$\boldsymbol v_{t+1} = \boldsymbol v_{t} + \Delta t \frac{\mathrm{D} \boldsymbol v}{\mathrm{D} t}$$

2. 位置更新：

   $$\boldsymbol x_{t+1} = \boldsymbol x_{t} + \Delta t \boldsymbol v_{t+1}$$

至此，系统完成了一次完整的时间步迭代。

效果展示

让我们来看看初步的实验结果：

emmm……从动图中可以看出，粒子虽然能够成功地向高灰度（高密度目标）区域聚集，但……总感觉哪里有点怪怪的？相对于CVT和误差扩散生成的点集视觉上并不均匀，点会更加倾向于进入高密度区域，并在边缘聚集。

尽管如此，这样的结果也表明了前面相关步骤是正确可行的，也算了对流体动力学有了一个初步的认知，是一次蛮有意思的尝试。

参考文献

1. Smoothed Particle Hydrodynamics Techniques for the Physics Based Simulation of Fluids and Solids
2. Games201