为图像处理初学者设计的 100 个问题 问题二十一至问题三十

问题二十一至问题三十

关于直方图的几个操作的中文翻译一直都十分混乱(成龙抓头.jpg)。下面是我查了资料做的一个对照表。我自己也不是很确定,仅供参考

中文 English 日本語 具体作用
直方图匹配(规定化) Histogram Matching ( Specification ) ヒストグラムマッチング 将两张图像的累积分布函数调为一样
直方图均衡化 Histogram Equalization ヒストグラム均等化(平坦化) 拉开灰度差别,增强对比度
直方图归一化 Histogram Normalization ヒストグラム正規化 将直方图所有分量限制在一定范围

问题二十一:直方图归一化( Histogram Normalization )

来归一化直方图吧。

有时直方图会偏向一边。

比如说,数据集中在$0$处(左侧)的图像全体会偏暗,数据集中在$255$处(右侧)的图像会偏亮。

如果直方图有所偏向,那么其动态范围( dynamic range )就会较低。

为了使人能更清楚地看见图片,让直方图归一化、平坦化是十分必要的。

这种归一化直方图的操作被称作灰度变换(Grayscale Transformation)。像素点取值范围从$[c,d]$转换到$[a,b]$的过程由下式定义。这回我们将MARKDOWN_HASHcddf1f424e6ef89f3c7af3e98cedc79aMARKDOWNHASH的灰度扩展到$[0, 255]$范围:
$$
x
{out}=
\begin{cases}
a& (\text{if}\quad x{in}<c)\
\frac{b-a}{d-c}\ (x
{in}-c)+a&(\text{else if}\quad c\leq x_{in}<d)\
b&(\text{else})
\end{cases}
$$

Answer

void HistogramNormalization(Mat& src,int start=0,int end=255) {
    if (!src.data) return;
    Mat chs[3];
    split(src, chs);
    double maxItem[3], minItem[3];
    for(int i=0;i<3;i++)
        minMaxLoc(chs[i],&minItem[i],&maxItem[i],0,0);
    src.forEach<Vec3b>([&](Vec3b& pix, const int *position) {
        int channal = src.channels();
        for (int c = 0; c < channal; c++) {
            if (pix[c]<minItem[c]) pix[c]=start;
            else if(pix[c] > maxItem[c]) pix[c] = end;
            else pix[c] = ((end - start) / (maxItem[c] - minItem[c])*pix[c] + start);
        }
    });
}

Show

image_1e51nibifn6omb1tg5te018ks9.png-26.4kBimage_1e51nio33f5k1482cft1k4q4iqm.png-39kB
以及两个直方图比较
image_1e51nqrj41cta1vbq1fi6shufnt13.png-32.8kB

Note

  • 在一个Mat中寻找最大值最小值可以使用minMaxLoc()

问题二十二:直方图操作

让直方图的平均值$m_0=128$,标准差$s_0=52$​吧!

这里并不是变更直方图的动态范围,而是让直方图变得平坦。

可以使用下式将平均值为$m$标准差为$s$的直方图变成平均值为$m_0$标准差为$s0$的直方图:
$$
x
{out}=\frac{s0}{s}\ (x{in}-m)+m_0
$$

Answer

void modifyMeanAndStandardDeviation(Mat& src,double mean, double sigma) {
    Mat chs[3];
    split(src, chs);
    Mat meanItem[3], devItem[3];
    for (int i = 0; i < 3; i++)
        meanStdDev(chs[i], meanItem[i], devItem[i]);
    src.forEach<Vec3b>([&](Vec3b& pix, const int *position) {
        int channal = src.channels();
        for (int c = 0; c < channal; c++) {
            pix[c] = (sigma / devItem[c].at<double>(0, 0))*(pix[c] - meanItem[c].at<double>(0, 0)) + mean;
        }
    });
}

Show

image_1e51pmcmq1qs51etcu1inpmpr1g.png-67.5kB
image_1e51pn4aaa991r7q109c1igp5ll1t.png-31.6kB

Note

求Mat的均值可以使用如下代码求

Mat mean,dev;
meanStdDev(src,mean,dev);
meam.at<double>(0,0);//均值
dev.at<double>(0,0);//标准差

问题二十三:直方图均衡化( Histogram Equalization )

来让均匀化直方图吧!

直方图均衡化是使直方图变得平坦的操作,是不需要计算上面的问题中的平均值、标准差等数据使直方图的值变得均衡的操作。

均衡化操作由以下式子定义。$S$是总的像素数;$Z{max}$是像素点的最大取值(在这里是$255$);$h(z)$表示取值为$z$的累积分布函数:
$$
Z' = \frac{Z
{max}}{S} \ \sum\limits_{i=0}^z\ h(i)
$$

Answer

void HistogramEqualization(Mat& src, uchar maxValue = 255) {
    int gross = src.channels()*src.cols*src.rows;
    double hist[255] = { 0.0 };
    int channal = src.channels();
    src.forEach<Vec3b>([&](Vec3b& pix, const int* position) {
        for (int c = 0; c < channal; c++){
            hist[pix[c]] += 1;
        }
    });
    src.forEach<Vec3b>([&](Vec3b& pix, const int * position) {
        for (int c = 0; c < channal; c++) {
            int cur_vaule = pix[c];
            int hist_sum = 0;
            for (int i = 0; i < cur_vaule; i++) hist_sum += hist[i];
            pix[c] = (int)((float)maxValue *(float)hist_sum / (float)gross);
        }
    });
}

Show

image_1e55b5alkp1b92qlhreaqm349.png-31.4kBimage_1e55b5svr2h915c11j571gn21unm.png-44.8kB
image_1e55b6kejtj75fu1tno5fpqg13.png-32.2kB

Note

  • $\sum\limits_{i=0}^z\ h(i)$可以提前在循环外面计算好,通过查表的方式得到结果可以有效提升计算速度

问题二十四:伽玛校正(Gamma Correction)1

imori_gamma.jpg进行伽马校正($c=1$,$g=2.2$)吧!

伽马校正用来对照相机等电子设备传感器的非线性光电转换特性进行校正。如果图像原样显示在显示器等上,画面就会显得很暗。伽马校正通过预先增大 RGB 的值来排除显示器的影响,达到对图像修正的目的。

由于下式引起非线性变换,在该式中,$x$被归一化,限定在$[0,1]$范围内。$c$是常数,$g$为伽马变量(通常取$2.2$):
$$
x' = c\ {I{in}}^ g
$$
因此,使用下面的式子进行伽马校正:
$$
I
{out} ={\frac{1}{c}\ I_{in}} ^ {\frac{1}{g}}
$$

Principle

  • Why
    人眼对外界光源的感光值与输入光强不是呈线性关系的,而是呈指数型关系的。在低照度下,人眼更容易分辨出亮度的变化,随着照度的增加,人眼不易分辨出亮度的变化。
    而摄像机感光与输入光强呈线性关系。
  • How
    假设图像中有一个像素,值是 200 ,那么对这个像素进行校正必须执行如下步骤: 
      1. 归一化 :将像素值转换为  0 ~ 1  之间的实数。 算法如下 : ( i + 0. 5)/256  这里包含 1 个除法和 1 个加法操作。对于像素  A  而言  , 其对应的归一化值为  0. 783203 。 

  2. 预补偿 :根据公式  , 求出像素归一化后的 数据以  1 /gamma  为指数的对应值。这一步包含一个 求指数运算。若  gamma  值为  2. 2 ,  则  1 /gamma  为  0. 454545 , 对归一化后的  A  值进行预补偿的结果就 是  0. 783203 ^0. 454545 = 0. 894872 。 

  3. 反归一化 :将经过预补偿的实数值反变换为  0  ~  255  之间的整数值。具体算法为 : f*256 - 0. 5  此步骤包含一个乘法和一个减法运算。续前 例  , 将  A  的预补偿结果  0. 894872  代入上式  , 得到  A  预补偿后对应的像素值为  228 , 这个  228  就是最后送 入显示器的数据。

Answer

void GammaCorrection(Mat& src, int coefficient = 1, double gama = 2.2) {
    int channel = src.channels();

    double normalizationMap[256];//通过查表的方式归一化
    for (int i = 0; i < 256; i++) {
        normalizationMap[i] = (i + 0.5) / 256;
    }

    //补偿
    src.forEach<Vec3b>([&](Vec3b& pix, const int * position) {
        for (int c = 0; c < channel; c++) {
            double after_normal = normalizationMap[pix[c]];//归一化结果
            double effect = coefficient*pow(after_normal, 1/gama);
            pix[c] = (uchar)(effect * 256 - 0.5);
        }
    });

}

Show

image_1e55dv3jm166j1k7412jh17051l591g.png-46.6kBimage_1e55dveha1l5p11nkdec2m919pa1t.png-48.7kB
image_1e55e06bc1ct1198j1603k5k3p92a.png-32.6kB

问题二十五:最邻近插值( Nearest-neighbor Interpolation )

使用最邻近插值将图像放大$1.5$倍吧!

最近邻插值在图像放大时补充的像素取最临近的像素的值。由于方法简单,所以处理速度很快,但是放大图像画质劣化明显。

使用下面的公式放大图像吧!$I'$为放大后图像,$I$为放大前图像,$a$为放大率,方括号是四舍五入取整操作:

image_1e55ee8kbhdp1hbv1s1d1vgohmh9.png-307.7kB
$$
I'(x,y) = I([\frac{x}{a}], [\frac{y}{a}])
$$

Answer

//最邻近插值的方法放大图像
void Nearest_neighborInterpolation(Mat&src, double times) {
    Mat result = Mat::zeros(src.rows*times, src.cols*times, CV_8UC3);
    result.forEach<Vec3b>([&](Vec3b& pix, const int * position) {
        int x = position[0];
        int y = position[1];
        int target_x = (double)x / (times);
        int target_y = (double)y / (times);
        for (int c = 0; c < 3; c++) {
            pix[c] = src.at<Vec3b>(target_x, target_y)[c];
        }
    });
    src=result;
}

Show

image_1e55f2bkm1o5iq9qk5tv1919ejm.png-103.9kB

Note

放大后有明显马赛克痕迹,可以尝试加入之前习题中的滤波器
以下为放大6倍并均值滤波效果
image_1e55fkq9r1besh1q18pe2nnr8713.png-832.8kB

问题二十六:双线性插值( Bilinear Interpolation )

使用双线性插值将图像放大$1.5$倍吧!

双线性插值考察$4$邻域的像素点,并根据距离设置权值。虽然计算量增大使得处理时间变长,但是可以有效抑制画质劣化。

  1. 放大后图像的座标$(x',y')$除以放大率$a$,可以得到对应原图像的座标$(\lfloor \frac{x'}{a}\rfloor , \lfloor \frac{y'}{a}\rfloor)$。

  2. 求原图像的座标$(\lfloor \frac{x'}{a}\rfloor , \lfloor \frac{y'}{a}\rfloor)$周围$4$邻域的座标$I(x,y)$,$I(x+1,y)$,$I(x,y+1)$,$I(x+1, y+1)$:

    image_1e56hv0thc9a83hak37pnvqj1g.png-348.7kB

  3. 分别求这4个点与$(\frac{x'}{a}, \frac{y'}{a})$的距离,根据距离设置权重:$w = \frac{d}{\sum\ d}$

  4. 根据下式求得放大后图像$(x',y')$处的像素值:
    $$
    d_x = \frac{x'}{a} - x\
    d_y = \frac{y'}{a} - y\
    I'(x',y') = (1-d_x)\ (1-d_y)\ I(x,y) + d_x\ (1-d_y)\ I(x+1,y) + (1-d_x)\ d_y\ I(x,y+1) + d_x\ d_y\ I(x+1,y+1)
    $$

Answer

//双线性差值
void BilinearInterpolation(Mat& src, double times) {
    Mat result = Mat::zeros(src.rows*times, src.cols*times, CV_8UC3);
    result.forEach<Vec3b>([&](Vec3b& pix, const int * position) {
        int x = position[0];
        int y = position[1];
        int y_before = (float)y / times;
        y_before = std::min(y_before, src.cols - 1);
        int x_before = (float)x / times;
        x_before = std::min(x_before, src.rows-1);
        double dx = x / times - x_before;
        double dy = y / times - y_before;
        for (int c = 0; c < 3; c++) {
            double val = (1 - dx)*(1 - dy)*src.at<Vec3b>(x_before, y_before)[c] +
                dx * (1 - dy)*src.at<Vec3b>(x_before + 1, y_before)[c] +
                (1 - dx)*(dy)*src.at<Vec3b>(x_before, y_before + 1)[c] +
                dx * dy*src.at<Vec3b>(x_before + 1, y_before + 1)[c];
            pix[c] = (uchar)val;
        }
    });
    src = result;
}

Show

这里是放大了6倍
image_1e57mnthm17ualdi1s701luj1c7d9.png-847kB

问题二十七:双三次插值( Bicubic Interpolation )

使用双三次插值将图像放大$1.5$倍吧!

双三次插值是双线性插值的扩展,使用邻域$16$像素进行插值。

各自像素间的距离由下式决定:
$$
\begin{align}
d_{x1} = |\frac{x'}{a\ x} - (x-1)|\quad
d
{x2} = |\frac{x'}{a\ x}- x| \quad
d
{x3} = |\frac{x'}{a\ x}- (x+1)|\quad
d
{x4} = |\frac{x'}{a\ x} - (x+2)|\
d
{y1} = |\frac{x'}{a\ y} - (y-1)|\quad
d
{y2} = |\frac{x'}{a\ y} - y| \quad
d
{y3} = |\frac{x'}{a\ y} - (y+1)| \quad
d
{y_4} = |\frac{x'}{a\ y} - (y+2)|
\end{align
}
$$
权重由基于距离的函数取得。$a$在大部分时候取$-1$。大体上说,图中蓝色像素的距离$|t|\leq 1$,绿色像素的距离$1<|t|\leq 2$:
$$
h(t)=
\begin{cases}
(a+2)\ |t|^3 - (a+3)\ |t|^2 + 1 &\text{when}\quad |t|\leq 1 \
a\ |t|^3 - 5\ a\ |t|^2 + 8\ a\ |t| - 4\ a&\text{when}\quad 1<|t|\leq 2\
0&\text{else}
\end{cases}
$$
利用上面得到的权重,通过下面的式子扩大图像。将每个像素与权重的乘积之和除以权重的和。
$$
I'(x', y')=\frac{1}{\sum\limits{j=1}^4\ \sum\limits{i=1}^4\ h(d{xi})\ h(d{yj})}\ \sum\limits{j=1}^4\ \sum\limits{i=1}^4\ I(x+i-2,y+j-2)\ h(d{xi})\ h(d{yj})
$$

问题二十八:仿射变换( Afine Transformations )——平行移动

以下公式中有些矩阵乘法左右是反的,原文中便是如此,这里未做更改,无法相乘的部分需要大家自己做纠正.

利用仿射变换让图像在$x$方向上$+30$,在$y$方向上$-30$吧!

仿射变换利用$3\times3$的矩阵来进行图像变换。

变换的方式有平行移动(问题二十八)、放大缩小(问题二十九)、旋转(问题三十)、倾斜(问题三十一)等。

原图像记为$(x,y)$,变换后的图像记为$(x',y')$。

图像放大缩小矩阵为下式:
$$
\left(
\begin{matrix}
x'\
y'
\end{matrix}
\right)=
\left(
\begin{matrix}
a&b\
c&d
\end{matrix}
\right)\
\left(
\begin{matrix}
x\
y
\end{matrix}
\right)
$$
另一方面,平行移动按照下面的式子计算:
$$
\left(
\begin{matrix}
x'\
y'
\end{matrix}
\right)=
\left(
\begin{matrix}
x\
y
\end{matrix}
\right)+
\left(
\begin{matrix}
t_x\
t_y
\end{matrix}
\right)
$$
把上面两个式子盘成一个:
$$
\left(
\begin{matrix}
x'\
y'\
1
\end{matrix}
\right)=
\left(
\begin{matrix}
a&b&t_x\
c&d&t_y\
0&0&1
\end{matrix}
\right)\
\left(
\begin{matrix}
x\
y\
1
\end{matrix}
\right)
$$
但是在实际操作的过程中,如果一个一个地计算原图像的像素的话,处理后的像素可能没有在原图像中有对应的坐标。2

因此,我们有必要对处理后的图像中各个像素进行仿射变换逆变换,取得变换后图像中的像素在原图像中的坐标。仿射变换的逆变换如下:
$$
\left(
\begin{matrix}
x\
y
\end{matrix}
\right)=
\frac{1}{a\ d-b\ c}\
\left(
\begin{matrix}
d&-b\
-c&a
\end{matrix}
\right)\
\left(
\begin{matrix}
x'\
y'
\end{matrix}
\right)-
\left(
\begin{matrix}
t_x\
t_y
\end{matrix}
\right)
$$
这回的平行移动操作使用下面的式子计算。$t_x$和$t_y$是像素移动的距离。
$$
\left(
\begin{matrix}
x'\
y'\
1
\end{matrix}
\right)=
\left(
\begin{matrix}
1&0&t_x\
0&1&t_y\
0&0&1
\end{matrix}
\right)\
\left(
\begin{matrix}
x\
y\
1
\end{matrix}
\right)
$$

Answer

先根据放射变换定义工具函数,原理就是根据放射变换的逆变换公式得出目标图像某一位置的像素在原图像中对应的位置

$$
\left(
\begin{matrix}
x\
y
\end{matrix}
\right)=
\frac{1}{a\ d-b\ c}\
\left(
\begin{matrix}
d&-b\
-c&a
\end{matrix}
\right)\
\left(
\begin{matrix}
x'\
y'
\end{matrix}
\right)-
\left(
\begin{matrix}
t_x\
t_y
\end{matrix}
\right)
$$

initSingleChannelMat()在之前的题目中有定义过,这里使用该方式初始化矩阵或者用opencv自带的方式都可以

//仿射变换的逆变换
Mat AfineTransformations(Mat& src, double a, double b, double c, double d, double tx, double ty) {
    double  coefficient = 1 / (a*d - b * c);
    double t[4] = {
        d,-b,
        -c,a
    };
    Mat a_ = initSingleChannelMat<double>(t,2, CV_64FC1);

    double t_[2] = {
        tx,
        ty
    };
    Mat b_(Size(1,2), CV_64FC1, t_);
    src.convertTo(src, CV_64FC1);
    Mat step_1 = a_ * src;
    return coefficient * step_1+ b_;
}

在根据平移需要的仿射变换进行变换

//通过放射变换实现平移
void transAT(Mat& src, double tx, double ty) {
    Mat result = Mat::zeros(src.rows, src.cols, CV_8UC3);
    result.forEach<Vec3b>([&](Vec3b& pix, const int * position) {
        int x = position[0];
        int y = position[1];
        int cur_pos[2][17] = { x,y };
        Mat m(Size(1, 2), CV_32SC1, cur_pos);
        Mat ori_pos = AfineTransformations(m, 1, 0, 0, 1, tx, ty);
        //printMat(m);
        for (int i = 0; i < 3; i++) {
            int ori_pos_x = ori_pos.at<Vec<double, 1>>(0, 0)[0];
            int ori_pos_y = ori_pos.at<Vec<double, 1>>( 1,0)[0];
            if (ori_pos_x >= src.rows || ori_pos_y >= src.cols || ori_pos_x < 0 || ori_pos_y < 0) {
                pix[i] = 0;
            }
            else
                pix[i] = src.at<Vec3b>(ori_pos_x, ori_pos_y)[i];
        }
    });
    src = result;
}

Show

image_1e59ijg9vs6l1j7t1a1d419lf3m.png-45.2kBimage_1e59ij54k165og69f6g4pjgl59.png-33.4kB

Note

  • 在进行矩阵乘法运算时要注意两个Mat 的数据类型,比如说CV_64FC1CV_32SC1之间无法直接相乘,会产生如下错误:
    image_1e59iorg61c6s17tn1mq918ds8os13.png-10.2kB
    可以使用src.convertTo(src, CV_64FC1);函数来进行主句类型之间的转换.其中src是一个Mat对象的实例.

问题二十九:仿射变换( Afine Transformations )——放大缩小

  1. 使用仿射变换,将图片在$x$方向上放大$1.3$倍,在$y$方向上缩小至原来的$\frac{4}{5}$。
  2. 在上面的条件下,同时在$x$方向上向右平移$30$($+30$),在$y$方向上向上平移$30$($-30$)。

Answer

//通过放射变换实现放大和缩小
void amplifyAT(Mat& src, double dx, double dy) {
    Mat result = Mat::zeros(src.rows*dx, src.cols*dy, CV_8UC3);

    result.forEach<Vec3b>([&](Vec3b& pix, const int * position) {
        int x = position[0];
        int y = position[1];
        int cur_pos[2][21] = { x,y };
        Mat m(Size(1, 2), CV_32SC1, cur_pos);
        Mat ori_pos = AfineTransformations(m, dx, 0, 0, dy, 0,0);
        //printMat(m);
        for (int i = 0; i < 3; i++) {
            int ori_pos_x = ori_pos.at<Vec<double, 1>>(0, 0)[0];
            int ori_pos_y = ori_pos.at<Vec<double, 1>>(1, 0)[0];
            if (ori_pos_x >= src.rows || ori_pos_y >= src.cols || ori_pos_x < 0 || ori_pos_y < 0) {
                pix[i] = 0;
            }
            else
                pix[i] = src.at<Vec3b>(ori_pos_x, ori_pos_y)[i];
        }
    });
    src = result;
}

Show

该变换是横向变为2倍,纵向变为0.5倍
image_1e59jc0k3pvlcd0bqqij1flh1s.png-73.2kB

问题三十:仿射变换( Afine Transformations )——旋转

  1. 使用仿射变换,逆时针旋转$30$度。
  2. 使用仿射变换,逆时针旋转$30$度并且能让全部图像显现(也就是说,单纯地做仿射变换会让图片边缘丢失,这一步中要让图像的边缘不丢失,需要耗费一些工夫)。

使用下面的式子进行逆时针方向旋转$A$度的仿射变换:

$$
\left(
\begin{matrix}
x'\
y'\
1
\end{matrix}
\right)=
\left(
\begin{matrix}
\cos(A)&-\sin(A)&t_x\
\sin(A)&\cos(A)&t_y\
0&0&1
\end{matrix}
\right)\
\left(
\begin{matrix}
x\
y\
1
\end{matrix}
\right)
$$

Answer

//通过放射变换实现旋转
void rollAT(Mat& src, double angle) {
    Mat result = Mat::zeros(src.rows, src.cols, CV_8UC3);
    angle = 3.1415926*angle /  180;
    result.forEach<Vec3b>([&](Vec3b& pix, const int * position) {
        int x = position[0];
        int y = position[1];
        int cur_pos[2][1] = { x,y };
        Mat m(Size(1,2), CV_32SC1, cur_pos);
        Mat ori_pos = AfineTransformations(m, cos(angle), sin(angle), -sin(angle), cos(angle), 0, 0);
        for (int i = 0; i < 3; i++) {
            int ori_pos_x = ori_pos.at<Vec<double, 1>>(0, 0)[0];
            int ori_pos_y = ori_pos.at<Vec<double, 1>>(1, 0)[0];
            if (ori_pos_x >= src.rows || ori_pos_y >= src.cols || ori_pos_x < 0 || ori_pos_y < 0) {
                pix[i] = 0;
            }
            else
                pix[i] = src.at<Vec3b>(ori_pos_x, ori_pos_y)[i];
        }
    });
    src = result;
}

Show

这里进行逆时针30度的变换
image_1e59m9ocv71e1mfgj70111kofc2c.png-45.2kBimage_1e59mamsb1asb1l0k3hnnvr1ou32p.png-28.9kB

Note

  • #include<math.h>中的sin()和cos()函数都是弧度制,可以根据需要进行转换

  1. 这里是一篇写伽马校正比较好的文章,我觉得可以作为背景知识补充。 

  2. 这句话原文是“処理後の画像で値が割り当てられない可能性がでてきてしまう。”直译大概是”处理后的图像可能没有被分配到值。“我也不知道该怎么翻译才好……你们看输出图像左下角黑色的那一块,就是这种没有被”分配“到的情况。 

为图像处理初学者设计的 100 个问题 问题三十一至问题四十

为图像处理初学者设计的 100 个问题 问题三十一至问题四十

标签(空格分隔): OpenCV

问题三十一至四十

问题三十一:仿射变换(Afine Transformations)——倾斜

  1. 使用仿射变换,输出(1)那样的$x$轴倾斜$30$度的图像($t_x=30$),这种变换被称为X-sharing。
  2. 使用仿射变换,输出(2)那样的y轴倾斜$30$度的图像($t_y=30$),这种变换被称为Y-sharing。
  3. 使用仿射变换,输出(3)那样的$x$轴、$y$轴都倾斜$30$度的图像($t_x = 30$,$t_y = 30$)。

原图像的大小为$h\ w$,使用下面各式进行仿射变换:

file

Answer

函数AfineTransformations()见上一章

void leanAt(Mat& src, double x_sharing, double y_sharing) {
    Mat result = Mat::zeros(src.rows, src.cols, CV_8UC3);
    result.forEach<Vec3b>([&](Vec3b& pix, const int * position) {
        int x = position[0];
        int y = position[1];
        int cur_pos[2][1] = { x,y };
        Mat m(Size(1, 2), CV_32SC1, cur_pos);
        double a_x = x_sharing / src.rows;
        double a_y = y_sharing / src.cols;
        Mat ori_pos = AfineTransformations(m, 1,a_x,a_y,1, 0, 0);
        //printMat(ori_pos);
        for (int i = 0; i < 3; i++) {
            int ori_pos_x = ori_pos.at<Vec<double, 1>>(0, 0)[0];
            int ori_pos_y = ori_pos.at<Vec<double, 1>>(1, 0)[0];
            if (ori_pos_x >= src.rows || ori_pos_y >= src.cols || ori_pos_x < 0 || ori_pos_y < 0) {
                pix[i] = 0;
            }
            else
                pix[i] = src.at<Vec3b>(ori_pos_x, ori_pos_y)[i];
        }
    });
    src = result;
}

Show

image_1e5bofanhfs54qe1c7t11jjq0r9.png-45.8kBimage_1e5bofkmk8nf10lp1ioq11921e3rm.png-40.9kB

Note

  • 利用OpenCV的API进行仿射变换

    • 第一步:获取仿射映射矩阵(两种)

      getAffineTransform()       //三点法
      Mat M1=getAffineTransform(const Point2f* src, const Point2f* dst)
        //参数const Point2f* src:原图的三个固定顶点

      //参数const Point2f dst:目标图像的三个固定顶点
      //返回值:Mat型变换矩阵,可直接用于warpAffine()函数
      //注意,顶点数组长度超过3个,则会自动以前3个为变换顶点;数组可用Point2f[]或Point2f
      表示

      
      ```cpp
      getRotationMatrix2D       //直接指定比例和角度
      Mat M2=getRotationMatrix2D (CvPoint2D32f  center,double angle,double scale)
          //参数CvPoint2D32f  center,表示源图像旋转中心
          //参数double angle,旋转角度,正值表示逆时针旋转
          //参数double scale,缩放系数
    • 第二步:进行仿射变换
      C++ void warpAffine(InputArray src, OutputArray dst, InputArray M, Size dsize, int flags=INTER_LINEAR, int borderMode=BORDER_CONSTANT, const Scalar& borderValue=Scalar())
      //参数InputArray src:输入变换前图像
      //参数OutputArray dst:输出变换后图像,需要初始化一个空矩阵用来保存结果,不用设定矩阵尺寸
      //参数InputArray M:变换矩阵,用另一个函数getAffineTransform()计算
      //参数Size dsize:设置输出图像大小
      //参数int flags = INTER_LINEAR:设置插值方式,默认方式为线性插值(另一种WARP_FILL_OUTLIERS)
      //参数int borderMode=BORDER_CONSTANT:边界像素模式,默认值BORDER_CONSTANT
      //参数const Scalar& borderValue=Scalar(),在恒定边界情况下取的值,默认值为Scalar(),即0
      第五个参数插值方式可选方式: 标识符 含义
      INTER_NEAREST 最邻近插值
      INTER_LINEAR 线性插值(默认)
      INTER_CUBIC 三次样条插值
      INTER_AREA 区域插值
      INTER_LANCZOS4 LancZoc4 插值
      INTER_LINEAR_EXACT 双线性插值
      INTER_MAX 插补码掩码
      WARP_FILL_OUTLIERS 填充所有目标图像像素。如果它们中的一些对应于源图像,它们被设置为0
      WARP_INVERSE_MAP 逆变换

问题三十二:傅立叶变换(Fourier Transform)

使用离散二维傅立叶变换(Discrete Fourier Transformation),将灰度化的imori.jpg表示为频谱图。然后用二维离散傅立叶逆变换将图像复原。

二维离散傅立叶变换是傅立叶变换在图像处理上的应用方法。通常傅立叶变换用于分离模拟信号或音频等连续一维信号的频率。但是,数字图像使用$[0,255]$范围内的离散值表示,并且图像使用$H\times W$的二维矩阵表示,所以在这里使用二维离散傅立叶变换。

二维离散傅立叶变换使用下式计算,其中$I$表示输入图像:
$$
G(k,l)=\frac{1}{H\ W}\ \sum\limits{y=0}^{H-1}\ \sum\limits{x=0}^{W-1}\ I(x,y)\ e^{-2\ \pi\ i\ (\frac{k\ x}{W}+\frac{l\ y}{H})}
$$
在这里让图像灰度化后,再进行离散二维傅立叶变换。

频谱图为了能表示复数$G$,所以图上所画长度为$G$的绝对值。这回的图像表示时,请将频谱图缩放至$[0,255]$范围。

二维离散傅立叶逆变换从频率分量$G$按照下式复原图像:
$$
I(x,y)=\frac{1}{H\ W}\ \sum\limits{l=0}^{H-1}\ \sum\limits{k=0}^{W-1}\ G(l,k)\ e^{2\ \pi\ i\ (\frac{k\ x}{W}+\frac{l\ y}{H})}
$$

上式中$\exp(j)$是个复数,实际编程的时候请务必使用下式中的绝对值形态1

如果只是简单地使用for语句的话,计算量达到$128^4$,十分耗时。如果善用NumPy的化,则可以减少计算量(答案中已经减少到$128^2$)。

Principe

本节Principe中的内容适合不了解傅里叶变换的读者阅读

 傅里叶变换可以将一副图片分解为正弦和余弦两个分量,换言之,它可以将一幅图像从空间域(spatial domain)转换为频域(frequency domain)。这种变换的思想是任何函数可以很精确的接近无穷个sin()函数和cos()函数的和。
我们来梳理一下概念:

  • 空间域
    一般情况下,空间域的图像是f(x, y) = 灰度级(0-255),形象一点就是一个二维矩阵,每一个坐标对应一个颜色值
  • 频率域
    对于一个正弦信号,如果它的幅度变化很快,我们称之为高频信号,如果变化非常慢,我们称之为低频信号。迁移到图像中,图像哪里的幅度变化非常大呢?边界点或者噪声。所以我们说边界点和噪声是图像中的高频分量(这里的高频是指变化非常快,不是出现的次数多),图像的主要部分集中在低频分量。

由于图像变换的结果原点在边缘部分,不容易显示,所以将原点移动到中心部分。
那么结果便是中间一个亮点朝着周围发散开来,越远离中心位置的能量越低(越暗)。

对于二位图像其傅里叶变换公式如下:
$$
F(k,l)= \sum\limits{y=0}^{N-1}\ \sum\limits{x=0}^{N-1}\ f(x,y)\ e^{-2\ \pi\ j\ (\frac{k\ x}{N}+\frac{l\ y}{N})}
$$
$$
e^{ix}=cos x+i sin x
$$
式中 f(i, j)是图像空间域的值,F(k ,l)是频域的值。傅里叶转换的结果是复数,这也显示了傅里叶变换是一幅实数图像和虚数图像叠加或者是幅度图像和相位图像叠加的结果。

Answer

toGray()函数在之第一章声明过,也可是使用OpenCV自带的API来转换灰度图
  • 傅里叶变换

    //// Discrete Fourier transformation
    //这个奇怪的vector就是用来存转换后的频域的
    std::vector<
    std::vector<
    std::complex<double>>> DFT(Mat& src) {
    int width = src.cols;
    int height = src.rows;
    std::complex<double> val;
    Mat grey = toGray(src);
    std::vector<std::vector<std::complex<double>>> FTMap;
    Mat result = Mat::zeros(width, height, CV_8UC1);
    
    //j,k用来表示频谱图中的坐标
    for (int j = 0; j < width; j++) {
        std::vector<std::complex<double>> MapLine;
        for (int k = 0; k < height; k++) {
            val.real(0);
            val.imag(0);
    
            //x,y表示时域中的坐标
            for (int x = 0; x < width; x++) {
                for (int y = 0; y < height; y++) {
                    uchar I = grey.at<Vec<uchar, 1>> (x, y)[0];
                    double theta = (double)-2 * 3.1415926*(((double)k*(double)y / (double)height) + ((double)j*(double)x / (double)width));
                    val += (double)I*std::complex<double>(cos(theta), sin(theta));
                }
            }
            val /= sqrt((double)width*(double)height);
            result.at<Vec<uchar, 1>>(j, k)[0] = (uchar)sqrt(val.imag()*val.imag() + val.real()*val.real());
            MapLine.push_back(val);
        }
        FTMap.push_back(MapLine);
    }
    //输出以下现在是什么样的图案
    imshow("FDT", result);
    return FTMap;
    }
  • 傅里叶逆变换

    Mat IDFT(std::vector<std::vector<std::complex<double>>> FTMap) {
    Mat result = Mat::zeros(128, 128, CV_8UC1);
    int width=128, height = 128;
    
    std::complex<double> val;
    
    for (int j = 0; j < width; j++) {
        std::vector<std::complex<double>> MapLine;
        for (int k = 0; k < height; k++) {
            val.real(0);
            val.imag(0);
    
            //x,y表示时域中的坐标
            for (int x = 0; x < width; x++) {
                for (int y = 0; y < height; y++) {
                    std::complex<double> G = FTMap[x][y];
                    double theta = (double)2 * 3.1415926*(((double)k*(double)y / (double)height) + ((double)j*(double)x / (double)width));
                    val += G*std::complex<double>(cos(theta), sin(theta));
                }
            }
            result.at<Vec<uchar, 1>>(j, k)[0] = (uchar)(std::abs(val) / sqrt(height * width));
        }
    }
    return result;
  • 在main函数中这样调用就可以了

    Mat img = imread("G:\\cv\\ImageProcessing100Wen-master\\Question_21_30\\imori.jpg");
    imshow("原图", img);
    imshow("effect",IDFT(DFT(img)));

SHow

  • 原图
    image_1e5fh45mu3ltgrr1qql1dc3lrl9.png-45.2kB

  • 转换后的频域图
    image_1e5fh4v6vtfq8pm1d91bideoim.png-23.4kB

  • 逆变换的复原图
    image_1e5fh5kvh18811qsk92hcu1rv113.png-24.7kB

Note

  • 可以看到这种方法来求傅里叶变换以及傅里叶逆变换的时间复杂度非常高,达到$O(n^4)$,所以常用的方法是FFT,ji快速傅里叶变换,在之后的题目中会介绍到

问题三十三:傅立叶变换——低通滤波

imori.jpg灰度化之后进行傅立叶变换并进行低通滤波,之后再用傅立叶逆变换复原吧!

通过离散傅立叶变换得到的频率在左上、右上、左下、右下等地方频率较低,在中心位置频率较高。2

!\[lpf\](assets/lpf.png)

在图像中,高频成分指的是颜色改变的地方(噪声或者轮廓等),低频成分指的是颜色不怎么改变的部分(比如落日的渐变)。在这里,使用去除高频成分,保留低频成分的低通滤波器吧!

在这里,假设从低频的中心到高频的距离为$r$,我们保留$0.5\ r$​的低频分量。

Answer

傅里叶变换和逆变换在前面已经实现过了,这一题直接调用

//低频滤波
Mat LowPassFilter(Mat& src, double pass_r) {
auto map = DFT(src);
int width = src.cols;
int height = src.rows;
for (int i = 0; i < width/2; i++) {
for (int j = 0; j < height/2; j++) {
if (sqrt(i * i + j * j) >= pass_r * height / 2) {
map[i][j] = 0;
map[width - i-1][j] = 0;
map[i][height - j-1] = 0;
map[width - i - 1][height - j - 1] = 0;
}
}
}
return IDFT(map);
}

Show

image_1e5fikpgu8jc1231a96dupo4g23.png-22kB

问题三十四:傅立叶变换——高通滤波

imori.jpg灰度化之后进行傅立叶变换并进行高通滤波,之后再用傅立叶逆变换复原吧!

在这里,我们使用可以去除低频部分,只保留高频部分的高通滤波器。假设从低频的中心到高频的距离为$r$,我们保留$0.2\ r$​的低频分量。

Answer

//高通滤波
Mat HighPassFilter(Mat& src, double pass_r) {
    auto map = DFT(src);
    int width = src.cols;
    int height = src.rows;
    for (int i = 0; i < width / 2; i++) {
        for (int j = 0; j < height / 2; j++) {
            if (sqrt(i * i + j * j) <= pass_r * height / 2) {
                map[i][j] = 0;
                map[width - i - 1][j] = 0;
                map[i][height - j - 1] = 0;
                map[width - i - 1][height - j - 1] = 0;
            }
        }
    }
    return IDFT(map);
}

Show

image_1e5fj0jir1g6u1r5q3r915hv1jae36.png-24.1kB

问题三十五:傅立叶变换——带通滤波

imori.jpg灰度化之后进行傅立叶变换并进行带通滤波,之后再用傅立叶逆变换复原吧!

在这里,我们使用可以保留介于低频成分和高频成分之间的分量的带通滤波器。在这里,我们使用可以去除低频部分,只保留高频部分的高通滤波器。假设从低频的中心到高频的距离为$r$,我们保留$0.1\ r$至$0.5\ r$的分量。

Answer

//带通滤波
Mat BandPassFilter(Mat& src, double pass_start, double pass_end) {
    auto map = DFT(src);
    int width = src.cols;
    int height = src.rows;
    for (int i = 0; i < width / 2; i++) {
        for (int j = 0; j < height / 2; j++) {
            if (sqrt(i * i + j * j) <= pass_start * height / 2|| sqrt(i * i + j * j)>= pass_end * height / 2) {
                map[i][j] = 0;
                map[width - i - 1][j] = 0;
                map[i][height - j - 1] = 0;
                map[width - i - 1][height - j - 1] = 0;
            }
        }
    }
    return IDFT(map);
}

Show

image_1e5fjd71kn791kbh1vgsm0k15k03s.png-23.5kB

问题三十六:JPEG 压缩——第一步:离散余弦变换(Discrete Cosine Transformation)

imori.jpg灰度化之后,先进行离散余弦变换,再进行离散余弦逆变换吧!

离散余弦变换(Discrete Cosine Transformation)是一种使用下面式子计算的频率变换:
file
离散余弦逆变换(Inverse Discrete Cosine Transformation)是离散余弦变换的逆变换,使用下式定义。

在这里,$K$是决定图像复原时分辨率高低的参数。$K=T$时,DCT的系数全被保留,因此IDCT时分辨率最大。$K=1$或$K=2$时,图像复原时的信息量(DCT系数)减少,分辨率降低。如果适当地设定$K$,可以减小文件大小。
file
在这里我们先将图像分割成$8\times8$的小块,在各个小块中使用离散余弦变换编码,使用离散余弦逆变换解码,这就是 JPEG的编码过程。现在我们也同样地,把图像分割成$8\times8$的小块,然后进行离散余弦变换和离散余弦逆变换。

这一整段我整体都在瞎**译,原文如下:

ここでは画像を8x8ずつの領域に分割して、各領域で以上のDCT, IDCTを繰り返すことで、JPEG符号に応用される。 今回も同様に8x8の領域に分割して、DCT, IDCTを行え。

——gzr

问题三十七:PSNR

离散余弦逆变换中如果不使用$8$作为系数,而是使用$4$作为系数的话,图像的画质会变差。来求输入图像和经过离散余弦逆变换之后的图像的峰值信噪比吧!再求出离散余弦逆变换的比特率吧!

峰值信噪比(Peak Signal to Noise Ratio)缩写为PSNR,用来表示信号最大可能功率和影响它的表示精度的破坏性噪声功率的比值,可以显示图像画质损失的程度。

峰值信噪比越大,表示画质损失越小。峰值信噪比通过下式定义。MAX表示图像点颜色的最大数值。如果取值范围是$[0,255]$的话,那么MAX的值就为255。MSE表示均方误差(Mean Squared Error),用来表示两个图像各个像素点之间差值平方和的平均数:
$$
\text{PSNR}=10\ \log{10}\ \frac{{v{max}}^2}{\text{MSE}}\
\text{MSE}=\frac{\sum\limits{y=0}^{H-1}\ \sum\limits{x=0}^{W-1}\ [I_1(x,y)-I_2(x,y)]^2}{H\ W}
$$
如果我们进行$8\times8$的离散余弦变换,离散余弦逆变换的系数为$KtimesK$的话,比特率按下式定义:
$$
\text{bit rate}=8\ \frac{K^2}{8^2}
$$


  1. 这里应该有个公式的,但是它不知道去哪儿了。 

  2. 下面图里文字的意思:高周波=高频;低周波=低频;入れ替え=替换。 

为图像处理初学者设计的 100 个问题 问题十一至问题二十

问题十一至问题二十

问题十一:均值滤波器

使用$3\times3$的均值滤波器来进行滤波吧!

均值滤波器使用网格内像素的平均值。

Answer

void AvgFilter(Mat src, Mat& dst,int pad_size) {
    if (!src.data) return;
    int height = src.rows;
    int width = src.cols;
    int channel = src.channels();
    dst = Mat::zeros(height+height, width+width, CV_8UC3);
    for (int x = 0; x < height - pad_size; x++) {
        for (int y = 0; y < width - pad_size; y++) {
            for (int c = 0; c < channel; c++) {
                Mat pad = src(Rect(x, y, pad_size, pad_size));
                int avg = sum(pad)[c] / (pad_size*pad_size);
                dst.at<Vec3b>( y,x)[c] = avg;
            }   
        }
    }
}

Show

image_1e4r1j7f6117o1akpldredm2l9.png-83.1kB

Note

  • 可以使用函数
    void cv::blur( InputArray src, OutputArray dst,Size ksize, Point anchor = Point(-1,-1),int borderType = BORDER_DEFAULT );j

问题十二:Motion Filter

使用$3\times3$的Motion Filter来进行滤波吧。

Motion Filter取对角线方向的像素的平均值,像下式这样定义:
$$
\left[
\begin{matrix}
\frac{1}{3}&0&0\
0&\frac{1}{3}&0\
0 & 0& \frac{1}{3}
\end{matrix}
\right]
$$

Answer

void MotionFilter(Mat src, Mat& dst, int pad_size) {
    if (!src.data) return;
    int height = src.rows;
    int width = src.cols;
    int channel = src.channels();
    dst = Mat::zeros(height + height, width + width, CV_8UC3);
    //initPad
    Mat pad = Mat::eye(pad_size, pad_size, CV_32F)/pad_size;
    pad.at<float>(0, 0);
    for (int x = 0; x < height-1; x++) {
        for (int y = 0; y < width-1; y++) {
            for (int c = 0; c < src.channels(); c++) {
                if (x - 1 >= 0 && y - 1 >= 0) {
                    for (int idx = 0; idx < pad_size; idx++)
                        dst.at<Vec3b>(x, y)[c] += src.at<Vec3b>(x - 1 + idx, y - 1 + idx)[c] * (pad.at<float>(idx, idx));
                }
                else {
                    dst.at<Vec3b>(x, y)[c] = src.at<Vec3b>(x, y)[c];
                }
            }
        }
    }
}

Note

  • OpenCV的forEach实现采用了并行的方法,所以可以在单独使用一个像素时使用,如果使用该像素时,需要依赖旁边的其他像素,那么可能造成访问冲突

Show

image_1e4r44gfjfpo1ckt1fmp1jme1arim.png-87.2kB

问题十三:MAX-MIN滤波器

使用MAX-MIN滤波器来进行滤波吧。

MAX-MIN滤波器使用网格内像素的最大值和最小值的差值对网格内像素重新赋值。通常用于边缘检测

边缘检测用于检测图像中的线。像这样提取图像中的信息的操作被称为特征提取

边缘检测通常在灰度图像上进行。

Principle

image_1e4r4t8s41lf710jm13cd113prqt13.png-19.2kB
我们知道,图像的细节属于低频信息,图像的边缘属于高频信息。我们使用一定大小的 Max-Min 滤波器作用于图像,当滤波器作用于图像细节时,输出结果往往趋向于0(黑色);而滤波器作用于图像边缘时,Max-Min 输出结果往往趋向于255(白色)。所以 最大-最小滤波器 能有效地用于检测图像的边缘和轮廓。

Answer

void Max_MinFilter(Mat src, Mat& dst, int pad_size) {
    if (!src.data) return;
    int height = src.rows;
    int width = src.cols;
    int channel = src.channels();
    dst = Mat::zeros(height + height, width + width,CV_8UC1);
    src=toGray(src);//第二问,将图像转换为灰度图像
    for (int x = 0; x < height - 1; x++) {
        for (int y = 0; y < width - 1; y++) {
            if (x - 1 >= 0 && y - 1 >= 0) {
                Mat pad = src(Rect( x - 1, y - 1, pad_size, pad_size)).clone().reshape(1, 1);
                Mat afterSort;
                sort(pad, afterSort, CV_SORT_EVERY_ROW + CV_SORT_ASCENDING);
                uchar max = afterSort.at<uchar>(0, pad_size*pad_size - 1);
                uchar min = afterSort.at<uchar>(0, 0);
                if (min == 0) {
                    std::cout << "?";
                }
                dst.at<Vec<uchar,1>>(y,x)[0] = max - min;
            }
            else {
                dst.at<Vec<uchar, 1>>(y, x)[0] = src.at<Vec<uchar, 1>>(y, x)[0];
            }

        }
    }
}

Show

  • pad_size为3
    image_1e4sd49prao1mrt1h41e1ul5j1g.png-67kB
  • pad_size为5
    image_1e4sdlbpc1cs512vf1b0k1lc4j3s2a.png-64.3kB

Note

  • 正如Principe所说,常用在边缘检测.虽然叫滤波器,但实际上并不能像之前的均值滤波,高斯滤波那样消除掉噪点,反而会是椒盐噪点扩大.如图:
    image_1e4sdbhc8qimfj4h1p1tt7r7j1t.png-72.5kB

问题十四:差分滤波器1(Differential Filter)

使用$3\times3$的差分滤波器来进行滤波吧。

差分滤波器对图像亮度急剧变化的边缘有提取效果,可以获得邻接像素的差值。

纵向:
$$
K=\left[
\begin{matrix}
0&-1&0\
0&1&0\
0&0&0
\end{matrix}
\right]
$$
横向:
$$
K=\left[
\begin{matrix}
0&0&0\
-1&1&0\
0&0&0
\end{matrix}
\right]
$$

Answer

void DifferentialFilter(Mat src, Mat& dst, int pad_size=3) {
    if (!src.data) return;
    int height = src.rows;
    int width = src.cols;
    int channel = src.channels();
    dst = Mat::zeros(height + height, width + width, CV_8UC1);
    src = toGray(src);

    //initKernel
    Mat pad = Mat::zeros(pad_size, pad_size, CV_32SC1);
    pad.at<Vec<int, 1>>(0, 1) = -1;
    pad.at<Vec<int, 1>>(1, 1)[0] = 1;
    printMat(pad);
    // filtering
    for (int x = 0; x < height; x++) {
        for (int y = 0; y < width; y++) {
            if (x - pad_size / 2 >= 0 && y - pad_size / 2 >= 0 && y + pad_size / 2 < width&& x + pad_size / 2 < height) {
                Mat rect = src(Rect(x - pad_size / 2, y - pad_size / 2, pad_size, pad_size));
                dst.at<Vec<uchar, 1>>( y, x)[0] =abs(rect.at<Vec<uchar, 1>>(1, 1)[0] * pad.at<Vec<int, 1>>(1, 1)[0] + rect.at<Vec<uchar, 1>>(0, 1)[0] * pad.at<Vec<int, 1>>(0, 1)[0]);

            }
            else {
                dst.at<Vec<uchar, 1>>(x, y)[0] = src.at<Vec<uchar, 1>>(x,y)[0];
            }
        }
    }
}

Note

Mat初始化的时候要注意每个像素的数据类型,否则会导致使用at函数时定位出现偏差.猜测内部实现原理是是通过offset*sizeof()来定位的.

矩阵数据类型:
    – CV_(S|U|F)C
S = 符号整型 U = 无符号整型 F = 浮点型

CV_8UC1 是指一个8位无符号整型单通道矩阵,
CV_32FC2是指一个32位浮点型双通道矩阵
CV_8UC1 CV_8SC1 CV_16U C1 CV_16SC1
CV_8UC2 CV_8SC2 CV_16UC2 CV_16SC2
CV_8UC3 CV_8SC3 CV_16UC3 CV_16SC3
CV_8UC4 CV_8SC4 CV_16UC4 CV_16SC4
CV_32SC1 CV_32FC1 CV_64FC1
CV_32SC2 CV_32FC2 CV_64FC2
CV_32SC3 CV_32FC3 CV_64FC3
CV_32SC4 CV_32FC4 CV_64FC4
其中,通道表示每个点能存放多少个数,类似于RGB彩色图中的每个像素点有三个值,即三通道的。
图片中的深度表示每个值由多少位来存储,是一个精度问题,一般图片是8bit(位)的,则深度是8.
1--bit_depth---比特数---代表8bite,16bites,32bites,64bites---举个例子吧--比如说,如
        如果你现在创建了一个存储--灰度图片的Mat对象,这个图像的大小为宽100,高100,那么,现在这张
        灰度图片中有10000个像素点,它每一个像素点在内存空间所占的空间大小是8bite,8位--所以它对
        应的就是CV_8
2--S|U|F--S--代表---signed int---有符号整形
               U--代表--unsigned int--无符号整形
               F--代表--float---------单精度浮点型
3--C----代表---一张图片的通道数,比如:
         1--灰度图片--grayImg---是--单通道图像
         2--RGB彩色图像---------是--3通道图像
         3--带Alph通道的RGB图像--是--4通道图像

所以在内部想要保存int类型整数是,要使用CV_32S*,32代表的是Bit数,一个int是4byte,32bit

Show

image_1e4sff5uq17o193m6jj1mbh18i12n.png-67.5kB

问题十五:Sobel滤波器

使用$3\times3$的Sobel滤波器来进行滤波吧。

Sobel滤波器可以提取特定方向(纵向或横向)的边缘,滤波器按下式定义:

纵向:
$$
K=\left[
\begin{matrix}
1&2&1\
0&0&0\
-1&-2&-1
\end{matrix}
\right]
$$
横向:
$$
K=\left[
\begin{matrix}
1&0&-1\
2&0&-2\
1&0&-1
\end{matrix}
\right]
$$

Answer

因为后面由多到题目用到给定模板进行滤波,所以编写几个函数

  • 根据kernel计算需要做filter的Mat中间值

    //_T 是图像的数据类型,8通道的即为Vec<uchar,1>
    //_Tk是Kernel的数据类型,因为kernel中需要出现负数,所以我采用了Vec<int,1>
    template<typename _T,typename _Tk>
    _T& filterMul(const Mat& src, const Mat& kernal) {
    _T result;
    int num = 0;
    if (!src.data) return result;
    int height = src.rows;
    int width = src.cols;
    if (height != kernal.rows || width != kernal.cols) return result;//两个形状要一样
    int channel = src.channels();
    for (int x = 0; x < height; x++) {
        for (int y = 0; y < width; y++) {
            _T a = src.at<_T>(x, y);
            _Tk b = kernal.at<_Tk>(x, y);
            num += a[0] * b[0];
        }
    }
    if (num < 0) num = 0;
    result[0] = num;
    
    return result;
    }
  • 将kernel与整个图像做卷积求计算中间像素值

    Mat filterConv(Mat& src, Mat& kernel) {
    Mat dst;
    if (!src.data) return dst;
    int height = src.rows;
    int width = src.cols;
    int pad_size = kernel.rows;
    dst = Mat::zeros(height*2, width*2, CV_8UC1);
    
    // filtering
    for (int x = 0; x < height; x++) {
        for (int y = 0; y < width; y++) {
            if (x - pad_size / 2 >= 0 && y - pad_size / 2 >= 0 && y + pad_size / 2 < width&& x + pad_size / 2 < height) {
                Mat rect = src(Rect(x - pad_size / 2, y - pad_size / 2, pad_size, pad_size));
                dst.at<Vec<uchar, 1>>(  y, x) = filterMul<Vec<uchar, 1>,Vec<int,1>>(rect, kernel);
    
            }
            else {
                dst.at<Vec<uchar, 1>>(x, y) = src.at<Vec<uchar, 1>>(x, y);
            }
        }
    }
    return dst;
    }
  • 初始化kernel

    template<typename _T>
    Mat initSingleChannelMat(_T* inputArray,int rows, int type) {
    Mat r = Mat::zeros(rows, rows, type);
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < rows; j++) {
            r.at<Vec<_T,1>>(i, j) = *(inputArray + j + rows * i);
        }
    }
    return r;
    }

所以最后的答案长这个样子

void SobelFilter(Mat src, Mat& dst, int pad_size = 3) {
    if (!src.data) return;
    src = toGray(src);
    //initKernel
    int s[] = {
            1,2,1,
            0,0,0,
            -1,-2,-1
    };
    Mat pad=initSingleChannelMat<int>(s,pad_size, CV_32SC1);

    // filtering
    dst = filterConv(src, pad);

}

之后的几个答案都会用到上述的几个工具函数

Show

image_1e4tjottfunr1cug12t21gs11d26i.png-66.2kB

Note

  • 在使用at()函数的时候,如果没有取到对应位置的值一定记得检查一下Mat的数据类型

问题十六:Prewitt滤波器

使用$3\times3$的Prewitt滤波器来进行滤波吧。

Prewitt滤波器是用于边缘检测的一种滤波器,使用下式定义:

纵向:
$$
K=\left[
\begin{matrix}
-1&-1&-1\
0&0&0\
1&1&1
\end{matrix}
\right]
$$
横向:
$$
K=\left[
\begin{matrix}
-1&0&-1\
-1&0&1\
-1&0&1
\end{matrix}
\right]
$$

Answer

有了上面几个工具函数,这样的问题只需要更换kernel即可

void PrewittFilter(Mat src, Mat& dst, int pad_size = 3) {
    if (!src.data) return;
    src = toGray(src);
    //initKernel
    int s[] = {
            -1,-1,-1,
            0,0,0,
            1,1,1
    };
    Mat pad = initSingleChannelMat<int>(s, pad_size, CV_32SC1);

    // filtering
    dst = filterConv(src, pad);

}

Show

image_1e4tk0kl0192v1qnb19kv16b116s5v.png-65.7kB

问题十七:Laplacian滤波器

使用Laplacian滤波器来进行滤波吧。

Laplacian滤波器是对图像亮度进行二次微分从而检测边缘的滤波器。由于数字图像是离散的,$x$方向和$y$方向的一次微分分别按照以下式子计算:
$$
I_x(x,y)=\frac{I(x+1,y)-I(x,y)}{(x+1)-x}=I(x+1,y)-I(x,y)\
Iy(x,y) =\frac{I(x, y+1) - I(x,y)}{(y+1)-y}= I(x, y+1) - I(x,y)
$$
因此二次微分按照以下式子计算:
$$
\begin{align*}
&I
{xx}(x,y) \
=& \frac{I_x(x,y) - I_x(x-1,y)}{(x+1)-x} \
=& I_x(x,y) - Ix(x-1,y)\
=&[I(x+1, y) - I(x,y)] - [I(x, y) - I(x-1,y)]\
=& I(x+1,y) - 2\ I(x,y) + I(x-1,y)
\end{align*}
$$
同理:
$$
I
{yy}(x,y)=I(x,y+1)-2\ I(x,y)+I(x,y-1)
$$
特此,Laplacian 表达式如下:
$$
\begin{align}
&\nabla^2\ I(x,y)\
=&I{xx}(x,y)+I{yy}(x,y)\
=&I(x-1,y) + I(x,y-1) - 4
I(x,y) + I(x+1,y) + I(x,y+1)
\end{align*}
$$
如果把这个式子表示为卷积核是下面这样的:
$$
K=
\left[
\begin{matrix}
0&1&0\
1&-4&1\
0&1&0
\end{matrix}
\right]
$$

Answer

void LaplacianFilter(Mat src, Mat& dst, int pad_size = 3) {
    if (!src.data) return;
    src = toGray(src);
    //initKernel
    int s[] = {
            0,1,0,
            1,-4,1,
            0,1,0
    };
    Mat pad = initSingleChannelMat<int>(s, pad_size, CV_32SC1);

    // filtering
    dst = filterConv(src, pad);

}

Show

image_1e4tk7sff3uh52d1irg3tltdt1l.png-69.7kB
image_1e4tkpoq21gao95717ni1i3v32c2b.png-478.6kB

问题十八:Emboss滤波器

使用Emboss滤波器来进行滤波吧。

Emboss滤波器可以使物体轮廓更加清晰,按照以下式子定义:
$$
K=
\left[
\begin{matrix}
-2&-1&0\
-1&1&1\
0&1&2
\end{matrix}
\right]
$$

Answer

void EmbossFilter(Mat src, Mat& dst, int pad_size = 3) {
    if (!src.data) return;
    src = toGray(src);
    //initKernel
    int s[] = {
            -2,-1,0,
            -1,1,1,
            0,1,2
    };
    Mat pad = initSingleChannelMat<int>(s, pad_size, CV_32SC1);

    // filtering
    dst = filterConv(src, pad);

}

Show

image_1e4tkumd93mulik51f17fh3s938.png-506.4kB

问题十九:LoG滤波器

使用LoG 滤波器,来对imori_noise.jpg检测边缘吧!

LoG即高斯-拉普拉斯(Laplacian of Gaussian)的缩写,使用高斯滤波器使图像平滑化之后再使用拉普拉斯滤波器使图像的轮廓更加清晰。

为了防止拉普拉斯滤波器计算二次微分会使得图像噪声更加明显,所以我们首先使用高斯滤波器来抑制噪声。

LoG 滤波器使用以下式子定义:
$$
\text{LoG}(x,y)=\frac{x^2 + y^2 - s^2}{2 \ \pi \ s^6} \ e^{-\frac{x^2+y^2}{2\ s^2}}
$$

Answer

void LoGFilter(Mat src, Mat& dst, int pad_size = 3) {
    if (!src.data) return;
    GaussianFilter(src, dst, 3, 10);
    Mat newSrc = Mat(dst);
    LaplacianFilter(newSrc, dst);
}

问题二十:直方图

使用Matplotlib来绘制imori_dark.jpg的直方图吧!

直方图显示了不同数值的像素出现的次数。在Matplotlib中有hist()函数提供绘制直方图的接口。

Answer

C++并没有可以直接绘制直方图的工具,所以利用Mat自己写了一个

/*
绘制直方图
*/
template<typename _T>
void paintHistogram(std::vector<_T> t,int height=480,int colWidth=5) {
    Mat m = Mat::zeros(height, t.size() *colWidth, CV_8UC3);
    _T maxItem = *std::max_element(t.begin(), t.end());
    Vec3b color = Vec3b(255, 255, 255);
    for (int item = 0; item < t.size(); item++) {
        int cur_height = ((float)t[item] / (float)maxItem)*height;
        for (int x = item * colWidth; x < (item + 1)*colWidth-3; x++) {
            for (int y = height-1; y >height-cur_height; y--) {
                m.at<Vec3b>( y,x)[1] = color[1];
                m.at<Vec3b>(y, x)[0] = color[0];
                m.at<Vec3b>(y, x)[2] = color[2];
            }
        }
    }
    imshow("Histogram", m);
}

有了上述工具绘制直方图就很简单

void paintGrayHistogram(Mat src) {
    Mat gray = toGray(src);
    std::vector<uchar> v;
    for (int i = 0; i < 256; i++) {
        v.push_back(0);
    }
    gray.forEach< Vec<uchar, 1>>([&](Vec<uchar, 1> pix, const int* position) {
        v[pix[0]] += 1;
    });
    paintHistogram(v);
}

Show

image_1e50actoh1c0s11ckl3p11lu1s70m.png-74.2kB

Note

  • 如果在运行过程中,出现图像展示时发生"堆已损坏","无法访问0xFFFFFF"之类的,很有可能是在给Mat赋值是越界了,仔细检查边界条件,看看是不是在位于width,和height的位置赋了值,x,y要小于等于width-1,height-1

  1. 原文是“微分フィルタ”,对应的英文应该是“Differential Filter”。考虑到这里处理离散的变量,应该是“差分滤波器”?问题十七:Laplacian滤波器中的“二次微分”没有什么太大歧义就没作处理。 

图像处理初学者设计的 100 个问题 原题+解析+笔记

问题一至问题十

问题一:通道交换

读取图像,然后将 '$ \text{RGB} $' 通道替换成$\text{BGR}$通道。

Answer

#include<opencv2\core\core.hpp>
#include<opencv2\highgui\highgui.hpp>
#include<opencv2\opencv.hpp>

using namespace cv;

Mat RGBToBGR(Mat mat) {
    int height = mat.rows;
    int width = mat.cols;
    Mat result = Mat::zeros(height, width, CV_8UC3);
    for (int y = 0; y < height; y++) {
        for (int x = 0; x < width; x++) {
            result.at<Vec3b>(y, x)[0] = mat.at<Vec3b>(y, x)[2];
            result.at<Vec3b>(y, x)[1] = mat.at<Vec3b>(y, x)[1];
            result.at<Vec3b>(y, x)[2] = mat.at<Vec3b>(y, x)[0];
        }
    }
    return result;

}
int main()
{
    Mat img = imread("C:\\Users\\PC\\Pictures\\1.jpg");
    imshow("原图", img);
    imshow("question_1", RGBToBGR(img));
    waitKey(0);
    destroyAllWindows();
}

Show

image_1e4j4465l1rco913mj8kl1180c9.png-561.7kB

Note

  • Vecused as element of multi-channel images

问题二:灰度化(Grayscale)

将图像灰度化吧!

灰度是一种图像亮度的表示方法,通过下式计算:

$$Y = 0.2126\ R + 0.7152\ G + 0.0722\ B$$

Answer

#include<opencv2\core\core.hpp>
#include<opencv2\highgui\highgui.hpp>
#include<opencv2\opencv.hpp>

Mat toGray(Mat mat) {
    Mat result = Mat::zeros(mat.rows, mat.cols, CV_8UC1);
    for (int y = 0; y < mat.rows; y++) {
        for (int x = 0; x < mat.cols; x++) {
            result.at<Vec<uchar, 1>>(y, x)[0] = mat.at<Vec3b>(y, x)[0] * 0.2126 + mat.at<Vec3b>(y, x)[1] * 0.7152 + mat.at<Vec3b>(y, x)[2] * 0.0722;
        }
    }
    return result;
}

int main()
{
    Mat img = imread("C:\\Users\\PC\\Pictures\\1.jpg");
    imshow("原图", img);
    imshow("question_2", toGray(img));
    waitKey(0);
    destroyAllWindows();
}

Show

image_1e4j4j5mjjv41qai2d4f59ia5m.png-466.4kB

Note

可以直接调用封装好的函数

    Mat img = imread("C:\\Users\\PC\\Pictures\\1.jpg");
    Mat imgGray;
    cvtColor(img, imgGray, CV_RGB2GRAY);

问题三:二值化(Thresholding)

把图像进行二值化吧。

二值化是将图像使用黑和白两种颜色表示的方法。

我们将灰度的阈值设置为$128$来进行二值化,即:

$$
y=
\begin{cases}
0& (\text{if}\quad y < 128) \
255& (\text{else})
\end{cases}
$$

#include<opencv2\core\core.hpp>
#include<opencv2\highgui\highgui.hpp>
#include<opencv2\opencv.hpp>

using namespace cv;

Mat Thresholding(Mat mat) {
    Mat result = Mat::zeros(mat.rows, mat.cols, CV_8UC1);
    Mat grayImg = toGray(mat);//接上一问灰度化
    for (int y = 0; y < mat.rows; y++) {
        for (int x = 0; x < mat.cols; x++) {
            result.at<Vec<uchar, 1>>(y, x)[0] = mat.at<Vec<uchar, 1>>(y, x)[0] > 128 ? 255 : 0;
        }
    }
    return result;
}

Show

image_1e4j5899m1co91ej51pl61u881dqi13.png-317.7kB

Note

可以使用 double cv::threshold( InputArray src, OutputArray dst,
double thresh, double maxval, int type );
函数
@param src Source 8-bit single-channel image.
@param dst Destination image of the same size and the same type as src.
@param maxValue Non-zero value assigned to the pixels for which the condition is satisfied
@param adaptiveMethod Adaptive thresholding algorithm to use, see #AdaptiveThresholdTypes.
The #BORDER_REPLICATE | #BORDER_ISOLATED is used to process boundaries.
@param thresholdType Thresholding type that must be either #THRESH_BINARY or #THRESH_BINARY_INV,
see #ThresholdTypes.
@param blockSize Size of a pixel neighborhood that is used to calculate a threshold value for the
pixel: 3, 5, 7, and so on.
@param C Constant subtracted from the mean or weighted mean (see the details below). Normally, it
is positive but may be zero or negative as well.

问题四:大津二值化算法(Otsu's Method)

使用大津算法来二值化图像吧。

大津算法,也被称作最大类间方差法,是一种可以自动确定二值化中阈值的算法。

类内方差类间方差的比值计算得来:

  • 小于阈值$t$的类记作$0$,大于阈值$t$的类记作$1$;
  • $w_0$和$w_1$是被阈值$t$分开的两个类中的像素数占总像素数的比率(满足$w_0+w_1=1$);
  • ${S_0}^2$, ${S_1}^2$是这两个类中像素值的方差;
  • $M_0$,$M_1$是这两个类的像素值的平均值;

即:

  • 类内方差:${S_w}^2=w_0\ {S_0}^2+w_1\ {S_1}^2$
  • 类间方差:${S_b}^2 = w_0 \ (M_0 - M_t)^2 + w_1\ (M_1 - M_t)^2 = w_0\ w_1\ (M_0 - M_1) ^2$
  • 图像所有像素的方差:${S_t}^2 = {S_w}^2 + {S_b}^2 = \text{常数}$

根据以上的式子,我们用以下的式子计算分离度$X$:1

$$ X = \frac{{S_b}^2}{{S_w}^2} = \frac{{S_b}^2}{{S_t}^2 - {S_b}^2} $$

也就是说:

\arg\max\limits_{t}\ X=\arg\max\limits_{t}\ {S_b}^2

换言之,如果使${S_b}^2={w_0}\ {w_1}\ (M_0 - M_1)^2$最大,就可以得到最好的二值化阈值$t$。

Answer

Mat thresholding_OTSU(Mat mat) {
    int height = mat.rows;
    int width = mat.cols;
    Mat result = Mat::zeros(height, width, CV_8UC1);
    Mat grayImg = toGray(mat);
    double S_temp = 0.0;
    int degree_temp = 0;
    for (int grayDegree = 0; grayDegree < 255; grayDegree++) {
        double downNum = 0.0;//小于阈值的像素个数
        double downSum = 0.0;//小于阈值的像素灰度值总和
        double upNum = 0.0;//大于阈值的像素个数
        double upSum = 0.0;//大于阈值的像素灰度值总和
        for (int y = 0; y < mat.rows; y++) {
            for (int x = 0; x < mat.cols; x++) {
                uchar pix = grayImg.at<Vec<uchar, 1>>(y, x)[0];
                if ( pix< grayDegree) {
                    downNum++;
                    downSum += pix;
                }
                else{
                    upNum++;
                    upSum += pix;
                }
            }
        }
        double w_0 = downNum / (downNum + upNum);
        double w_1 = upNum / (downNum + upNum);
        double M_0,M_1 = 0.0;
        if(downNum>0)  M_0 = downSum / downNum; else  M_0 = 0.0;
        if (upNum > 0) M_1 = upSum / upNum; else M_1 = 0.0;
        double S = w_0 * w_1*(M_0 - M_1)*(M_0 - M_1);
        if (S > S_temp) {
            S_temp = S;
            degree_temp = grayDegree;
        }
    }
    return Thresholding(grayImg, degree_temp);
}
}

Show

image_1e4jag20v8rplbgnv25kt1fum1t.png-317.6kB

Note

问题五:$\text{HSV}$变换

将使用$\text{HSV}$表示色彩的图像的色相反转吧!

$\text{HSV}$即使用色相(Hue)、饱和度(Saturation)、明度(Value)来表示色彩的一种方式。

  • 色相:将颜色使用$0^{\circ}$到$360^{\circ}$表示,就是平常所说的颜色名称,如红色、蓝色。色相与数值按下表对应:

    绿 青色 蓝色 品红
    $0^{\circ}$ $60^{\circ}$ $120^{\circ}$ $180^{\circ}$ $240^{\circ}$ $300^{\circ}$ $360^{\circ}$
  • 饱和度:是指色彩的纯度,饱和度越低则颜色越黯淡($0\leq S < 1$);

  • 明度:即颜色的明暗程度。数值越高越接近白色,数值越低越接近黑色($0\leq V < 1$);

从$\text{RGB}$色彩表示转换到$\text{HSV}$色彩表示通过以下方式计算:

$\text{RGB}$的取值范围为$[0, 1]$,令:
$$
\text{Max}=\max(R,G,B)\
\text{Min}=\min(R,G,B)
$$
色相:
$$
H=\begin{cases}
0&(\text{if}\ \text{Min}=\text{Max})\
60\ \frac{G-R}{\text{Max}-\text{Min}}+60&(\text{if}\ \text{Min}=B)\
60\ \frac{B-G}{\text{Max}-\text{Min}}+180&(\text{if}\ \text{Min}=R)\
60\ \frac{R-B}{\text{Max}-\text{Min}}+300&(\text{if}\ \text{Min}=G)
\end{cases}
$$
饱和度:
$$
S=\text{Max}-\text{Min}
$$
明度:
$$
V=\text{Max}
$$
从$\text{HSV}$色彩表示转换到$\text{RGB}$色彩表示通过以下方式计算:
image_1e4javjc71kft1mk0opi1jg6132c2a.png-53.7kB
请将色相反转(色相值加$180$),然后再用$\text{RGB}$色彩空间表示图片。

Answer

void RGBToHSV(Mat src, Mat& dst) {
    int height = src.rows;
    int width = src.cols;

    dst = Mat::zeros(src.rows, src.cols, CV_32FC3);
    for (int x = 0; x < height; x++) {
        for (int y = 0; y < width; y++) {
            float r = (float)src.at<Vec3b>(x, y)[0]/255;
            float g = (float)src.at<Vec3b>(x, y)[1]/255;
            float b = (float)src.at<Vec3b>(x, y)[2]/255;
            float max = std::max(r, std::max(g, b)); 
            float min = std::min(r, std::min(g, b)); 
            float H = 0.0;
            if (max == min) H = 0.0;
            else if (min == b) H = (60.0*(g - r) / (max - min)) + 60.0;
            else if (min == r) H = (60.0*(b - g) / (max - min)) + 180.0;
            else if (min == g) H = (60.0*(r - b) / (max - min)) + 300.0;

            float S = (float)max - (float)min;

            float V = (float)max;

            dst.at<Vec3f>(x, y)[0] = H;
            dst.at<Vec3f>(x, y)[1] = S;
            dst.at<Vec3f>(x, y)[2] = V;
        }
    }
}

Show

image_1e4l63s72gnfo5i14aj115i1m3m2n.png-524.2kB

Note

可以使用自带的函数cv::cvtColor(src,dst,CV_RGB2HSV)处理

问题六:减色处理2

我们将图像的值由$256^3$压缩至$4^3$,即将$\text{RGB}$的值只取${32, 96, 160, 224}$。这被称作色彩量化。色彩的值按照下面的方式定义:
$$
\text{val}=
\begin{cases}
32& (0 \leq \text{var} < 64)\
96& (64\leq \text{var}<128)\
160&(128\leq \text{var}<192)\
224&(192\leq \text{var}<256)
\end{cases}
$$

Answer

void subColor(Mat src, Mat& dst) {
    int height = src.rows;
    int width = src.cols;
    dst = Mat::zeros(height, width, CV_MAKETYPE(8, 3));
    for (int x = 0; x < height; x++) {
        for (int y = 0; y < width; y++) {
            uchar r = src.at<Vec3b>(x, y)[0];
            uchar g = src.at<Vec3b>(x, y)[1];
            uchar b = src.at<Vec3b>(x, y)[2];
            dst.at <Vec<uchar, 3 >>(x, y)[0] = (r > 64) ? (r > 128 ? (r > 192 ? 224 : 160) : 96) : 32;
            dst.at <Vec<uchar, 3 >>(x, y)[1] = (g > 64) ? (g > 128 ? (g > 192 ? 224 : 160) : 96) : 32;
            dst.at <Vec<uchar, 3 >>(x, y)[2] = (b > 64) ? (b > 128 ? (b > 192 ? 224 : 160) : 96) : 32;
        }
    }
}

Show

image_1e4m12v3qg4u1v2n13sm1uuv18dv34.png-378.3kB

Note

预定义类型的结构如下所示:
CV_(S|U|F)C
1–bit_depth—比特数—代表8bite,16bites,32bites,64bites—举个例子吧–比如说,如
如果你现在创建了一个存储–灰度图片的Mat对象,这个图像的大小为宽100,高100,那么,现在这张
灰度图片中有10000个像素点,它每一个像素点在内存空间所占的空间大小是8bite,8位–所以它对
应的就是CV_8
2–S|U|F–S--代表—signed int—有符号整形
U–代表–unsigned int–无符号整形
F–代表–float---------单精度浮点型
3–C----代表—一张图片的通道数,比如:
1–灰度图片–grayImg—是–单通道图像
2–RGB彩色图像---------是–3通道图像
3–带Alph通道的RGB图像–是--4通道图像

问题七:平均池化(Average Pooling)

将图片按照固定大小网格分割,网格内的像素值取网格内所有像素的平均值。

我们将这种把图片使用均等大小网格分割,并求网格内代表值的操作称为池化(Pooling)

池化操作是卷积神经网络(Convolutional Neural Network)中重要的图像处理方式。平均池化按照下式定义:
$$
v=\frac{1}{|R|}\ \sum\limits_{i=1}^R\ v_i
$$
请使用$8\times8$的网格做平均池化。

Answer

void AveragePooling(Mat src,Mat& dst) {
    int height = src.rows;
    int width = src.cols;
    int size = 8;
    dst = Mat::zeros(height, width, CV_8UC3);
    struct AssignColor
    {
    Scalar src;
    public:AssignColor(Scalar& src_) { 
        src = src_; 
    }

        void operator()(Vec3b& pixel, const int* position) const{
            pixel[0] = src[0];
            pixel[1] = src[1];
            pixel[2] = src[2];
        }
    };
    for (int x = 0; x < height- size; x+= size) {
        for (int y = 0; y < width- size; y+= size) {
            Mat conv = src(Rect( y, x, size, size)).clone();
            Scalar avg = mean(conv);
            Mat dstConv = dst(Rect(y, x, size, size));
            dstConv.forEach< Vec3b>(AssignColor(avg));
        }
    }
 }

Show

image_1e4n5et531fkj17j416r51lb71o23h.png-317.1kB

Note

  • 对于提取Mat中的矩阵可以使用mat(Rect(x,y,height,width)),mat为Mat对象.该提取后的对象做修改会影响原始图像,新建副本要使用clone()

  • Mat对象中有forEach()方法,可以新建结构体,构造void operator()函数;函数定义如下:

    template<typename _Tp, typename Functor> inline
    void Mat::forEach(const Functor& operation) {
    this->forEach_impl<_Tp>(operation);
    }

    _Tp是每个像素的数据类型,本题中使用的是Vec3b
    结构体构建可以参考如下

    public:AssignColor(Scalar& src_) { 
        src = src_; 
    }
        void operator()(Vec3b& pixel, const int* position) const{
            pixel[0] = src[0];
            pixel[1] = src[1];
            pixel[2] = src[2];
        }
    };

    因为我要给每个循环传递参数所以给该结构体一个构造函数

  • 注意上述结构体中void operator()函数后面加的const关键词,非静态成员函数后面加const(加到非成员函数或静态成员后面会产生编译错误),表示成员函数隐含传入的this指针为const指针,决定了在该成员函数中,任意修改它所在的类的成员的操作都是不允许的(因为隐含了对this指针的const引用);唯一的例外是对于mutable修饰的成员。加了const的成员函数可以被非const对象和const对象调用,但不加const的成员函数只能被非const对象调用。

问题八:最大池化(Max Pooling)

网格内的值不取平均值,而是取网格内的最大值进行池化操作。

Answer

void MaxPooling(Mat src, Mat& dst) {
    int height = src.rows;
    int width = src.cols;
    int size = 8;
    dst = Mat::zeros(height, width, CV_8UC3);
    struct AssignColor
    {
        Vec3b src;
    public:AssignColor(Vec3b& src_) {
        src = src_;
    }

           void operator()(Vec3b& pixel, const int* position) const {
               pixel[0] = src[0];
               pixel[1] = src[1];
               pixel[2] = src[2];
           }
    };
    for (int x = 0; x < height - size; x += size) {
        for (int y = 0; y < width - size; y += size) {
            Mat channal[3];
            Mat conv = src(Rect(y, x, size, size)).clone();
            split(conv, channal);
            Point* maxIndex=new Point();
            minMaxLoc(channal[0], NULL, NULL, NULL, maxIndex);
            Mat dstConv = dst(Rect(y, x, size, size));
            dstConv.forEach< Vec3b>(AssignColor(conv.at<Vec3b>(maxIndex->x,maxIndex->y)));
        }
    }
}

Show

image_1e4n7hqbigub1fj7kql11ql1gf047.png-314.4kB

Note

这一题偷懒了,按照第一个通道的最大值位置为每个通道的最大值位置.

问题九:高斯滤波(Gaussian Filter)

使用高斯滤波器($3\times3$大小,标准差$\sigma=1.3$)来对imori_noise.jpg进行降噪处理吧!

高斯滤波器是一种可以使图像平滑的滤波器,用于去除噪声。可用于去除噪声的滤波器还有中值滤波器(参见问题十),平滑滤波器(参见问题十一)、LoG滤波器(参见问题十九)。

高斯滤波器将中心像素周围的像素按照高斯分布加权平均进行平滑化。这样的(二维)权值通常被称为卷积核(kernel)或者滤波器(filter)

但是,由于图像的长宽可能不是滤波器大小的整数倍,因此我们需要在图像的边缘补$0$。这种方法称作Zero Padding。并且权值$g$(卷积核)要进行归一化操作($\sum\ g = 1$)。

按下面的高斯分布公式计算权值:

$$
g(x,y,\sigma)=\frac{1}{2\ \pi\ \sigma^2}\ e^{-\frac{x^2+y^2}{2\ \sigma^2}}
$$

标准差$\sigma=1.3$的$8-$近邻高斯滤波器如下:

$$
K=\frac{1}{16}\ \left[
\begin{matrix}
1 & 2 & 1 \
2 & 4 & 2 \
1 & 2 & 1
\end{matrix}
\right]
$$

Answer

void GaussianFilter(Mat src, Mat& dst, int kernel_size,float sigma) {
    int hetght = src.rows;
    int width = src.cols;
    dst = Mat::zeros(hetght, width, CV_8UC3);
    Mat kernel = Mat::zeros(kernel_size, kernel_size, CV_32F);
    const float pi = 3.1415926;

    //initKernel
    float kernelSum = 0.0;
    kernel.forEach<float>([&](float& item, const int* position) {
        int x = position[0]-(int)kernel_size/2;
        int y = position[1]- (int)kernel_size / 2;
        item = (1 / (2 * pi*sigma*sigma))*exp(-(x*x+y*y)/(2*sigma*sigma));
        kernelSum += item;
    });
    kernel.forEach<float>([&](float& item, const int* position) {
        item /= kernelSum;
    });
    std::cout << "\n" << sum << "\n";
    // filtering
    int channel = src.channels();
    for (int x = 0; x < hetght; x++) {
        for (int y = 0; y < width; y++) {
            for (int c = 0; c < channel; c++) {
                int v = 0;
                int pad = kernel_size / 2;
                for (int dx = -pad; dx < pad + 1; dx++) {
                    for (int dy = -pad; dy < pad + 1; dy++) {
                        if (((x + dx) >= 0) && ((y + dy) >= 0)) {
                            v += (double)src.at<cv::Vec3b>( x + dx, y + dy)[c] * kernel.at <float>(dx + pad, dy + pad);
                        }
                    }
                }
                dst.at<Vec3b>(x, y)[c] = v;
            }
        }
    }
}

Show

image_1e4pncc4gjon12ta1mpr1ghj1pl74k.png-85.6kB
$\sigma=1.3$时的滤波效果
image_1e4pnfsi9pvib4tdgfo5n1rfe51.png-84.9kB
$\sigma=10$时的滤波效果
image_1e4pnhdj81qmogtt2ei1dmu14445e.png-54.6kB
$\sigma=0.1$时的滤波效果

Note

  • 在C++中的lambda表达式中,要使用外面的任意变量,可以用&{}的方式,&代表捕获所有局部变量。

  • 可以调用函数GaussianBlur(src, GaussianBlurImg, Size(5, 5), 1, 1);

问题十:中值滤波(Median Filter)

使用中值滤波器($3\times3$大小)来对imori_noise.jpg进行降噪处理吧!

中值滤波器是一种可以使图像平滑的滤波器。这种滤波器用滤波器范围内(在这里是$3\times3$)像素点的中值进行滤波,请在这里也采用Zero Padding。

Answer

void MedianFilter(Mat src, Mat& dst, int pad_size) {
    if (!src.data) return;
    int height = src.rows;
    int width = src.cols;
    int channel = src.channels();
    dst = Mat::zeros(height + height, width + width, CV_8UC3);
    for (int x = 1; x < height - pad_size; x++) {
        for (int y = 1; y < width - pad_size; y++) {
            Mat pad = src(Rect(x-1, y-1, pad_size, pad_size));
            Mat chs[3];
            split(pad, chs);
            for (int c = 0; c < channel; c++) {
                Mat flut;
                chs[c].reshape(1, 1).copyTo(flut);
                sort(flut, flut,CV_SORT_EVERY_ROW + CV_SORT_ASCENDING);
                int mid = flut.at<Vec<uchar,1>>(0, (pad_size*pad_size) / 2)[0];
                dst.at<Vec3b>(y, x)[c] = mid;
            }
        }
    }
}

Show

image_1e4pq6no39cl61p26l41f5od6h.png-82.6kB

Note

  • 可以调用函数medianBlur(src, MedianBlurImg, 5);进行中值滤波
  • 要对Mat进行降维可以使用mat.reshape()


  1. 这里原repo配图里的公式好像打错了。 

  2. 这里没有找到"減色処理"准确的中文翻译,所以直译了。