图像处理之空间滤波

1 原理

1.1 空间滤波简介

  滤波器即只让一部分频率的波形通过来达到波形过滤目的的器件。空间域指一张图像像素平面一定范围内的像素域，相对的是时间域，即多帧图像之间的关系，主要在处理视频帧时描述。在图像处理中，滤波分为两种：

• 频域滤波，将图像转换到频域进行相应的滤波；
• 空间滤波，更准确的说法是掩膜，即在图像的像素平面上利用空间滤波器对图像进行处理。

  空间滤波器是一个$m\times n$大小的矩阵（通常为$3\times 3$，一般都为奇数，偶数核无法完全对齐，需要进行但方向上的padding就会产生图像移位的问题），矩阵定义了希望对当前像素邻域执行的操作和具体的参数，该滤波器会掠过待处理图像的每个像素点来达到，得到对应的输出图像。空间滤波器根据滤波器的特性分为：

• 线性空间滤波：即对邻域内像素的操作为线性。比如当前像素的输出为滤波器与邻域像素点乘积之和；
• 非线性空间滤波：即对邻域内像素的操作为非线性。比如输出像素只取当前像素邻域的最大值。

  比如，对于输入为$M\times N$的图像，滤波器的核大小为$m\times n$的线性滤波的数学表示为：
$$g(x,y)=\sum _{i=-a}^a\sum _{j=-b}^{b}w(i,j)f(x+i,y+j)$$
  其中$g(x,y)$为输出像素，区间$[-a,-b]$的长度为$m$，区间$[-n,n]$的长度为$n$，$w(i,j)$为对应$i,j$位置的滤波器的权重参数。
  在进行空间滤波时如果不对原图像进行填充（padding）的话输出的图像的大小就是$(M-m+1)\times (N-n+1)$，如果希望滤波器的中心能够掠过每个像素或者得到和输入大小箱通风的图像就需要对原图进行填充。填充的方式根据具体的需要来定，最为常用的为在边界上填充固定像素，如果$m\ne n$时，则水平和垂直方向padding的像素值不同。当然也可以在原图像素间填充固定像素值。

1.2 空间相关和卷积

  线性空间滤波中相关和卷积的概念相近：

• 相关：滤波器掠过图像中每个像素计算乘积之和；
• 卷积：首先将滤波器旋转$180^{。}$滤波器略过图像中每个像素计算每个像素的乘积之和。

  二者的区别是，假如输入图像是单位冲激图像，则相关得到的是滤波器输入的180度旋转的值；而卷积得到的是原滤波器的参数值。

2 不同类型的滤波器

2.1 均值滤波器

  均值滤波器，顾名思义，对邻域邻域内像素取均值，属于低通滤波器。均值滤波器可以将预期内的像素值平均化，降低图像中原来比较尖锐的部分抑制图像的高频。均值滤波的数学表达比较简单：
$$g(x,y)=\frac{{\sum }_{s=-a}^a{\sum }_{t=-b}^{b}w(s,t)f(x+s,y+t)}{{\sum }_{s=-a}^a{\sum }_{t=-b}^{b}w(s,t)}$$
  公式中，$w,f,g$分别为滤波器参数，输入像素，输出像素。输出像素的值为领域内每个像素的像素值和滤波器权重逐一相乘然后求和再除以滤波器权重之和求平均。
  影响滤波器的效果有两部分，滤波器的参数和滤波器的大小。滤波器的核越大所能作用的邻域越大，均质化的效果越明显。下图为核尺寸从$[3,33]$步进为2，权重参数为全1的效果。

  如果均值滤波器的权重参数全部为1就是盒装滤波器。下面两个矩阵左边是盒装滤波器的参数，右边是根据像素距离作为权重的矩阵。
$$\left[\begin{array}{ccc}1& 1& 1\\ 1& 1& 1\\ 1& 1& 1\end{array}\right]\left[\begin{array}{ccc}1& 2& 1\\ 2& 4& 2\\ 1& 2& 1\end{array}\right]$$

  在实践中，由于这些模板在一幅图像中跨过的区域很小，即便权重参数不同，通常也很难看出差异。——《数字图像处理》

  另外需要注意的是，从上面的图像中能够明显的看出边缘有黑边，而且滤波器的核越大越明显。原因很简单就是因为边缘填充的是0导致的，实际处理中只需要填充最近的像素值即可。

static Mat avgFilterTransformOneChannel(const Mat &gray, const Mat &weight) {
    int weiHeight = weight.rows, weiWidth = weight.cols;
    assert( gray.channels() == 1 && weight.channels() == 1 && weiHeight % 2 == 1 && weiWidth % 2 == 1);
    double wsum = cv::sum(weight)[0];
    //填充padding像素保证每个像素都能被滤波器掠过
    Mat paddingMat, retMat(cv::Size(gray.rows, gray.cols), CV_8UC1);
    int horPaddPixels = floor(weiHeight / 2), verPaddPixels = floor(weiHeight / 2);
    copyMakeBorder(gray, paddingMat, verPaddPixels, verPaddPixels, horPaddPixels, horPaddPixels, BORDER_CONSTANT, Scalar(0));
    unsigned int imgWidth = gray.cols, imgHeight = gray.rows;
    for (int i = 0; i < imgHeight; i++) {
        for (int j = 0; j < imgWidth; j++) {
            Mat currentRoi = paddingMat(cv::Range(i, i + 2 * horPaddPixels + 1), cv::Range(j, j + 2 * verPaddPixels + 1)), multMat;
            cv::multiply(weight, currentRoi, multMat);
            double currentValue = sum(multMat)[0] / (wsum * 1.0);
            retMat.at<uchar>(i, j) = floor(currentValue);
        }
    }

    return retMat;
}
Mat GrayTransform::avgFilterTransform(const cv::Mat &img, const cv::Mat &weight) {
    std::vector<Mat> chanelMats;
    split(img, chanelMats);
    for (int i = 0; i < chanelMats.size(); i++) {
        chanelMats[i] = avgFilterTransformOneChannel(chanelMats[i], weight);
    }

    Mat returnMat;
    merge(chanelMats, returnMat);
    return returnMat;
}

2.2 统计排序滤波

  统计排序滤波就是利用邻域内像素的统计特性的一种非线性滤波。比如输出的像素取邻域内像素最大值的最大值滤波，取最小值的最小值滤波，取中值的中值滤波。
  下图从左到右分别为：原图，最大值滤波，最小值滤波，中值滤波。

  中值滤波可以用来消除图像中的椒盐噪声。

static double matMaxValue(const Mat &img) {
    assert(img.elemSize() == sizeof(uchar));
    std::vector<uchar> vec(img.begin<uchar>(), img.end<uchar>());
    std::sort(vec.begin(), vec.end());
    return vec.at(vec.size() - 1);
}

static double matMinValue(const Mat &img) {
    assert(img.elemSize() == sizeof(uchar));
    std::vector<uchar> vec(img.begin<uchar>(), img.end<uchar>());
    std::sort(vec.begin(), vec.end());
    return vec.at(0);
}

static double matMedValue(const Mat &img) {
    assert(img.elemSize() == sizeof(uchar));
    std::vector<uchar> vec(img.begin<uchar>(), img.end<uchar>());
    std::sort(vec.begin(), vec.end());
    return vec.size() / 2 == 1 ? vec[int(vec.size() / 2)] : (vec[int(vec.size() / 2)] + vec[int(vec.size() / 2) + 1]) / 2.0;
}

using StaticFilterFunc = double (*)(const Mat &area);
static Mat statsticFilterTransformOneChannel(const cv::Mat gray, StaticFilterFunc func, const Size &filterSize) {
    int weiHeight = filterSize.height, weiWidth = filterSize.width;
    assert(gray.channels() == 1 && weiHeight % 2 == 1 && weiWidth % 2 == 1);
    //填充padding像素保证每个像素都能被滤波器掠过
    Mat paddingMat, retMat(cv::Size(gray.rows, gray.cols), CV_8UC1);
    int horPaddPixels = floor(weiHeight / 2), verPaddPixels = floor(weiHeight / 2);
    copyMakeBorder(gray, paddingMat, verPaddPixels, verPaddPixels, horPaddPixels, horPaddPixels, BORDER_CONSTANT, Scalar(0));
    unsigned int imgWidth = gray.cols, imgHeight = gray.rows;
    for (int i = 0; i < imgHeight; i++) {
        for (int j = 0; j < imgWidth; j++) {
            Mat currentRoi = paddingMat(cv::Range(i, i + 2 * horPaddPixels + 1), cv::Range(j, j + 2 * verPaddPixels + 1)), multMat;
            double currentValue = func(currentRoi);
            retMat.at<uchar>(i, j) = floor(currentValue);
        }
    }

    return retMat;
}

Mat GrayTransform::statsticFilterTransform(const cv::Mat &img, const FilterStatsticType type, const Size &filterSize) {
    std::vector<Mat> chanelMats;
    split(img, chanelMats);
    StaticFilterFunc func;
    switch (type)
    {
    case FILTER_STATSTIC_MIN:
        func = matMinValue;
        break;
    case FILTER_STATSTIC_MAX:
        func = matMaxValue;
        break;
    case FILTER_STATSTIC_MED:
        func = matMedValue;
        break;
    default:
        abort();
    }

    for (int i = 0; i < chanelMats.size(); i++) {
        chanelMats[i] = statsticFilterTransformOneChannel(chanelMats[i], func, filterSize);
    }

    Mat returnMat;
    merge(chanelMats, returnMat);
    return returnMat;
}

2.3 锐化空间滤波器

  如果将卷积相关看做是对邻域内像素和权重进行积分的话，那么相对的说对应像素的微分也可以用用于滤波处理。
  在数学上，微分最直接的含义就是$\frac{\mathrm{d}f}{\mathrm{d}x}={\lim }_{x\to 0}\frac{f(x+\sigma )-f(x)}{\sigma }$。但是在数字图像处理中数据都是离散的，数字量都是有限的，且变换都发生在两个像素之间，因此根据一阶导数的定义，图像上近似的一阶偏微分为：
$$\begin{array}{c}\begin{array}{rl}\frac{\partial f}{\partial x}& =f(x+1,y)-f(x,y)\\ \frac{\partial f}{\partial y}& =f(x,y+1)-f(x,y)\end{array}\end{array}$$

  二阶偏导，我们很容易能够推导出：

$$\begin{array}{c}\begin{array}{rl}\frac{{\partial }^{2}f}{\partial x^2}& =f(x+1,y)+f(x-1)-2f(x,y)\\ \frac{{\partial }^{2}f}{\partial y^2}& =f(x,y+1)+f(x,y-1)-2f(x,y)\end{array}\end{array}$$

  数字图像中，对于有边缘的区域，一阶导数很容易形成一个比较粗略的边缘，而二阶导数可以产生由0分开的一个像素宽的双边缘。因此，二阶微分在增强图形的细节上要比一阶微分要好。下图原图在为$x$和$y$方向上的一阶微分和二阶微分结果，可以明显的看到一阶微分的边缘要比二阶微分宽（从左到右分别为原图，x方向上的一阶微分，x方向上的二阶微分，y方向上的一阶微分，y方向上的二阶微分）

2.3.1 拉普拉斯算子

  拉普拉斯算子是一个各向同性算子，滤波器的响应和滤波器作用的图像的方向无关，即即便图像旋转后进行滤波和滤波后旋转效果相同。拉普拉斯算子的定义为：
$${▽}^{2}f=\frac{{\partial }^{2}f}{\partial x^2}+\frac{{\partial }^{2}f}{\partial y^2}$$
  根据上面求得的二阶微分，将其应用到离散图像上：
$${▽}^{2}f(x,y)=f(x+1,y)+f(x-1,y)+f(x,y-1)+f(x,y+1)-4f(x,y)$$
  上面的公式描述了以$90^{。}$为增量进行旋转的各向同性拉普拉斯算子，如果希望扩展到$45^{。}$只需要将45度方向上的二阶偏导引入即可：
$${▽}^{2}f(x,y)=f(x+1,y)+f(x-1,y)+f(x,y-1)+f(x,y+1)+f(x-1,y-1)+f(x-1,y+1)+f(x+1,y-1)+f(x+1,y+1)-8f(x,y)$$
  应用到空间滤波中，对应的矩阵参数分别为
$$\left[\begin{array}{ccc}0& 1& 0\\ 1& -4& 1\\ 0& 1& 0\end{array}\right]\left[\begin{array}{ccc}1& 1& 1\\ 1& -8& 1\\ 1& 1& 1\end{array}\right]$$
  通过上面的方式我们能够得图像的边缘，直接将计算出的结果通过下面的公式和原图相加即可。
$$g(x,y)=f(x,y)+c{▽}^{2}f(x,y)$$
  下图中展示了两种拉普拉斯空间滤波的结果，上一张图是90度各向同性的结果，下面的图示45度各向同性的结果。（下图中从左到右为1,2,3,4四张图，分别为原图，拉普拉斯变化值，标定后的拉普拉斯变化值，和原图线性相加后的值c=-1）

  下图为参数c对两种矩阵的影响。

//标定，实际上就是将计算出来的值映射到0-255
static Mat stardandOneChannel(const cv::Mat &gray) {
    Mat fltMat;
    Mat returnMat(gray.rows, gray.cols, CV_32FC1);
    constexpr int layer = 255;
    gray.convertTo(fltMat, CV_64FC1);
    double maxValue = matMaxValue<double>(fltMat);
    double minValue = matMinValue<double>(fltMat);
    for (int i = 0; i < fltMat.rows; i++) {
        for (int j = 0; j < fltMat.cols; j++) {
            double value = (fltMat.at<double>(i, j) - minValue * 1.0) / (maxValue - minValue) * layer;
            returnMat.at<float>(i, j) = static_cast<float>(value);
        }
    }

    return returnMat;
}

Mat GrayTransform::standardMat(const cv::Mat &img) {
    std::vector<Mat> chanelMats;
    split(img, chanelMats);
    for (int i = 0; i < chanelMats.size(); i++) {
        chanelMats[i] = stardandOneChannel(chanelMats[i]);
    }

    Mat returnMat;
    merge(chanelMats, returnMat);
    return returnMat;
}

2.3.2 非锐化掩蔽和高提升滤波

  非锐化掩蔽是利用模糊图像和原图之间的差值模板来突出原图中的边缘，然后将得到的模板加到原图来突出边缘达到锐化的目的。
$$\begin{array}{c}\begin{array}{rl}{g}_{mask}(x,y)& =f(x,y)-\overline{f}(x,y)\\ g(x,y)& =f(x,y)+k{g}_{mask}(x,y)\end{array}\end{array}$$
  其中$\overline{f}(x,y)$为经过空间滤波的模糊图像，$k$用来控制最终希望得到的图像的锐化程度。

• $k<1$，则不强调非锐化模板的影响；
• $k=1$，为非锐化掩蔽；
• $k>1$，为高提升滤波。

  下图从左到右分别为：原图，均值模糊，差值，k=0.5，k=1，k=2。

2.3.3 一阶微分图像锐化

  对于一张图像有多个方向，一般很容易求得$x，y$两个方向的偏导$g_x=\frac{\partial f}{\partial x}$和g_y=$\frac{\partial f}{\partial y}$，该偏导只表示了当前点在对应方向上的变化程度。而值$M(x,y)=mag(▽)=\sqrt{g_x^2+g_y^2}$表示了该点的梯度。再具体实现中可用$M(x,y)=|g_x|+|g_y|$求得梯度的近似。
  罗伯特交叉梯度算子利用当前点的对角值计算梯度，能够比直接在x或者y方向上计算梯度多近似0.5个像素的偏移，且$g_x,g_y$都包含了x和y方向上的梯度信息。罗伯特交叉梯度算子中：
$$\begin{array}{c}\begin{array}{rl}{g}_x(x,y)& =f(x+1,y+1)-f(x,y)\\ g_y(x,y)& =f(x,y+1)-f(x+1,y)\\ M(x,y)& \approx |f(x+1,y+1)-f(x,y)|+|f(x,y+1)-f(x+1,y)|\end{array}\end{array}$$
  对应的算子矩阵为
$$\left[\begin{array}{cc}-1& 0\\ 0& 1\end{array}\right]\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right]$$

  罗伯特交叉梯度算子为偶数算子没有对称中心，不容易实现。Soble算子采用一阶微分，奇数核。其中第二个权重系数为2是为了突出中心点来平滑图像。
$$\begin{array}{c}\begin{array}{rl}{g}_x(x,y)& =(f(x+1,y-1)+2f(x+1,y)+f(x+1,y-1))-(f(x-1,y-1)+2f(x-1,y)+f(x-1,y-1))\\ g_y(x,y)& =(f(x+1,y+1)+2f(x,y+1)+f(x-1,y+1))-(f(x+1,y-1)+2f(x,y-1)+f(x-1,y-1))\\ M(x,y)& \approx |g_x|+|g_y|\end{array}\end{array}$$
  对应的矩阵为
$$g_x=\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right]g_y=\left[\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right]$$

  下图分别为原图，$g_x$,$g_y$,$0.3M(x,y)$,处理后的图。