3. 双边滤波保持边缘的ToF深度图去噪

任务说明

ToF（Time-of-Flight）深度相机在弱反射、多路径干扰等场景下，深度图常含有高斯噪声与离群脉冲噪声，同时物体边缘处深度跳变明显。本案例实现一种保持边缘的深度图去噪算法：采用双边滤波（Bilateral Filter）对深度图进行处理，在平滑噪声的同时保留深度不连续边缘。输入为单张16位深度图（单位mm），输出为滤波后的16位深度图，并给出可量化的边缘保持评价指标。

算法原理与实现

双边滤波的权重由空间邻近度（高斯核）与灰度/深度相似度（高斯核）共同决定：

I'(p) = (1/W_p) · Σ_{q∈N(p)} I(q) · exp(-||p-q||²/(2σ_s²)) · exp(-|I(p)-I(q)|²/(2σ_r²))

其中σ_s控制空间平滑范围，σ_r控制深度值差异的敏感度。当像素位于边缘两侧时，深度差异大，对应权重趋近于0，从而避免边缘被模糊。本实现采用分离式加速：先计算空间权重（仅与坐标差有关，可预计算），再对每个像素遍历邻域累加。为处理深度图中的无效像素（值为0），滤波时跳过这些点并重新归一化权重。

代码实现


#include <opencv2/opencv.hpp>
#include <vector>
#include <cmath>
#include <iostream>

// 保持边缘的双边滤波（针对16位深度图）
// src: 输入CV_16UC1深度图（单位mm，0表示无效）
// dst: 输出滤波后深度图
// sigmaSpace: 空间高斯标准差（像素）
// sigmaRange: 深度值高斯标准差（mm）
// radius: 邻域半径（实际窗口为2*radius+1）
void bilateralDepthFilter(const cv::Mat& src, cv::Mat& dst,
                          double sigmaSpace, double sigmaRange, int radius) {
    CV_Assert(src.type() == CV_16UC1);
    dst = cv::Mat::zeros(src.size(), CV_16UC1);

    // 预计算空间权重表
    int winSize = 2 * radius + 1;
    std::vector<float> spaceWeight(winSize * winSize);
    float spaceCoeff = -0.5f / (sigmaSpace * sigmaSpace);
    for (int dy = -radius; dy <= radius; ++dy) {
        for (int dx = -radius; dx <= radius; ++dx) {
            int idx = (dy + radius) * winSize + (dx + radius);
            spaceWeight[idx] = std::exp(spaceCoeff * (dx * dx + dy * dy));
        }
    }

    float rangeCoeff = -0.5f / (sigmaRange * sigmaRange);
    int rows = src.rows, cols = src.cols;

    for (int y = 0; y < rows; ++y) {
        const uint16_t* srcRow = src.ptr<uint16_t>(y);
        uint16_t* dstRow = dst.ptr<uint16_t>(y);
        for (int x = 0; x < cols; ++x) {
            uint16_t centerVal = srcRow[x];
            if (centerVal == 0) continue; // 无效像素保持0

            double sum = 0.0, weightSum = 0.0;
            for (int dy = -radius; dy <= radius; ++dy) {
                int ny = y + dy;
                if (ny < 0 || ny >= rows) continue;
                const uint16_t* neighborRow = src.ptr<uint16_t>(ny);
                for (int dx = -radius; dx <= radius; ++dx) {
                    int nx = x + dx;
                    if (nx < 0 || nx >= cols) continue;
                    uint16_t neighborVal = neighborRow[nx];
                    if (neighborVal == 0) continue; // 跳过无效点

                    int idx = (dy + radius) * winSize + (dx + radius);
                    float w = spaceWeight[idx];
                    if (w == 0.0f) continue;

                    int diff = (int)centerVal - (int)neighborVal;
                    w *= std::exp(rangeCoeff * diff * diff);
                    sum += w * neighborVal;
                    weightSum += w;
                }
            }
            if (weightSum > 0) {
                dstRow[x] = (uint16_t)std::round(sum / weightSum);
            } else {
                dstRow[x] = centerVal; // 邻域全无效时保留原值
            }
        }
    }
}

// 计算边缘保持率：边缘像素中滤波前后深度差小于阈值的比例
double edgePreservationRate(const cv::Mat& src, const cv::Mat& dst, int edgeThreshold = 50) {
    cv::Mat srcEdge, dstEdge;
    cv::Sobel(src, srcEdge, CV_32F, 1, 0, 3);
    cv::Sobel(dst, dstEdge, CV_32F, 1, 0, 3);
    cv::Mat srcMag, dstMag;
    cv::magnitude(srcEdge, srcEdge, srcMag);
    cv::magnitude(dstEdge, dstEdge, dstMag);

    // 边缘像素：原图梯度幅值大于阈值
    cv::Mat edgeMask = srcMag > edgeThreshold;
    int edgeCount = cv::countNonZero(edgeMask);
    if (edgeCount == 0) return 1.0;

    // 边缘处深度变化量
    cv::Mat diff;
    cv::absdiff(src, dst, diff);
    cv::Mat edgeDiff;
    diff.copyTo(edgeDiff, edgeMask);
    int preserved = cv::countNonZero(edgeDiff <= 30); // 30mm内视为保持
    return (double)preserved / edgeCount;
}

int main() {
    // 生成模拟深度图：200x200，含两个深度平面（1000mm和2000mm）及高斯噪声
    cv::Mat depth(200, 200, CV_16UC1, cv::Scalar(1000));
    for (int y = 0; y < 200; ++y) {
        for (int x = 0; x < 200; ++x) {
            if (x > 100) depth.at<uint16_t>(y, x) = 2000;
        }
    }
    // 添加高斯噪声（sigma=15mm）
    cv::Mat noise = cv::Mat::zeros(depth.size(), CV_16SC1);
    cv::randn(noise, 0, 15);
    cv::Mat depthNoisy;
    depth.convertTo(depthNoisy, CV_16SC1);
    depthNoisy += noise;
    depthNoisy.convertTo(depth, CV_16UC1);

    // 执行双边滤波
    cv::Mat filtered;
    bilateralDepthFilter(depth, filtered, 5.0, 30.0, 3);

    // 验证：计算边缘保持率与平滑度
    double edgeRate = edgePreservationRate(depth, filtered);
    cv::Mat smoothRegion = depth(cv::Rect(20, 20, 60, 60)); // 左平面平滑区
    cv::Mat smoothFiltered = filtered(cv::Rect(20, 20, 60, 60));
    cv::Scalar meanSrc, stdSrc, meanDst, stdDst;
    cv::meanStdDev(smoothRegion, meanSrc, stdSrc);
    cv::meanStdDev(smoothFiltered, meanDst, stdDst);

    std::cout << "边缘保持率: " << edgeRate << std::endl;
    std::cout << "平滑区噪声标准差(滤波前): " << stdSrc[0] << " mm" << std::endl;
    std::cout << "平滑区噪声标准差(滤波后): " << stdDst[0] << " mm" << std::endl;

    // 保存结果
    cv::imwrite("depth_noisy.png", depth);
    cv::imwrite("depth_filtered.png", filtered);
    return 0;
}

编译与运行


# 依赖：OpenCV 4.x（含core、imgproc、imgcodecs模块）
g++ -o bilateral_depth bilateral_depth.cpp $(pkg-config --cflags --libs opencv4)
./bilateral_depth

验证方法与判定标准

程序输出两个关键指标：

边缘保持率：应大于0.85（即85%以上的边缘像素深度变化不超过30mm），证明边缘未被过度平滑；
平滑区噪声标准差：滤波后应显著低于滤波前（本模拟场景中滤波前约15mm，滤波后应降至5mm以下），证明噪声被有效抑制。

同时可打开输出的depth_noisy.png与depth_filtered.png，目视检查垂直边缘是否清晰锐利，而平坦区域是否更均匀。若边缘保持率偏低，可增大σ_r；若噪声残留较多，可增大σ_s或radius。

📚 课程目录