她的性格像一片乡土之地,初次相识的时候觉得挺广大,又很荒凉,不过很快你便能在里面发现连绵起伏的巍峨群山中的,被果树包围着的怡人的小村子和惬意地流过茂密牧场的欢快的小河。尽管这些美景让你惊叹,甚至能够让你打消疑虑,但站在这被大风侵袭的褐土高坡上,却不足以让你有家的感觉。
作者: 避堇
终于终于,从事了小时候梦想的职业,坐进了宽敞明亮的办公室,当上了都市白领,拿到了远超均值的薪资,却并没有曾经期望中的那么开心,甚至有些难过。
不过过的似乎还挺滋润的。。。每天睡到自然醒!!!
CPython源码解读——前言
开一个大大的坑,不知道能不能写完。加油!
文章结构
本文主要分为以下几个部分
- 一、基础知识介绍
- 二、Python对象实现
- 三、Python虚拟机
- 四、Pyhton标准库的实现
其中基础知识介绍中主要介绍与C相关的前置基础知识,以及Python的基础知识。该部分均为零散的知识点,阅读后面部分时有看不懂的地方可以在该部分查找。
Python对象实现:在 Python 的世界一切皆对象,不论是整数,还是字符串,甚至连类型、函数等都是一种对象。这一节从源码的角度,介绍这些对象的实现原理。
Python虚拟机 在介绍完对象的实现原理后,介绍Python如何解释你的Py文件中的代码。该部分将包括虚拟机的编译处理,内存管理,线程管理等部分。
Pyhton标准库的实现 包含了所有用 C 语言编写的模块,比如 math、hashlib 等。它们都是那些对速度要求非常严格的模块。也包括Lib 目录下则是存放一些对速度没有太严格要求的模块,比如 os。
当然,整个文章还包括这篇前言。前言介绍了文章结构和一些开始前的准备工作。读者可以跟着我一起做好准备工作,这样我们可以一起读源码了。
对了,本文基于的系统环境时Window 10。部分内容与Linux会稍有区别
本文阅读需要以下前置知识:
- 【了解】编译原理,计算机组织原理
- 【熟悉】C语言,操作系统
- 【掌握】Python语法
下载CPython源码
https://github.com/python/cpython
源码里有什么
cpython/
│
├── Doc ← 源代码文档说明
├── Grammar ← 计算机可读的语言定义
├── Include ← C 语言头文件(头文件中一般放一些重复使用的代码)
├── Lib ← Python 写的标准库文件
├── Mac ← Mac 支持的文件
├── Misc ← 杂项
├── Modules ← C 写的标准库文件
├── Objects ← 核心类型和对象模块
├── Parser ← Python 解析器源码
├── PC ← Windows 编译支持的文件
├── PCbuild ← 老版本的 Windows 系统 编译支持的文件
├── Programs ← Python 可执行文件和其他二进制文件的源代码
├── Python ← CPython 解析器源码
└── Tools ← 用于构建或扩展 Python 的独立工具
其中:
- Include 目录:包含了 Python 提供的所有头文件,如果用户需要自己用 C 或 C++来编写自定义模块扩展 Python,那么就需要用到这里提供的头文件。
- Lib 目录:包含了 Python 自带的所有标准库,且都是用 Python 语言编写的。
- Modules 目录:包含了所有用 C 语言编写的模块,比如 math、hashlib 等。它们都是那些对速度要求非常严格的模块。而相比而言,Lib 目录下则是存放一些对速度没有太严格要求的模块,比如 os。
- Parser 目录:包含了 Python 解释器中的 Scanner 和 Parser 部分,即对 Python 源代码进行词法分析和语法分析的部分。除此以外,此目录还包含了一些有用的工具,这些工具能够根据 Python 语言的语法自动生成 Python 语言的词法和语法分析器,与 YACC 非常类似。
- Objects 目录:包含了所有 Python 的内建对象,包括整数、list、dict 等。同时,该目录还包括了 Python 在运行时需要的所有的内部使用对象的实现。
- Python 目录:包含了 Python 解释器中的 Compiler 和执行引擎部分,是 Python 运行的核心所在。
- PCbuild 目录:包含了 Visual Studio 2003 的工程文件,研究 Python 源代码就从这里开始(本文采用 Visual Studio 2017 对 Python 进行编译)。
- Programs 目录:包含了 Python 二进制可执行文件的源码。
太阳坟
一
今天是星与父亲退役前的最后一晚.
四年前星才刚刚十六岁,和四十岁的父亲一起参了军,日夜坚守在边塞。戍边之城建立在一道深渊面前,对面便是敌人的堡垒。据说两座城池如此僵持已有百年。这四年里,星时时刻刻都在想象着深渊的对面有什么,那里的的人长什么样子,说什么语言,是不是也有一对父子在城里久久不能回家。
两座城没有真正意义上打过仗,更多地只是威慑对方,这点星见识过。白天,对面架出炮筒朝这边开炮,只可惜这个深渊又宽又深,连炮弹都无法飞跃。开上一轮炮后就各回各家了。晚上,对面趁着天地安静下来,开始齐声歌唱他们的国歌与军歌,开始文化与信仰输出。当然,星这边也没少唱。有时候,星甚至怀疑自己参加的是军事体验夏令营而不是战场。四年里,他没有参加过一场真正意义上的战争,终于迎来了和平退役的日子。
这四年里,星逐渐熟悉了这样的生活,他明白这座城存在的意义,也明白了明明所有的攻击都无济于事却为什么还要输出。唯一令他困惑的是,这四年里,每天都需要派“重兵”来守夜,有时候甚至值守人员数量甚至超过了白天。星没少执过“夜班”,他总会和父亲一起在夜里值守。星在这漫长夤夜里,数过天上的星星,跟父亲一起认星座,隔岸看对面的灯火。深夜湿漉漉的寒风吹得衣服裹在身上。有时候,星幻想自己变成一条巨龙,一只只在东方微亮时出现的巨龙,借助氤氲的水汽,在广袤的平原上飞驰,划破天空,划破大地,划破这漫长的黑夜,在太阳出来的那一刹那,消散成一股灵气,唤醒新的一天。
今天是星与父亲退役前的最后一晚。只要东方泛起白光,新的一天到来时,星的军人生涯就宣告结束了。星想起了自己的母亲,姐姐还有弟弟,已经整整四年没有见他们了。今天的夜晚似乎特别漫长,星望着东方苦苦等着,激动又忐忑。希望今晚赶紧过去,又希望留住这段时光的最后一刻。
可是日出依旧没有到来。
人们发现了情况的不对。沉睡中的士兵发现今天居然睡了个自然醒,外面还是没有亮。守夜的士兵肚子饿得咕咕叫,可是今晚的执勤还没有结束。报时博士看着刻漏开始陷入了慌张,无法理解到底发生了什么。城中骚动了起来,人们惊恐又无助,不知道要做些什么。有人说是深渊里不干净的东西爬了出来挡住了太阳,很快就要攻击人类的城池。有人收拾行李马匹,准备逃离此地。还有人跪在地上,双目紧闭双手合十,潜心祈祷着什么。守夜军也开始骚动了起来。
"Rage!Rage against the dying of the light!"
军营中逐渐亮起了火把。
二
桑一名是太阳神教的传教士,游历各国传教布道,宣扬太阳神教。桑认为自己的教意不言而喻。每天太阳必然升起也必然落下,阳光普照大地滋养万物,没有太阳就没有人间的一切。所有的生命都要向太阳朝拜,谨以诚意。桑把自己生命献给了这份自己喜爱的事业,从入教第一天起,发誓虔心修行,弘扬教义,永不叛教,愿意用一生的时间去守护神圣的太阳教。
太阳神教已经积累了几千万的教徒,他们无一不崇拜着太阳神,视太阳为神圣的图腾。每到周一清早,大家一同举行朝觐仪式,祭祀伟大的太阳。每到晚上,教堂内的教徒齐聚一堂,为落下的太阳的太阳讼唱经文,祈愿明日新日冉冉升起。
“我们深爱着太阳,太阳用他的光辉爱着万物。”,桑主持着今晚的讼经,“我愿意将自己全部的生命奉献给伟大的太阳,我愿意将自己毕生的事业奉献给传播神教,即使有一天太阳降下天谴,惩罚卑微弱小的人类,我也无怨无悔,依旧潜心朝觐,永远热爱......祝愿明早的太阳依旧照常升起。”
三
星从没觉得星河有什么可滚烫的。知道今天才见识了,不但滚烫而且煎熬。星意识到,从此时此刻开始,自己的历史,人类的历史就要变得不一样了。
军营中各单位紧急集合,灯火通明,说是为了防止敌人趁乱偷袭。其实士兵们都知道,哪有什么敌人,不过是自己心中的恐惧与不安罢了。主事的将领也慌了神,让手下清点粮草燃料,集中起来严格监管,禁止滥用,并赶紧派了人马赴京报告情况。此去京城不过五百里里,快马加鞭日夜兼程,不出两日便可回还。但是将领很快意识到,就算京城得到消息也无济于事,如果是全天下的太阳都落了,京城也必是一片恐慌,那有功夫搭理自己的信使呢。
同样意识到这点的,还有星的父亲。他站在人群中,双目紧闭,面色凝重。星问父亲怎么了,父亲摇摇头说,“既然眼睛看不见,那就用别的感官去看,从此恐怕没有生存,只有幸存。我多想回到家乡,再看一眼你娘,只恐怕回不去了”。
星看看双目紧闭的父亲,欲言又止,再看看天空,天空居然没有半点乌云,星河璀璨。月亮落山了,没有月亮的晚上星星居显得尤为明亮,甚至有些耀眼。星意识到,从此时此刻开始,自己的历史,人类的历史就要变得不一样了。
当太阳没有出现的时候,其他怪异的事情就出现了。
为图像处理初学者设计的 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
以及两个直方图比较
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
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
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
问题二十五:最邻近插值( Nearest-neighbor Interpolation )
使用最邻近插值将图像放大$1.5$倍吧!
最近邻插值在图像放大时补充的像素取最临近的像素的值。由于方法简单,所以处理速度很快,但是放大图像画质劣化明显。
使用下面的公式放大图像吧!$I'$为放大后图像,$I$为放大前图像,$a$为放大率,方括号是四舍五入取整操作:
$$
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
Note
放大后有明显马赛克痕迹,可以尝试加入之前习题中的滤波器
以下为放大6倍并均值滤波效果
问题二十六:双线性插值( Bilinear Interpolation )
使用双线性插值将图像放大$1.5$倍吧!
双线性插值考察$4$邻域的像素点,并根据距离设置权值。虽然计算量增大使得处理时间变长,但是可以有效抑制画质劣化。
-
放大后图像的座标$(x',y')$除以放大率$a$,可以得到对应原图像的座标$(\lfloor \frac{x'}{a}\rfloor , \lfloor \frac{y'}{a}\rfloor)$。
-
求原图像的座标$(\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)$:
-
分别求这4个点与$(\frac{x'}{a}, \frac{y'}{a})$的距离,根据距离设置权重:$w = \frac{d}{\sum\ d}$
-
根据下式求得放大后图像$(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倍
问题二十七:双三次插值( 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
Note
- 在进行矩阵乘法运算时要注意两个Mat 的数据类型,比如说CV_64FC1与CV_32SC1之间无法直接相乘,会产生如下错误:
可以使用src.convertTo(src, CV_64FC1);
函数来进行主句类型之间的转换.其中src是一个Mat对象的实例.
问题二十九:仿射变换( Afine Transformations )——放大缩小
- 使用仿射变换,将图片在$x$方向上放大$1.3$倍,在$y$方向上缩小至原来的$\frac{4}{5}$。
- 在上面的条件下,同时在$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倍
问题三十:仿射变换( Afine Transformations )——旋转
- 使用仿射变换,逆时针旋转$30$度。
- 使用仿射变换,逆时针旋转$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度的变换
Note
#include<math.h>
中的sin()和cos()函数都是弧度制,可以根据需要进行转换
为图像处理初学者设计的 100 个问题 问题三十一至问题四十
为图像处理初学者设计的 100 个问题 问题三十一至问题四十
标签(空格分隔): OpenCV
问题三十一至四十
问题三十一:仿射变换(Afine Transformations)——倾斜
- 使用仿射变换,输出(1)那样的$x$轴倾斜$30$度的图像($t_x=30$),这种变换被称为X-sharing。
- 使用仿射变换,输出(2)那样的y轴倾斜$30$度的图像($t_y=30$),这种变换被称为Y-sharing。
- 使用仿射变换,输出(3)那样的$x$轴、$y$轴都倾斜$30$度的图像($t_x = 30$,$t_y = 30$)。
原图像的大小为$h\ w$,使用下面各式进行仿射变换:
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
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
-
原图
-
转换后的频域图
-
逆变换的复原图
Note
- 可以看到这种方法来求傅里叶变换以及傅里叶逆变换的时间复杂度非常高,达到$O(n^4)$,所以常用的方法是FFT,ji快速傅里叶变换,在之后的题目中会介绍到
问题三十三:傅立叶变换——低通滤波
将imori.jpg
灰度化之后进行傅立叶变换并进行低通滤波,之后再用傅立叶逆变换复原吧!
通过离散傅立叶变换得到的频率在左上、右上、左下、右下等地方频率较低,在中心位置频率较高。2
在图像中,高频成分指的是颜色改变的地方(噪声或者轮廓等),低频成分指的是颜色不怎么改变的部分(比如落日的渐变)。在这里,使用去除高频成分,保留低频成分的低通滤波器吧!
在这里,假设从低频的中心到高频的距离为$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
问题三十四:傅立叶变换——高通滤波
将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
问题三十五:傅立叶变换——带通滤波
将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
问题三十六:JPEG 压缩——第一步:离散余弦变换(Discrete Cosine Transformation)
imori.jpg
灰度化之后,先进行离散余弦变换,再进行离散余弦逆变换吧!
离散余弦变换(Discrete Cosine Transformation)是一种使用下面式子计算的频率变换:
离散余弦逆变换(Inverse Discrete Cosine Transformation)是离散余弦变换的逆变换,使用下式定义。
在这里,$K$是决定图像复原时分辨率高低的参数。$K=T$时,DCT的系数全被保留,因此IDCT时分辨率最大。$K=1$或$K=2$时,图像复原时的信息量(DCT系数)减少,分辨率降低。如果适当地设定$K$,可以减小文件大小。
在这里我们先将图像分割成$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}
$$
为图像处理初学者设计的 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
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
问题十三:MAX-MIN滤波器
使用MAX-MIN滤波器来进行滤波吧。
MAX-MIN滤波器使用网格内像素的最大值和最小值的差值对网格内像素重新赋值。通常用于边缘检测。
边缘检测用于检测图像中的线。像这样提取图像中的信息的操作被称为特征提取。
边缘检测通常在灰度图像上进行。
Principle
我们知道,图像的细节属于低频信息,图像的边缘属于高频信息。我们使用一定大小的 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
- pad_size为5
Note
- 正如Principe所说,常用在边缘检测.虽然叫滤波器,但实际上并不能像之前的均值滤波,高斯滤波那样消除掉噪点,反而会是椒盐噪点扩大.如图:
问题十四:差分滤波器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
问题十五: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
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
问题十七: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
问题十八: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
问题十九: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
Note
- 如果在运行过程中,出现图像展示时发生"堆已损坏","无法访问0xFFFFFF"之类的,很有可能是在给Mat赋值是越界了,仔细检查边界条件,看看是不是在位于width,和height的位置赋了值,x,y要小于等于width-1,height-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
Note
- 类
Vec
used 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
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
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
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}$色彩表示通过以下方式计算:
请将色相反转(色相值加$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
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
Note
预定义类型的结构如下所示:
CV_
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
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
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
$\sigma=1.3$时的滤波效果
$\sigma=10$时的滤波效果
$\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
Note
- 可以调用函数medianBlur(src, MedianBlurImg, 5);进行中值滤波
- 要对Mat进行降维可以使用mat.reshape()
错题集-计算机组成原理
计算机系统概述
- 【2009】冯诺依曼计算机中指令和数据均以二进制形式存放在存储器中,CPU区分它们的依据是()
A 指令操作码的译码结果
B 指令和数据的寻址方式
C 指令周期的不同
D 指令和数据所在的存储单元
C
- 在运算器中,不包含()
A 状态寄存器
B 数据总线
C ALU
D 地址寄存器
D
- 下列不属于系统软件的是():
A 数据库系统
B 操作系统
C 编译程序
D 以上3个都是
A 数据库系统包括数据库,数据库管理系统,用户程序等,只有数据库管理系统是系统软件
- 【2018】冯诺依曼结构计算机数据采用二进制编码表示,其主要原因是()
I 二进制运算规则简单
II 制造两个稳态的物理期间比较容易
III 便于用逻辑门电路实现算数运算
A I,II
B I,III
C II,III
D I,II,III
D,运算规则简单居然也是原因.....
- 32位机指的是该计算机所用的CPU()
A 具有32位寄存器
B 能同时处理32位的二进制
C 具有32个寄存器
D 能处理32个字符
B 机器字长进行一次整数运算(即定点整数运算)所能处理的二进制数据的位数,通常与CPU寄存器位数,加法器有关.
- 在CPU的寄存器中,()对于用户是完全透明的.
A 程序计数器
B 指令寄存器
C 状态寄存器
D 通用寄存器
B
IR用来存放当前指令,其内容来自主存的MAR
取指令:PC->MAR->M->MDR->IR
分析指令:OP(IR)->CU 将IR中的操作码送到控制淡雅UN
执行指令:Ad(IR)->MAR->M->MDR->ACC
- 下列书法正确的是
I 在微信计算机的广泛应用中,会计电算化属于科学计算方面的应用
II 决定计算机计算精度的重要技术是计算机的字长
III 计算机"运行速度"指标的含义是每秒能执行多少条操作系统的命令
IV 利用大规模集成电路技术把计算机的运算器部件和控制部件做在一块集成电路芯片上,这样的一块芯片叫做单片机
A I III
B II.IV
C II
D I,III,IV
C
决定计算机计算精度的重要技术就是计算机的字长
单片机还集成了RAM,ROM,I/O等,题中说的那是CPU
设计模式
工厂模式
工厂模式(Factory Pattern)是 Java 中最常用的设计模式之一。这种类型的设计模式属于创建型模式,它提供了一种创建对象的最佳方式。
在工厂模式中,我们在创建对象时不会对客户端暴露创建逻辑,并且是通过使用一个共同的接口来指向新创建的对象。
主要解决:主要解决接口选择的问题。
优点: 1、一个调用者想创建一个对象,只要知道其名称就可以了。 2、扩展性高,如果想增加一个产品,只要扩展一个工厂类就可以。 3、屏蔽产品的具体实现,调用者只关心产品的接口。
缺点:每次增加一个产品时,都需要增加一个具体类和对象实现工厂,使得系统中类的个数成倍增加,在一定程度上增加了系统的复杂度,同时也增加了系统具体类的依赖。这并不是什么好事。
使用场景: 1、日志记录器:记录可能记录到本地硬盘、系统事件、远程服务器等,用户可以选择记录日志到什么地方。 2、数据库访问,当用户不知道最后系统采用哪一类数据库,以及数据库可能有变化时。 3、设计一个连接服务器的框架,需要三个协议,"POP3"、"IMAP"、"HTTP",可以把这三个作为产品类,共同实现一个接口。
注意事项:作为一种创建类模式,在任何需要生成复杂对象的地方,都可以使用工厂方法模式。有一点需要注意的地方就是复杂对象适合使用工厂模式,而简单对象,特别是只需要通过 new 就可以完成创建的对象,无需使用工厂模式。如果使用工厂模式,就需要引入一个工厂类,会增加系统的复杂度。
单例模式
单例模式(Singleton Pattern)是 Java 中最简单的设计模式之一。这种类型的设计模式属于创建型模式,它提供了一种创建对象的最佳方式。
这种模式涉及到一个单一的类,该类负责创建自己的对象,同时确保只有单个对象被创建。这个类提供了一种访问其唯一的对象的方式,可以直接访问,不需要实例化该类的对象。
意图:保证一个类仅有一个实例,并提供一个访问它的全局访问点。
主要解决:一个全局使用的类频繁地创建与销毁。
何时使用:当您想控制实例数目,节省系统资源的时候。
如何解决:判断系统是否已经有这个单例,如果有则返回,如果没有则创建。
关键代码:构造函数是私有的。
应用实例:
- 1、一个班级只有一个班主任。
- 2、Windows 是多进程多线程的,在操作一个文件的时候,就不可避免地出现多个进程或线程同时操作一个文件的现象,所以所有文件的处理必须通过唯一的实例来进行。
- 3、一些设备管理器常常设计为单例模式,比如一个电脑有两台打印机,在输出的时候就要处理不能两台打印机打印同一个文件。
优点:
- 1、在内存里只有一个实例,减少了内存的开销,尤其是频繁的创建和销毁实例(比如管理学院首页页面缓存)。
- 2、避免对资源的多重占用(比如写文件操作)。
缺点:没有接口,不能继承,与单一职责原则冲突,一个类应该只关心内部逻辑,而不关心外面怎么样来实例化。
使用场景:
- 1、要求生产唯一序列号。
- 2、WEB 中的计数器,不用每次刷新都在数据库里加一次,用单例先缓存起来。
- 3、创建的一个对象需要消耗的资源过多,比如 I/O 与数据库的连接等。
注意事项:getInstance() 方法中需要使用同步锁 synchronized (Singleton.class) 防止多线程同时进入造成 instance 被多次实例化。
建造者模式
主要解决:主要解决在软件系统中,有时候面临着"一个复杂对象"的创建工作,其通常由各个部分的子对象用一定的算法构成;由于需求的变化,这个复杂对象的各个部分经常面临着剧烈的变化,但是将它们组合在一起的算法却相对稳定。
适配器模式
适配器模式(Adapter Pattern)是作为两个不兼容的接口之间的桥梁。这种类型的设计模式属于结构型模式,它结合了两个独立接口的功能。
这种模式涉及到一个单一的类,该类负责加入独立的或不兼容的接口功能。举个真实的例子,读卡器是作为内存卡和笔记本之间的适配器。您将内存卡插入读卡器,再将读卡器插入笔记本,这样就可以通过笔记本来读取内存卡。
意图:将一个类的接口转换成客户希望的另外一个接口。适配器模式使得原本由于接口不兼容而不能一起工作的那些类可以一起工作。

观察者模式
当对象间存在一对多关系时,则使用观察者模式(Observer Pattern)。比如,当一个对象被修改时,则会自动通知依赖它的对象。观察者模式属于行为型模式。
意图:定义对象间的一种一对多的依赖关系,当一个对象的状态发生改变时,所有依赖于它的对象都得到通知并被自动更新。
主要解决:一个对象状态改变给其他对象通知的问题,而且要考虑到易用和低耦合,保证高度的协作。
何时使用:一个对象(目标对象)的状态发生改变,所有的依赖对象(观察者对象)都将得到通知,进行广播通知。
策略模式
意图:定义一系列的算法,把它们一个个封装起来, 并且使它们可相互替换。
主要解决:在有多种算法相似的情况下,使用 if...else 所带来的复杂和难以维护。
何时使用:一个系统有许多许多类,而区分它们的只是他们直接的行为。
如何解决:将这些算法封装成一个一个的类,任意地替换。
关键代码:实现同一个接口。