桑原滤波 (Kuwahara Filter) - 边缘保持的图像平滑算法

视频教程

本文配套视频讲解：Bilibili - 桑原滤波原理与实现

桑原滤波效果

示意图：桑原滤波在保持边缘的同时对图像进行平滑处理

一、什么是桑原滤波？

桑原滤波（Kuwahara Filter）是一种边缘保持的非线性滤波算法，由日本学者桑原幸夫（Yukio Kuwahara）于 1970 年代提出。

与高斯模糊等传统线性滤波不同，桑原滤波在平滑图像的同时，能够很好地保持图像的边缘信息，不会导致边缘模糊。这使得它在风格化渲染、艺术化图像处理、卡通渲染等领域有着广泛的应用。

核心思想

桑原滤波的核心思想非常巧妙：在图像的局部区域中，选择方差最小的子区域来计算输出像素值。

为什么选择方差最小的子区域？

• 方差小表示该区域颜色/亮度变化小，说明区域内没有边缘
• 在没有边缘的区域进行平滑，自然就不会模糊边缘

---

二、数学原理

1. 基础统计量

在介绍算法之前，先回顾两个基本的统计量：

均值（Mean）：表示区域内像素的平均亮度/颜色

$$\bar{x} = \frac{1}{n} \sum_{i=1}^{n} x_i$$

方差（Variance）：表示区域内像素的变化程度

$$\sigma^2 = \frac{1}{n} \sum_{i=1}^{n} (x_i - \bar{x})^2$$

方差越小，区域内像素越相似，边缘存在的可能性越低。

2. 子区域划分

对于一个半径为 $r$ 的桑原滤波核，通常将其划分为 4 个重叠的方形子区域：

每个子区域的大小为 $(r+1) \times (r+1)$，四个子区域在中心处重叠。

3. 算法步骤

对于每个像素 $(x, y)$：

1. 划分区域：以该像素为中心，划分出 4 个重叠的子区域 $Q_1, Q_2, Q_3, Q_4$

2. 计算统计量：对每个子区域计算均值 $\mu_i$ 和方差 $\sigma_i^2$

   $$\mu_i = \frac{1}{n} \sum_{p \in Q_i} I(p)$$

   $$\sigma_i^2 = \frac{1}{n} \sum_{p \in Q_i} (I(p) - \mu_i)^2$$

3. 选择区域：找到方差最小的子区域

   $$i_{min} = \arg\min_{i=1..4} \sigma_i^2$$

4. 输出结果：将该子区域的均值作为输出像素值

   $$I_{out}(x, y) = \mu_{i_{min}}$$

---

三、GLSL 完整实现

1. 基础版本

// 桑原滤波 - 基础版本
// radius: 滤波半径，推荐值 2-4
vec3 kuwahara_filter(sampler2D tex, vec2 uv, float radius, vec2 tex_size) {
    vec2 texel = 1.0 / tex_size;
    
    float min_variance = 1e10;
    vec3 result = vec3(0.0);
    
    // 四个子区域的偏移
    vec2 offsets[4] = vec2[](
        vec2(-radius, -radius),  // 左上
        vec2(0.0,      -radius),  // 右上
        vec2(-radius,  0.0),      // 左下
        vec2(0.0,       0.0)       // 右下
    );
    
    // 遍历四个子区域
    for (int q = 0; q < 4; q++) {
        vec3 sum = vec3(0.0);
        vec3 sum_sq = vec3(0.0);
        float count = 0.0;
        
        // 遍历子区域内的像素
        for (float i = 0.0; i <= radius; i++) {
            for (float j = 0.0; j <= radius; j++) {
                vec2 offset = offsets[q] + vec2(i, j);
                vec3 color = texture(tex, uv + offset * texel).rgb;
                sum += color;
                sum_sq += color * color;
                count += 1.0;
            }
        }
        
        // 计算均值和方差
        vec3 mean = sum / count;
        vec3 variance = (sum_sq / count) - mean * mean;
        
        // 使用亮度方差来选择区域
        float luminance_variance = dot(variance, vec3(0.299, 0.587, 0.114));
        
        // 选择方差最小的区域
        if (luminance_variance < min_variance) {
            min_variance = luminance_variance;
            result = mean;
        }
    }
    
    return result;
}

2. 优化版本 - 使用分离的方差计算

基础版本需要在循环内计算方差，我们可以优化一下计算方式：

// 桑原滤波 - 优化版本
// 使用快速方差计算，减少运算量
vec3 kuwahara_filter_fast(sampler2D tex, vec2 uv, float radius, vec2 tex_size) {
    vec2 texel = 1.0 / tex_size;
    
    // 预先计算所有需要的像素
    int size = int(radius * 2.0 + 1.0);
    
    float min_var = 1e10;
    vec3 result = vec3(0.0);
    
    // 四个子区域
    for (int qy = 0; qy < 2; qy++) {
        for (int qx = 0; qx < 2; qx++) {
            vec3 sum = vec3(0.0);
            vec3 sum_sq = vec3(0.0);
            
            float start_x = float(qx) * radius - radius;
            float start_y = float(qy) * radius - radius;
            
            for (float y = 0.0; y <= radius; y++) {
                for (float x = 0.0; x <= radius; x++) {
                    vec2 offset = vec2(start_x + x, start_y + y);
                    vec3 c = texture(tex, uv + offset * texel).rgb;
                    sum += c;
                    sum_sq += c * c;
                }
            }
            
            float n = (radius + 1.0) * (radius + 1.0);
            vec3 mean = sum / n;
            vec3 var = (sum_sq / n) - mean * mean;
            
            float luma_var = dot(var, vec3(0.299, 0.587, 0.114));
            
            if (luma_var < min_var) {
                min_var = luma_var;
                result = mean;
            }
        }
    }
    
    return result;
}

3. 在 Unity Shader 中使用

// Unity URP 桑原滤波 Post Process
half4 Frag(Varyings i) : SV_Target {
    float2 uv = i.uv;
    
    // 参数控制
    float radius = _Radius;          // 滤波半径
    float strength = _Strength;      // 滤波强度
    
    // 原始颜色
    half3 original = SAMPLE_TEXTURE2D(_MainTex, sampler_MainTex, uv).rgb;
    
    // 桑原滤波结果
    half3 filtered = kuwahara_filter(_MainTex, uv, radius, _ScreenParams.xy);
    
    // 混合原始和滤波结果
    half3 result = lerp(original, filtered, strength);
    
    return half4(result, 1.0);
}

---

四、算法变种与改进

1. 广义桑原滤波（Generalized Kuwahara）

原始的桑原滤波只使用 4 个矩形区域，广义版本可以使用更多方向的区域：

// 广义桑原滤波 - 8 方向版本
vec3 generalized_kuwahara(sampler2D tex, vec2 uv, float radius, vec2 tex_size) {
    vec2 texel = 1.0 / tex_size;
    float min_var = 1e10;
    vec3 result = vec3(0.0);
    
    // 8 个方向的区域
    for (int dir = 0; dir < 8; dir++) {
        float angle = float(dir) * PI / 4.0;
        vec2 dir_vec = vec2(cos(angle), sin(angle));
        
        vec3 sum = vec3(0.0);
        vec3 sum_sq = vec3(0.0);
        float count = 0.0;
        
        // 在该方向扇形区域内采样
        for (float r = 1.0; r <= radius; r++) {
            for (float a = -PI/8.0; a <= PI/8.0; a += PI/16.0) {
                vec2 offset = r * vec2(
                    cos(angle + a), 
                    sin(angle + a)
                );
                vec3 c = texture(tex, uv + offset * texel).rgb;
                sum += c;
                sum_sq += c * c;
                count += 1.0;
            }
        }
        
        vec3 mean = sum / count;
        vec3 var = (sum_sq / count) - mean * mean;
        float luma_var = dot(var, vec3(0.299, 0.587, 0.114));
        
        if (luma_var < min_var) {
            min_var = luma_var;
            result = mean;
        }
    }
    
    return result;
}

2. 各向异性桑原滤波

根据图像的局部结构调整滤波方向：

// 各向异性桑原滤波
// 根据局部梯度方向自适应调整区域方向
vec3 anisotropic_kuwahara(sampler2D tex, vec2 uv, float radius, vec2 tex_size) {
    vec2 texel = 1.0 / tex_size;
    
    // 1. 计算局部梯度
    vec3 gx = texture(tex, uv + vec2(texel.x, 0)).rgb 
            - texture(tex, uv - vec2(texel.x, 0)).rgb;
    vec3 gy = texture(tex, uv + vec2(0, texel.y)).rgb 
            - texture(tex, uv - vec2(0, texel.y)).rgb;
    
    float dx = dot(gx, vec3(0.299, 0.587, 0.114));
    float dy = dot(gy, vec3(0.299, 0.587, 0.114));
    
    // 梯度方向
    float angle = atan(dy, dx);
    
    // 2. 沿梯度垂直方向进行滤波
    vec2 tangent = vec2(cos(angle + PI/2.0), sin(angle + PI/2.0));
    vec2 normal = vec2(cos(angle), sin(angle));
    
    // 3. 在垂直梯度的方向上选择最平滑的区域
    // ... (实现略)
    
    return vec3(0.0);
}

---

五、性能对比

不同半径的性能开销

滤波半径	像素采样数	相对性能	效果质量
r = 2	16	1.0x	⭐⭐⭐
r = 3	36	0.45x	⭐⭐⭐⭐
r = 4	64	0.25x	⭐⭐⭐⭐⭐
r = 5	100	0.16x	⭐⭐⭐⭐⭐

与其他滤波算法对比

算法	边缘保持	平滑效果	性能	适用场景
高斯模糊	❌ 差	✅ 好	⚡ 极快	通用平滑
中值滤波	⭐ 一般	⭐ 一般	🐌 慢	去除椒盐噪声
双边滤波	✅ 好	✅ 好	🐢 较慢	保边去噪
桑原滤波	✅✅ 极好	⭐ 一般	🐢 较慢	风格化、艺术渲染

性能优化建议

1. 降采样处理：先将图像降采样，滤波后再升采样，性能提升 4-16 倍
2. 分离滤波：桑原滤波不可分离，但可以先做快速的预滤波
3. 合理选择半径：r = 2 或 r = 3 通常是效果和性能的最佳平衡点

---

六、应用场景

1. 卡通渲染 / 风格化渲染

桑原滤波是实现 NPR（非真实感渲染）效果的常用算法：

// 卡通渲染管线中的桑原滤波
vec3 toon_rendering(vec3 color, vec2 uv) {
    // 1. 桑原滤波平滑颜色
    vec3 smoothed = kuwahara_filter(_MainTex, uv, 3.0, _ScreenSize);
    
    // 2. 颜色量化
    smoothed = floor(smoothed * _Steps) / _Steps;
    
    // 3. 添加边缘检测
    float edge = detect_edge(uv);
    smoothed = lerp(smoothed, _EdgeColor, edge);
    
    return smoothed;
}

2. 油画效果

桑原滤波可以模拟油画笔触的效果：

3. 图像处理 - 去噪

对于含有高频噪声但需要保持边缘的图像，桑原滤波是一个很好的选择：

噪声类型	桑原滤波效果
高斯噪声	⭐⭐⭐
椒盐噪声	⭐⭐⭐⭐
胶片颗粒	⭐⭐⭐⭐

4. 科学可视化

在医学影像、卫星图像等领域，需要平滑区域同时保持组织边界：

• MRI 图像去噪
• 卫星图像分割预处理
• 显微镜图像处理

---

七、常见问题与解决方案

1. 方块 artifacts

问题：滤波结果出现明显的方块状痕迹

原因：四个矩形区域的硬边界

解决方案：

// 使用加权平均替代硬选择
// 根据方差的 softmax 权重混合四个区域的结果
vec3 kuwahara_soft(sampler2D tex, vec2 uv, float radius, vec2 tex_size) {
    vec3 means[4];
    float vars[4];
    
    // ... 计算四个区域的均值和方差 ...
    
    // 计算 softmax 权重
    float sum_exp = 0.0;
    for (int i = 0; i < 4; i++) {
        vars[i] = exp(-vars[i] * _Sharpness);
        sum_exp += vars[i];
    }
    
    // 加权混合
    vec3 result = vec3(0.0);
    for (int i = 0; i < 4; i++) {
        result += means[i] * (vars[i] / sum_exp);
    }
    
    return result;
}

2. 计算太慢

问题：大半径桑原滤波在移动端性能太差

解决方案：

1. 使用降采样（推荐）
2. 降低半径到 2 或 3
3. 使用简化版本（只计算亮度通道的方差）

3. 色彩偏移

问题：滤波结果出现不自然的色彩偏移

原因：分别对 RGB 三个通道独立计算方差和选择

解决方案：

// 只使用亮度通道来选择区域，颜色使用该区域的原始均值
// 这样可以避免不同通道选择不同区域导致的色彩偏移
float luminance = dot(color, vec3(0.299, 0.587, 0.114));
// 只根据亮度方差来选择区域

---

八、参考资料

原始论文

1. "Edge Preserving Smoothing"
   • Kuwahara, M., Hachimura, K., & Eiho, S. (1976)
   • 桑原滤波的原始论文

推荐阅读

• "Real-Time Anisotropic Kuwahara Filtering"

  • Kyprianidis et al. (2009)
  • 各向异性桑原滤波的经典论文

• "Image and Video Abstraction by Anisotropic Kuwahara Filtering"

  • Kyprianidis et al. (2010)
  • 非常好的综述和扩展

开源实现参考

• Unity Kuwahara Filter：https://github.com/keijiro/Kuwahara
• FFmpeg Kuwahara：FFmpeg 内置的桑原滤波实现
• OpenCV：cv.ximgproc.createAnisotropicDiffusion() 类似思想

---

九、总结

桑原滤波是一个非常优雅的算法，用简单的思想实现了出色的边缘保持效果：

核心要点回顾

1. 非线性：不是简单的加权平均，而是根据局部内容自适应选择
2. 边缘保持：通过方差选择平滑区域，天然避免边缘模糊
3. 艺术化：产生的效果有油画质感，非常适合风格化渲染
4. 计算密集：O(n²) 的时间复杂度，需要注意性能优化

什么时候选择桑原滤波？

✅ 适合：

• 卡通渲染、风格化渲染
• 需要艺术化效果的后处理
• 边缘重要的去噪任务

❌ 不适合：

• 通用的图像平滑（用高斯）
• 对性能要求极高的场景
• 需要精确数值的科学计算

桑原滤波虽然是一个"古老"的算法，但它的思想至今仍然闪耀着智慧的光芒。理解它不仅能让你的渲染工具箱增加一件利器，更能让你体会到计算机图形学中"简单即美"的哲学。

---

写在最后

我第一次接触桑原滤波是在研究 NPR 渲染的时候，被它简单却巧妙的思想深深打动。有时候最有效的算法不一定最复杂，而是能够抓住问题本质的那一个。

如果你在实现过程中遇到任何问题，欢迎与我交流！