下沉市场模式解析

定义：五环以外的市场，以及三线城市和三线城市以下的市场。
特点：人口基数大（国内超过十亿）
例子：1. 拼多多，创立三年上市，市值和京东相仿；2. 趣头条，创立两年零三月上市；3.快手，估值超过200亿刀 4.水滴筹，估值超过10亿 5.云集（赚钱模式是拉下线，下线购物和销售自己能拿到佣金的百分之10）

第一个挑战：现在了解市场不再是通过人的经验，而是通过人工智能和数据分析
第二个挑战：自己不一定要是自己产品的用户（比如张一鸣不一定用今日头条）
底层逻辑本质不在于了解用户，而在于：1.数据分析 2. 建立商业模型
产品各行业举例：电商拼多多；搜索百度；游戏有专门针对的产品；社交微信；OYO连锁酒店
产品创始人举例：做下沉市场最牛的人都不是在三四线城市出生，而是在北京上海；都不是草根出生，而是高学历高智商。

  • 拼多多黄峥天才少年，
  • 快手宿华清华博士退学，曾任职google百度
  • 趣头条CEO谭思亮清华本科，曾任雅虎盛大高管。

大局上看共性：
下沉市场人口超10亿，大部分的家庭月收入不足3000，
但是

  1. 可支配收入不低，因为开支不多
  2. 广告价值不低
  3. 闲暇时间多达六小时
  4. 没有品牌概念买东西图实惠，很多人根本不知道自己的电视品牌，手机持有率0.5
  5. 手机拥有率和App下载率低，扫码下载不现实，多为朋友帮忙装
  6. 对品质要求不高，对价格敏感
  7. 通勤时间短，15min内（所以原本针对通勤时间听的知识付费类产品难以攻入）

面对一线城市做，需要帮他们节约时间，提高效率；
面对下沉市场做，需要帮他们消磨时间，过得快乐。
总之：

  1. 有红利，获客成本低，获得用户方式：拼团、收徒、人拉人给补贴。（但补贴只是降低获客成本的方法）
  2. 建立合适的商业模型变现。比如使用周期中看多少条广告（带来利润）+留存+拉人头+裂变；算得，如果用户在半年内如果能给公司赚6毛，那我就能拿出3毛来分给用户。

c++实现最大堆和最小堆

堆是具有以下特性的完全二叉树，每个结点的值都大于或等于其左右孩子结点的值，叫做最大堆；每个结点的值都小于或等于其左右孩子结点的值，叫做最小堆。

一、用vector和push_heap、pop_heap实现堆

（1）建堆

1
vector<int> nums = {9, 6, 2, 4, 7, 0, 1, 8, 3, 5};

1、如果使用nums构建最大堆：

1
2
3
make_heap(nums.begin(), nums.end());
//或
make_heap(nums.begin(), nums.end(), less<int>());

输出nums的结果为

1
2
//最大堆是按照层序遍历的顺序存入vector的
9 8 2 6 7 0 1 4 3 5

2、如果使用nums构建最小堆：

1
make_heap(nums.begin(), nums.end(), greater<int>());

输出nums的结果为

1
0 3 1 4 5 2 9 8 6 7

（2）调整堆

当使用上述的make_heap()建完堆后，如果vector使用push_back()插入数据或pop_back()删除数据后，会破坏最大堆/最小堆的性质，所以需要调整堆，常用push_heap()和pop_heap()两个方法

1、push_heap()用法是，vector先push_back()，后push_heap()：

1
2
nums.push_back(10);
push_heap(nums.begin(), nums.end(), less<int>());

输出nums的结果：

1
2
3
4
5
6
//原vector 
9 8 2 6 7 0 1 4 3 5 
//push_back()后
9 8 2 6 7 0 1 4 3 5 10
//push_heap()后
10 9 2 6 8 0 1 4 3 5 7

2、pop_heap()用法是，先pop_heap()，vector后pop_back()：

1
2
pop_heap(nums.begin(), nums.end(), less<int>());
nums.pop_back();

输出nums的结果：

1
2
3
4
5
6
//原vector
9 8 2 6 7 0 1 4 3 5 
//pop_heap()后
8 7 2 6 5 0 1 4 3 9 
//pop_back()后
8 7 2 6 5 0 1 4 3

为什么pop_heap()的用法要反过来呢？
要从我们的目的来考虑，使用pop_heap()的绝大部分目的是要把堆顶元素pop出堆中，因为它最大或最小。如果先用vector的pop_back()，它删除的不是堆顶元素（nums[0]），而是vector的最后一个元素。可见这不是我们想要的结果：对于最大堆，最后一个元素既不是最大，也不一定是最小；对于最小堆，最后一个元素既不是最小，也不一定是最大。pop出来没有意义。
观察pop_heap()对堆做了什么？
pop_heap()把堆顶元素放到了最后一位，然后对它前面的数字重建了堆。这样一来只要再使用pop_back()把最后一位元素删除，就得到了新的堆。

二、用priority_queue实现堆

priority_queue
对于这个模板类priority_queue，它是STL所提供的一个非常有效的容器。
作为队列的一个延伸，优先队列包含在头文件 中。

（1）简述

优先队列时一种比较重要的数据结构，它是有二项队列编写而成的，可以以O(log n) 的效率查找一个队列中的最大值或者最小值，其中是最大值还是最小值是根据创建的优先队列的性质来决定的。

（2）模板参数

优先队列有三个参数，其声明形式为：

1
priority_queue< type, container, function >

这三个参数，后面两个可以省略，第一个不可以。
其中：

  1. type：数据类型；
  2. container：实现优先队列的底层容器；
  3. function：元素之间的比较方式；

对于container，要求必须是数组形式实现的容器，例如vector、deque，而不能使list。
在STL中，默认情况下（不加后面两个参数）是以vector为容器，以 operator< 为比较方式，所以在只使用第一个参数时，优先队列默认是一个最大堆，每次输出的堆顶元素是此时堆中的最大元素。

（3）成员函数

假设type类型为int，则：

  1. bool empty() const
    返回值为true，说明队列为空；
  2. int size() const
    返回优先队列中元素的数量；
  3. void pop()
    删除队列顶部的元素，也即根节点
  4. int top()
    返回队列中的顶部元素，但不删除该元素；
  5. void push(int arg)
    将元素arg插入到队列之中；

（4）大顶堆和小顶堆

  • 大顶堆
1
2
3
4
5
//构造一个空的优先队列（此优先队列默认为大顶堆）
priority_queue<int> big_heap;   
//另一种构建大顶堆的方法
priority_queue<int,vector<int>,less<int> > big_heap2;
  • 小顶堆
1
2
//构造一个空的优先队列,此优先队列是一个小顶堆
priority_queue<int,vector<int>,greater<int> > small_heap;

需要注意的是，如果使用less<int>和less<int>，需要头文件：
#include <functional>

转载自文章
c++使用vector建立最大堆和最小堆

C++中priority_queue理解与使用

C++string的substr操作需要注意的地方

s.substr(pos,n)

返回一个string，包含s中从pos开始的n个字符的拷贝。pos的默认值为0，n的默认值为s.size()-pos，即拷贝从pos开始的所有字符。

需要注意的是pos不能越界，pos+n可以越界。如果开始位置pos超过了string的大小，则substr函数会抛出一个out_of_range异常。如果开始位置pos加上计数值n大于string的大小，则substr会调整计数值，只拷贝到string的末尾。

压缩感知SL0图像重建算法的opencv实现

最近研究生的项目中涉及到压缩感知的SL0图像重建算法，需要对其进行C++实现。

按照压缩感知的理论框架，可以将它分为三个部分：图像的稀疏表示，即寻找一个正交基使原始图像尽可能的稀疏，图像的稀疏表示越充分，就越有利于图像的重建；测量矩阵的设计，为了重构稀疏信号，编码采样测量矩阵必须满足约束等距性RIP(Restricted isometry property)条件， 即将编码采样视为对原有信号的一种映射变换；图像重建，即通过测量信号重构原始信号。

重建算法相当于信号测量过程的逆过程，设一个长度为N的原始信号经过测量矩阵得到一个长度为M(M<N)的测量值信号，高分辨率成像重建的过程是通过测量向量y重构原始信号x的过程，由于方程个数远少于求解未知量个数，必须解一个欠定方程，很难直接求解。但由于原始信号是稀疏的或者能够稀疏表示，所以通常会加一个稀疏表达作为正则项，而稀疏表达通常用求解最小L0范数的方式加以表示，如下式所示：

$$
\hat{x}=\arg \min |x|_{0} \quad \text { s.t. } \quad \Phi x=y
$$

但求解最小L0范数通常很困难，是一个NP难问题，比如长度为N的稀疏信号中有K个非零值，则稀疏信号有$\mathrm{C}_{\mathrm{N}}^{\mathrm{K}}$种排列可能，要找到最接近于原始信号的最优排列，计算复杂度非常高。因此为了高分辨率成像重建的高效进行，常常会规避零范数问题，比如会选择以最小L1范数法为代表的次最优解算法以及以正交匹配追踪算法为代表的的贪婪算法进行求解，此外还常用迭代阈值法和最小全变分法。也可以用一种渐近逼近L0范数的高分辨率成像重建算法，通过渐近化解非凸优化求解的困难，通过L0约束确保信号的稀疏度。

信号重构算法是指由长度为M测量向量的y重构长度为N(M<N)的稀疏信号s的过程。

$$
\hat{\mathrm{s}}=\operatorname{argmin}|\mathrm{x}|_{0} \quad \text { s.t. } \mathrm{y}=\Theta \mathrm{s}
$$

项目中涉及到的SL0算法是由最小L2范数向最小L0范数的逐渐逼近，上式中Θ是传感矩阵，$|\mathrm{x}|_{0}$表示的是L0范数，研究首先定义一个逼近函数模型：

$$
\mathrm{f}{\sigma}(s)=1-\exp \left(-s^{2} / 2 \sigma^{2}\right) \quad \mathrm{f}{\sigma}\left(s{i}\right)=1-\exp \left(-s{i}^{2} / 2 \sigma^{2}\right)
$$

从而有$\mathrm{F}{\sigma}(\mathrm{s})=\sum{\mathrm{i}=1}^{\mathrm{n}} \mathrm{f}{\sigma}\left(\mathrm{s}{\mathrm{i}}\right)$，其中，n为稀疏信号向量s的长度，si为稀疏信号向量s对应的第i个元素，σ为逼近L0范数的调制参数，$\mathrm{F}{\sigma}(\mathrm{s})$近似表示s非零较大项的数目。逼近函数$\mathrm{f}{\sigma}\left(\mathrm{s}{\mathrm{i}}\right)$和$\mathrm{s}{\mathrm{i}}$、σ的关系如下图所示：

当σ较大时，$\mathrm{F}{\sigma}(\mathrm{s})$可以近似表示为L2范数。根据渐进思想，当σ取值逐渐减小时，$\mathrm{F}{\sigma}(\mathrm{s})$逐渐逼近L0范数，向量s的L0范数准则的优化问题可近似表示为$|\mathrm{s}|{0}=\mathrm{F}{\sigma}(\mathrm{s})$。另外从上图可知，由于逼近的L0范数的函数曲线是光滑可导的，因此被称为逼近光滑L0范数算法（简称为SL0）。

由此，稀疏表示$|\mathrm{x}|{0}$极小问题，转化为连续函数$\mathrm{F}{\sigma}(\mathrm{s})$的极小。从而重建式转化为一种全新的稀疏表示模型：

$$
\hat{\mathbf{s}}=\arg \min \mathrm{F}_{\sigma}(\mathrm{s}) \text { s.t } \Theta \mathrm{s}=\mathrm{y}
$$

此模型通常只针对于普通的稀疏表示，对于图像重构这一类的逆问题，稀疏表示通常用来作为先验正则项，所以为了求解图像重构问题的最优解，增加了一个重构逼近项的约束优化问题。在模型的基础上，增加了重构的残差项，进而形成了光滑L0范数稀疏表示加误差逼近的完整的压缩感知图像重构模型：

$$
\mathrm{J}(\mathrm{s})=\arg \min \left{\mathrm{F}{\sigma}(\mathrm{s})+\lambda|\Theta \mathrm{s}-\mathrm{y}|{2}^{2}\right}
$$

式中λ为权重平衡参数。上述模型中的稀疏表示先验项以稀疏信号s为处理内容，而误差逼近项是重构图像投影测量值与实际值的残差极小为整幅图像的全局优化。由于是逼近最小L0范数的算法求解方式，$\mathrm{F}{\sigma}(\mathrm{s})$中的σ参数是逐步减小的，因此逼近最小L0范数稀疏表示的优化问题可采用梯度下降的迭代计算方式来解决，即对式中s进行求导，在梯度方向$\Delta \mathrm{J}(\mathrm{s})$进行迭代：
$$
\Delta J(s)=d=2 \lambda\left(\Theta^{T}(\Theta s-y)\right)+\left(1 / \sigma^{2}\right)^{*}\left[s{1} \exp \left(-s{i}^{2} \sigma^{2} / 2\right), \cdots, s{m} \exp \left(-s_{m}^{2} \sigma^{2} / 2\right)\right]^{\mathrm{T}}
$$

基于梯度下降的逼近最小L0范数算法的伪代码过程如下：
1) 变量名的物理含义：传感矩阵Θ，测量值y，原始信号的稀疏表示s；
2) 初始化：$\hat{\mathrm{s}}{0}=\Theta^{\perp} \mathrm{y}, \Theta^{\perp}$为Θ的伪逆阵，$\Theta^{\perp}=\left(\Theta^{\mathrm{T}} \Theta\right)^{-1} \Theta^{\mathrm{T}}$，设增长序列$\sigma=\left[\sigma{1}, \ldots, \sigma{k}\right]$、λ和梯度下降步长µ；
3) 循环σ序列k：
a) $\sigma=\sigma{\mathrm{k}}, \quad \widehat{\mathrm{s}}=\widehat{\mathrm{s}}{\mathrm{k}-1}$；
b) 梯度下降法迭代L次：
a.梯度下降方向：
$$
d=2 \lambda\left(\Theta^{T}(\Theta s-y)\right)+\left(1 / \sigma^{2}\right)^{*}\left[s{1} \exp \left(-s{i}^{2} \sigma^{2} / 2\right), \cdots, s{m} \exp \left(-S{m}^{2} \sigma^{2} / 2\right)\right]^{\mathrm{T}}
$$b.梯度方向更新：$\widehat{\mathrm{s}}=\widehat{\mathrm{s}}-\mu \mathrm{d}$
c.约束正交投影：$\widehat{\mathrm{s}}=\widehat{\mathrm{s}}-\Theta^{\perp}(\Theta \widehat{\mathrm{s}}-\mathrm{y})$
c) $\hat{\mathrm{s}}{\mathrm{k}}=\hat{\mathrm{s}}, \quad \mathrm{k}=\mathrm{k}+1$；
4) 结束循环：输出$\widehat{\mathrm{s}}_{\mathrm{k}}$为最终求得的稀疏信号s。

值得注意的是，在算法求解过程中，随着迭代的进行，误差$|\Theta s-y|{2}^{2}$越来越小，使稀疏表示正则项$\mathrm{F}{\sigma}(\mathrm{s})=\sum{\mathrm{i}=1}^{\mathrm{n}} \mathrm{f}{\sigma}\left(\mathrm{s}{\mathrm{i}}\right)$越来越逼近于L0范数，同时可以让误差$|\Theta s-y|{2}^{2}$所占权重逐渐减小，逼近L0范数的稀疏表示正则项$\mathrm{F}{\sigma}(\mathrm{s})$所占权重逐渐增大。因此λ设置成递增序列，σ设置成递减序列，从而在迭代过程中可以根据误差$|\Theta s-y|{2}^{2}$的大小能够自适应地去调整λ和σ的值。

由逼近图可知，随着σ的逐渐减小，稀疏表示正则项$\mathrm{F}_{\sigma}(\mathrm{s})$越来越光滑地逼近近似的L0范数，在逐渐降低了重构稀疏信号的稀疏度的同时以更大效率逐渐逼近我们需要的全局最优解，从而降低了算法复杂度和提高了图像重建的精度，从而更加利于高分辨率图像的重建。

简单的实验验证

例如对于矩阵的乘法$\Theta^{*} \mathrm{s}=y$

$$
\left[\begin{array}{ccccc}{0.245} & {-0.053} & {4.478} & {-3.426} & {3.487} & {2.85} \ {-1.205} & {3.42} & {-2.39} & {-2.752} & {4.49} & {11.676}\end{array}\right] *\left[\begin{array}{c}{0} \ {0} \ {-8.361} \ {0} \ {0} \ {-12.982}\end{array}\right]=\left[\begin{array}{c}{-74.4383} \ {-131.5950}\end{array}\right]
$$

其中Θ为传感矩阵，s为原始信号在变换域中的稀疏表达，y为维度更低的测量信号。由2维的y来恢复6维的s是一个欠定问题，但如测量矩阵满足RIP条件，将s的稀疏性作为正则项，就有较大的概率重构出原始信号。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
int main(int argc, char **argv)
{
	Mat A = (Mat_<double>(2, 6) << 0.245f, -0.053f, 4.478f, -3.426f, 3.487f, 2.85f, 
							       -1.205f, 3.42f, -2.39f, -2.752f, 4.49f, 11.676f);
	Mat y = (Mat_<double>(2, 1) << -74.4393f, -131.5950f);
	Mat s(6, 1, CV_64FC1);
	/*sigma_min为sigma迭代终止的条件，sigma_min越小，重构结果就越精确*/
	double sigma_min = 0.001;
	if (SL0(A, y, sigma_min, s))
	{
		cout << "传感矩阵为：" << endl << A << endl << endl;;
		cout << "生成向量为：" << endl << y << endl << endl;;
		cout << "重构出更高维的原始向量的稀疏表示为：" << endl << s << endl << endl;;
	}
	else
		fprintf(stderr, "Error in calling SL0 function.\nsigma_min = %f", 
				sigma_min);
	system("pause");
	return 0;
}

验证结果为

SL0算法的实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
bool MyDelta(const Mat s, const double sigma, Mat &delta)
{
	for (int i = 0; i < s.rows; ++i)
		for (int j = 0; j < s.cols; ++j)
			delta.at<double>(i, j)  = s.at<double>(i, j) *
			exp(-1 * s.at<double>(i, j) * s.at<double>(i, j) / (2 * sigma * sigma));
	return true;
}
bool SL0(const Mat A, const Mat y, const double sigma_min, Mat &s)
{
	double iteration_factor = 0.5;
	double lamdba = 0.1;
	double sigma;
	int mu_0 = 2;
	int L = 3;
	
	//A_pinv为矩阵A的伪逆
	Mat A_pinv;
	invert(A, A_pinv, cv::DECOMP_SVD);
	
	//向量s和增长序列sigma的初始化
	s = A_pinv * y;
	double minVal, maxVal;
	minMaxLoc(s, &minVal, &maxVal);
	sigma = abs(minVal) > abs(maxVal) ? abs(minVal) : abs(maxVal);
	sigma *= 2;
	//到达指定的停机准则就退出，siama_min越小，重建结果越精确
	while (sigma > sigma_min)
	{
		//迭代下降法迭代L次
		for (int i = 0; i != L; ++ i)
		{
			//delta为梯度下降方向
			Mat delta = Mat::zeros(s.size(), s.type());
			if (MyDelta(s, sigma, delta))
				delta += 2 * lamdba * (A.t() * (A * s - y));
			else
				fprintf(stderr, "Error in calculating delta");
			s = s - mu_0 * delta;		//梯度方向更新
			s = s - A_pinv * (A * s - y);		//约束正交投影
		}
		//迭代过程中，误差所占权重逐渐减小，自适应的调整lamdba和siama的值
		sigma = sigma * iteration_factor;
		lamdba = lamdba * (1 / iteration_factor);
	}
	return true;
}

矩阵伪逆的opencv实现

1、矩阵的逆
定义：
设A是数域上的一个n阶方阵，若在相同数域上存在另一个n阶矩阵B，使得： AB=BA=I。 则我们称B是A的逆矩阵，而A则被称为可逆矩阵。

可逆条件：
A是可逆矩阵的充分必要条件是，即可逆矩阵就是非奇异矩阵。（当 时，A称为奇异矩阵）

性质：
-矩阵A可逆的充要条件是A的行列式不等于0。
-可逆矩阵一定是方阵。
-如果矩阵A是可逆的，A的逆矩阵是唯一的。
-可逆矩阵也被称为非奇异矩阵、满秩矩阵。
-两个可逆矩阵的乘积依然可逆。
-可逆矩阵的转置矩阵也可逆。
-矩阵可逆当且仅当它是满秩矩阵。

2、矩阵的伪逆
伪逆矩阵：
伪逆矩阵是逆矩阵的广义形式。由于奇异矩阵或非方阵的矩阵不存在逆矩阵，但在matlab里可以用函数pinv(A)求其伪逆矩阵。基本语法为X=pinv(A),X=pinv(A,tol),其中tol为误差,pinv为pseudo-inverse的缩写：max(size(A))norm(A)eps。函数返回一个与A的转置矩阵A’ 同型的矩阵X，并且满足：AXA=A,XAX=X.此时，称矩阵X为矩阵A的伪逆，也称为广义逆矩阵。pinv(A)具有inv(A)的部分特性，但不与inv(A)完全等同。如果A为非奇异方阵，pinv(A)=inv(A)，但却会耗费大量的计算时间，相比较而言，inv(A)花费更少的时间。

3、opencv函数invert（）介绍

double invert(InputArray src, OutputArraydst, int flags=DECOMP_LU);

功能：用以求取一个矩阵的逆或者伪逆。

src： 输入，浮点型（32位或者64位）的M×N的矩阵，当参数3的使用方法为DECOMP_CHOLESKY DECOMP_LU DECOMP_EIG时函数功能为求逆，此时需保证M=N（参见参数flag）。

dst： 输出，与输入矩阵类型一致的N×M的矩阵。

flag：求逆方法，提供4种可选择的方法：DECOMP_CHOLESKY（基于CHOLESKY分解的方法）， DECOMP_LU（基于LU分解的方法）， DECOMP_EIG（基于特征值分解的方法）， DECOMP_SVD（基于奇异值分解的方法）。其中，前三种方法要求输入的矩阵必须为方阵，此时计算结果为矩阵的逆；最后一种方法为对非方阵的伪逆计算，对矩阵的形状没有要求。函数接口的默认参数为DECOMP_LU方法（应该是效率较高的一种方法）。

4、实现代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#include <iostream>
#include <vector>
#include <opencv2/opencv.hpp>
using namespace std;
using namespace cv;
int main()
{
	vector<vector<float>> vec{ { 0.68f, 0.597f, -0.211f },
						{ 0.823f, 0.566f, -0.605f } };
	const int rows{ 2 }, cols{ 3 };
	printf("source matrix:\n");
	for (int i = 0; i < vec.size(); i++)
	{
		for (int j = 0; j < vec[i].size(); j++)
			cout << vec[i][j] << '\t';
		cout << endl;
	}
	
	printf("\nopencv implement pseudoinverse:\n");
	Mat mat(rows, cols, CV_32FC1);
	for (int y = 0; y < rows; ++y) {
		for (int x = 0; x < cols; ++x) {
			mat.at<float>(y, x) = vec[y][x];
		}
	}
	Mat pinv2;
	invert(mat, pinv2, cv::DECOMP_SVD);
	
	for (int y = 0; y < pinv2.rows; ++y)
	{
		for (int x = 0; x < pinv2.cols; ++x)
			cout << pinv2.at<float>(y, x) << '\t';
		cout << endl;
	}
	system("pause");
	return 0;
}

分水岭算法及其实现

１ - 算法描述
1.1 分水岭算法的原理
　　分水岭的概念是以三维方式来形象化一幅图像为基础的：两个空间坐标再加上强度。在这种“地形学”解释中，考虑三种类型的点：（a）局部最小值点，该点对应一个盆地的最低点，当我们在盆地里滴一滴水的时候，由于重力作用，水最终会汇聚到该点。注意：可能存在一个最小值面，该平面内的都是最小值点。（b）盆地的其它位置点，该位置滴的水滴会汇聚到局部最小点。（c）盆地的边缘点，是该盆地和其它盆地交接点，在该点滴一滴水，会等概率的流向任何一个盆地。对于一个特定的区域最小值，满足条件（b）的点的集合称为该最小值的汇水盆地或分水岭。满足条件（c）的点形成地面的锋线，称为分割线或分水线。


图1 分水岭算法中三种类型的点

  基于这些概念的分割算法的主要目标是找出分水线。假设在每个区域的最小值上打一个洞，并且让水通过洞以均匀的速率上升，从低到高淹没整个地形。当不同汇水盆地中上升的水聚集时，修建一个水坝来阻止这种聚合。水将达到在水线上只能见到各个水坝的顶部的程度。这些大坝的边界对应于分水岭的分割线。因此，他们是由分水岭算法提取出来的边界。

１.2 水坝构建
　　在实际操作中，水位的上升常利用灰度0-255的分层处理来实现，当处理到第N层时，灰度小于N的像素就会被淹没。对于大多数情况，聚水盆的生长速度是不定的，也许在第N-1次处理时，两个聚水盆之间还有相当的距离，而在第N次处理中他们就已经相互连通了，这为堤坝的建立带来了一定的难度。


图2 水坝构造示意图

  图2中a显示了两个相邻的聚水盆，其中颜色越深的区域海拔越低，当水位逐渐上升到第N-1阶段时，黑色区域首先被淹没，如图2中b所示，这时图中存在两个彼此分开的水域，他们之间还有相当一段距离。然而当水面继续上升至第N阶段时，图中的水域已经连通，如图2中c所示。
  为了在水域之间适当的位置构建堤坝防止水域的连通，可以借助图像形态学处理中的膨胀方法。由于已知两个水域会在第N-1阶段至第N阶段间连通，回到第N-1阶段的图像如图2中b所示，尝试对图中不同水域分别进行膨胀操作，并在膨胀中寻找两个区域的交点。如图2中d所示为两个水域分别进行一次膨胀操作的结果，从结果中可以看出两个区域仍然没有交点，因此还不能形成堤坝。继续对两个水域进行膨胀操作，在对右边区域的膨胀操作中，两个区域出现了交点，标记这些交点并在交点处修建堤坝，如图2中e所示，其中堤坝位置用交叉图案表示。继续进行膨胀操作，并不断标记新的交点，最终就能够得到完整的水坝轮廓。值得注意的是，在算法实现中，应对每次膨胀结果与第N阶段的淹没区域进行比较，保证堤坝位置在第N阶段的淹没区域之内，算法一直进行直到经过膨胀的水域完全覆盖第N阶段中淹没区域为止。通过以上步骤就能够获得较为准确的堤坝位置，如图2中f所示。
　　由于实际中处理的图像往往比较复杂，图中灰度极小值点较多，这就导致聚水盆数目过多，从而造成图像的过度分割，比如下面的图3，这样的分割效果是毫无用处的。解决图像的过度分割问题最直接的方法就是进行区域的合并，如果在分水岭算法中定义聚水盆的最小深度，就能比较好地解决过度分割的问题。


图3 过度分割的结果

  聚水盆的深度可以定义为盆地中海拔最低点与当前水面高度的差值，在计算时可以用当前已淹没的灰度值减去区域的灰度最小值来得到。当水位上升造成两个不同的水域连接时，判断两个水域所在聚水盆的深度是否都大于定义的最小深度，若是则在两个区域间构建堤坝，否则对两个水域进行合并。运用这样的方法，就可以通过改变最小深度的值来控制图像分割的程度，从而获得图像最有效的分割结果。

２ - 分水岭算法的OpenCv实现
　　函数watershed实现的分水岭算法是基于标记的分割算法中的一种。在把图像传给函数之前，需要大致勾画标记出图像中的期望进行分割的区域，他们被标记为正指数。所以，每一个区域都会被标记为像素值1、2、3等，表示成为一个或者多个连接组件。这些标记的值可以使用findContours()函数和drawContours()函数由二进制的掩码检测出来。不难理解，这些标记就是即将绘制出来的分割区域的“种子”，而没有标记清楚的区域，被置为0。在函数输出中，每一个标记中的像素被设置为“种子”的值，而区域间的值被设置为-1。
　　C++: void watershed(InputArray image, InputOutputArray markers)
　　第一个参数，InputArray类型的src，输入图像，即原图像，填Mat类的对象即可，且需为8位三通道的彩色图像。
　　第二个参数，InputOutputArray类型的markers，函数调用后的运算结果存在这里，输入/输出32位单通道图像的标记结果。即这个参数用于存放函数调用后的输出结果，需和原图像有一样的尺寸和类型。


图4 原始图

  编写程序，对图4进行基于分水岭算法的图像分割。


图5 期望分割的区域

  直接使用图像处理工具箱的提供的现成函数进行分割往往很难达到预期效果，而预先勾画标记出图像中期望进行分割的区域，在应用分水岭算法就会取得较好的分割效果。如图5所示，其中左图显示了原始图中期望进行分割的各个部分，右图更清楚的显示了期望分割部分的空间位置。


图6 分割结果

  图6显示了分割的结果。分割的各个部分使用随机的颜色进行表示；若存在没有标记清楚的区域，则将此区域用纯黑色表示；用于划分区域的线条，即被设置为“种子”的像素值，用纯白色表示。左图显示了原始图像的分割结果，右图显示了分割结果和原始图按等比例进行混合的效果。

3 - OpenCv代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <iostream>
using namespace cv;
using namespace std;
//-----------------------------------【宏定义部分】-------------------------------------------- 
//  描述：定义一些辅助宏 
//------------------------------------------------------------------------------------------------ 
#define WINDOW_NAME1 "【程序窗口1】"        //为窗口标题定义的宏 
#define WINDOW_NAME2 "【分水岭算法效果图】"        //为窗口标题定义的宏
//-----------------------------------【全局函变量声明部分】--------------------------------------
//		描述：全局变量的声明
//-----------------------------------------------------------------------------------------------
Mat g_maskImage, g_srcImage;
Point prevPt(-1, -1);
static void on_Mouse(int event, int x, int y, int flags, void*);
int main(int argc, char** argv)
{
	//【0】改变console字体颜色
	system("color 6F");
	//【1】载入原图并显示，初始化掩膜和灰度图
	g_srcImage = imread("1.jpg", 1);
	imshow(WINDOW_NAME1, g_srcImage);
	Mat srcImage, grayImage;
	g_srcImage.copyTo(srcImage);
	cvtColor(g_srcImage, g_maskImage, COLOR_BGR2GRAY);
	cvtColor(g_maskImage, grayImage, COLOR_GRAY2BGR);
	g_maskImage = Scalar::all(0);
	//【2】设置鼠标回调函数
	setMouseCallback(WINDOW_NAME1, on_Mouse, 0);
	//【3】轮询按键，进行处理
	while (1)
	{
		//获取键值
		int c = waitKey(0);
		//若按键键值为ESC时，退出
		if ((char)c == 27)
			break;
		//按键键值为2时，恢复源图
		if ((char)c == '2')
		{
			g_maskImage = Scalar::all(0);
			srcImage.copyTo(g_srcImage);
			imshow("image", g_srcImage);
		}
		//若检测到按键值为1或者空格，则进行处理
		if ((char)c == '1' || (char)c == ' ')
		{
			//定义一些参数
			int i, j, compCount = 0;
			vector<vector<Point> > contours;
			vector<Vec4i> hierarchy;
			//寻找轮廓
			findContours(g_maskImage, contours, hierarchy, CV_RETR_CCOMP, CV_CHAIN_APPROX_SIMPLE);
			//轮廓为空时的处理
			if (contours.empty())
				continue;
			//拷贝掩膜
			Mat maskImage(g_maskImage.size(), CV_32S);
			maskImage = Scalar::all(0);
			//循环绘制出轮廓    //感觉这里绘制轮廓不重要，重要的是几个 因为后面又覆盖了maskImage watershed( srcImage, maskImage );
			for (int index = 0; index >= 0; index = hierarchy[index][0], compCount++)
				drawContours(maskImage, contours, index, Scalar::all(compCount + 1), -1, 8, hierarchy, INT_MAX);
			//compCount为零时的处理
			if (compCount == 0)
				continue;
			//生成随机颜色
			vector<Vec3b> colorTab;
			for (i = 0; i < compCount; i++)
			{
				int b = theRNG().uniform(0, 255);
				int g = theRNG().uniform(0, 255);
				int r = theRNG().uniform(0, 255);
				colorTab.push_back(Vec3b((uchar)b, (uchar)g, (uchar)r));
			}
			//计算处理时间并输出到窗口中
			double dTime = (double)getTickCount();
			watershed(srcImage, maskImage);
			dTime = (double)getTickCount() - dTime;
			printf("\t处理时间 = %gms\n", dTime*1000. / getTickFrequency());
			//双层循环，将分水岭图像遍历存入watershedImage中
			Mat watershedImage(maskImage.size(), CV_8UC3);
			for (i = 0; i < maskImage.rows; i++)
				for (j = 0; j < maskImage.cols; j++)
				{
					int index = maskImage.at<int>(i, j);           //不同的index用不同的index颜色
					if (index == -1)
						watershedImage.at<Vec3b>(i, j) = Vec3b(255, 255, 255);
					else if (index <= 0 || index > compCount)
						watershedImage.at<Vec3b>(i, j) = Vec3b(0, 0, 0);
					else
						watershedImage.at<Vec3b>(i, j) = colorTab[index - 1];
				}
			//混合灰度图和分水岭效果图并显示最终的窗口
			watershedImage = watershedImage * 0.5 + grayImage * 0.5;
			imshow(WINDOW_NAME2, watershedImage);
		}
	}
	return 0;
}
static void on_Mouse(int event, int x, int y, int flags, void*)
{
	//处理鼠标不在窗口中的情况
	if (x < 0 || x >= g_srcImage.cols || y < 0 || y >= g_srcImage.rows)
		return;
	//处理鼠标左键相关消息
	if (event == CV_EVENT_LBUTTONUP || !(flags & CV_EVENT_FLAG_LBUTTON))
		prevPt = Point(-1, -1);
	else if (event == CV_EVENT_LBUTTONDOWN)
		prevPt = Point(x, y);
	//鼠标左键按下并移动，绘制出白色线条
	else if (event == CV_EVENT_MOUSEMOVE && (flags & CV_EVENT_FLAG_LBUTTON))
	{
		Point pt(x, y);
		if (prevPt.x < 0)
			prevPt = pt;
		line(g_maskImage, prevPt, pt, Scalar::all(255), 5, 8, 0);
		line(g_srcImage, prevPt, pt, Scalar::all(0), 5, 8, 0);
		prevPt = pt;
		imshow(WINDOW_NAME1, g_srcImage);
	}
}

Canny算子c++实现

实现一个canny算子，opencv版本为3.4.5

在这里插入图片描述

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
#include <iostream>
#include "opencv2/opencv.hpp"
using namespace std;
using namespace cv;
/*
生成高斯卷积核 kernel
*/
void Gaussian_kernel(int kernel_size, int sigma, Mat &kernel)
{
	const double PI = 3.1415926;
	int m = kernel_size / 2;
	kernel = Mat(kernel_size, kernel_size, CV_32FC1);
	float s = 2 * sigma * sigma;
	//笔记： 注意二维高斯函数下面就是2*pi*sigma^2，没有根号
	for (int i = 0; i < kernel_size; i++)
	{
		for (int j = 0; j < kernel_size; j++)
		{
			int x = i - m;
			int y = j - m;
			kernel.at<float>(i, j) = exp(-(x * x + y * y) / s) / (PI * s);
		}
	}
}
/*
计算梯度值和方向
imageSource 原始灰度图
imageX X方向梯度图像
imageY Y方向梯度图像
gradXY 该点的梯度幅值
theta 梯度方向角度 theta=arctan(imageY/imageX)
*/
void GradDirection(const Mat imageSource, Mat &imageX, Mat &imageY, Mat &gradXY, Mat &theta)
{
	
	imageX = Mat::zeros(imageSource.size(), CV_32SC1);
	imageY = Mat::zeros(imageSource.size(), CV_32SC1);
	gradXY = Mat::zeros(imageSource.size(), CV_32SC1);
	theta = Mat::zeros(imageSource.size(), CV_32SC1);
	int rows = imageSource.rows;
	int cols = imageSource.cols;
	/*
	Mat.step参数指图像的一行实际占用的内存长度，以字节为基本单位，
	因为opencv中的图像会对每行的长度自动补齐（8的倍数），
	编程时尽量使用指针，指针读写像素是速度最快的，使用at函数最慢。
	*/
	int stepXY = imageX.step;
	int step = imageSource.step;
	/*
	Mat::data的默认类型为uchar*，但很多时候需要处理其它类型，如float、int，
	此时需要将data强制类型转换，如：
	Mat src(1000,1000,CV_32F);
	float* myptr = (float*)src.data;
	无论Mat的type是何种类型，Mat::data均为uchar*
	*/
	uchar *PX = imageX.data;
	uchar *PY = imageY.data;
	uchar *P = imageSource.data;
	uchar *XY = gradXY.data;
	for (int i = 1; i < rows - 1; i++)
	{
		for (int j = 1; j < cols - 1; j++)
		{
			int a00 = P[(i - 1)*step + j - 1];
			int a01 = P[(i - 1)*step + j];
			int a02 = P[(i - 1)*step + j + 1];
			int a10 = P[i*step + j - 1];
			int a11 = P[i*step + j];
			int a12 = P[i*step + j + 1];
			int a20 = P[(i + 1)*step + j - 1];
			int a21 = P[(i + 1)*step + j];
			int a22 = P[(i + 1)*step + j + 1];
			double gradX = double(a02 + 2 * a12 + a22 - a00 - 2 * a10 - a20);
			double gradY = double(a00 + 2 * a01 + a02 - a20 - 2 * a21 - a22);
			//遍历是从i,j=1开始到rows,cols-2的，图像像素都初始化为0，
			//所以最外边的一圈像素为0是黑的
			//PX[i*stepXY + j*(stepXY / step)] = abs(gradX);
			//PY[i*stepXY + j*(stepXY / step)] = abs(gradY);
			imageX.at<int>(i, j) = abs(gradX);
			imageY.at<int>(i, j) = abs(gradY);
			if (gradX == 0)
			{
				gradX = 0.000000000001;
			}
			//弧度数改为角度数,除以pi乘以180
			theta.at<int>(i, j) = atan(gradY / gradX) * 57.3;
			theta.at<int>(i, j) = (theta.at<int>(i, j) + 360) % 360;
			gradXY.at<int>(i, j) = sqrt(gradX*gradX + gradY * gradY);
			//XY[i * stepXY + j * (stepXY / step)] = sqrt(gradX * gradX + gradY * gradY);
		}
	}
	/*
	在经过处理后，需要用convertScaleAbs()函数将其转回原来的uint8形式，
	否则将无法显示图像，而只是一副灰色的窗口。 
	函数原型为
	void convertScaleAbs(InputArray src, OutputArray dst,
                                  double alpha = 1, double beta = 0);
	其中可选参数alpha是伸缩系数，beta是加到结果上的一个值，结果返回uint8类型的图片
	功能：实现将原图片转换为uint8类型
	*/
	convertScaleAbs(imageX, imageX);
	convertScaleAbs(imageY, imageY);
	convertScaleAbs(gradXY, gradXY);
}
/*
局部非极大值抑制
沿着该点梯度方向，比较前后两个点的幅值大小，若该点大于前后两点，则保留，
若该点小于前后两点任意一点，则置为0；
imageInput 输入得到梯度图像
imageOutput 输出的非极大值抑制图像
theta 每个像素点的梯度方向角度
imageX X方向梯度
imageY Y方向梯度
*/
void NonLocalMaxValue(const Mat imageInput, Mat &imageOutput, const Mat &theta, const Mat &imageX, const Mat &imageY)
{
	imageOutput = imageInput.clone();
	int cols = imageInput.cols;
	int rows = imageInput.rows;
	for (int i = 1; i < rows - 1; i++)
	{
		for (int j = 1; j < cols - 1; j++)
		{
			//执行continue时，会执行本次迭代的剩余部分，并开始下一次迭代
			//执行break时，会执行包含它的循环，并执行下一阶段
			if (0 == imageInput.at<uchar>(i, j))continue;
			int g00 = imageInput.at<uchar>(i - 1, j - 1);
			int g01 = imageInput.at<uchar>(i - 1, j);
			int g02 = imageInput.at<uchar>(i - 1, j + 1);
			int g10 = imageInput.at<uchar>(i, j - 1);
			int g11 = imageInput.at<uchar>(i, j);
			int g12 = imageInput.at<uchar>(i, j + 1);
			int g20 = imageInput.at<uchar>(i + 1, j - 1);
			int g21 = imageInput.at<uchar>(i + 1, j);
			int g22 = imageInput.at<uchar>(i + 1, j + 1);
			int direction = theta.at<int>(i, j); //该点梯度的角度值
			int g1 = 0;
			int g2 = 0;
			int g3 = 0;
			int g4 = 0;
			double tmp1 = 0.0; //保存亚像素点插值得到的灰度数
			double tmp2 = 0.0;
			//在C中，有abs(),labs(),fabs()分别计算int ,long ,double 类型的绝对值；
			//而在C++中，abs()已被重载，可适用于各种类型
			double weight = fabs((double)imageY.at<uchar>(i, j) / (double)imageX.at<uchar>(i, j));
			//去掉也可以，如果weight==0，那direction==0或180就不用差值，直接用原来的值
			if (weight == 0)
				weight = 0.0000001;
			/*
			关于这些公式的含义
			https://www.cnblogs.com/love6tao/p/5152020.html 有解释
			形如g10 * (1 - weight) + g20 * (weight)都是用两边的像素计算得到的亚像素
			不过他w和1-w应该写反了
			*/
			if (weight > 1)
			{
				weight = 1 / weight;
			}
			if ((0 <= direction && direction < 45) || 180 <= direction && direction < 225)
			{
				tmp1 = g10 * (1 - weight) + g20 * (weight);
				tmp2 = g02 * (weight)+g12 * (1 - weight);
			}
			if ((45 <= direction && direction < 90) || 225 <= direction && direction < 270)
			{
				tmp1 = g01 * (1 - weight) + g02 * (weight);
				tmp2 = g20 * (weight)+g21 * (1 - weight);
			}
			if ((90 <= direction && direction < 135) || 270 <= direction && direction < 315)
			{
				tmp1 = g00 * (weight)+g01 * (1 - weight);
				tmp2 = g21 * (1 - weight) + g22 * (weight);
			}
			if ((135 <= direction && direction < 180) || 315 <= direction && direction < 360)
			{
				tmp1 = g00 * (weight)+g10 * (1 - weight);
				tmp2 = g12 * (1 - weight) + g22 * (weight);
			}
			if (imageInput.at<uchar>(i, j) < tmp1 || imageInput.at<uchar>(i, j) < tmp2)
			{
				imageOutput.at<uchar>(i, j) = 0;
			}
		}
	}
}
/*
双阈值的机理是：
指定一个低阈值A，一个高阈值B，选取占直方图总数70%为B，且B为1.5到2倍大小的A；
灰度值小于A的，置为0,灰度值大于B的，置为255；
*/
void DoubleThreshold(Mat &imageInput, const double lowThreshold, const double highThreshold)
{
	int cols = imageInput.cols;
	int rows = imageInput.rows;
	for (int i = 0; i < rows; i++)
	{
		for (int j = 0; j < cols; j++)
		{
			double temp = imageInput.at<uchar>(i, j);
			temp = temp > highThreshold ? (255) : (temp);
			temp = temp < lowThreshold ? (0) : (temp);
			imageInput.at<uchar>(i, j) = temp;
		}
	}
}
/*
连接处理:
灰度值介于A和B之间的，考察该像素点临近的8像素是否有灰度值为255的，
若没有255的，表示这是一个孤立的局部极大值点，予以排除，置为0；
若有255的，表示这是一个跟其他边缘有“接壤”的可造之材，置为255，
之后重复执行该步骤，直到考察完最后一个像素点。
其中的邻域跟踪算法，从值为255的像素点出发找到周围满足要求的点，把满足要求的点设置为255，
然后修改i,j的坐标值，i,j值进行回退，在改变后的i,j基础上继续寻找255周围满足要求的点。
当所有连接255的点修改完后，再把所有上面所说的局部极大值点置为0；（算法可以继续优化）。
参数1，imageInput：输入和输出的梯度图像
参数2，lowTh:低阈值
参数3，highTh:高阈值
*/
void DoubleThresholdLink(Mat &imageInput, double lowTh, double highTh)
{
	int cols = imageInput.cols;
	int rows = imageInput.rows;
	for (int i = 1; i < rows - 1; i++)
	{
		for (int j = 1; j < cols - 1; j++)
		{
			double pix = imageInput.at<uchar>(i, j);
			if (pix != 255)continue;
			bool change = false;
			for (int k = -1; k <= 1; k++)
			{
				for (int u = -1; u <= 1; u++)
				{
					if (k == 0 && u == 0)continue;
					double temp = imageInput.at<uchar>(i + k, j + u);
					if (temp >= lowTh && temp <= highTh)
					{
						imageInput.at<uchar>(i + k, j + u) = 255;
						change = true;
					}
				}
			}
			//如果要连接像素值，则如果连接的是上边三个或者左边紧挨着的一个
			//可能会对之前的像素产生影响，即可能前面的邻接像素也要改变
			//这里采取的方法是不管要改变哪个像素，都将当前位置返回到左上角
			//最上面一条边和最左边一条边不能返回到左上角，就分别左移和上移一格
			//注意j马上还要+1，i不用+1
			//至于为什么是i>1而j>2，其实j>1也可以，只是因为j是减2，而2-2=0，不过还要+1所以没问题
			if (change)
			{
				if (i > 1)i--;
				if (j > 2)j -= 2;
			}
		}
	}
	for (int i = 0; i < rows; i++)
	{
		for (int j = 0; j < cols; j++)
		{
			if (imageInput.at<uchar>(i, j) != 255)
			{
				imageInput.at<uchar>(i, j) = 0;
			}
		}
	}
}
int main()
{
	//Mat image = imread("lena.png");
	//imshow("origin image", image);
	//转换为灰度图
	//Mat grayImage;
	//cvtColor(image, grayImage, CV_RGB2GRAY);
	Mat grayImage = imread("I.bmp", 0);
	imshow("gray image", grayImage);
	//高斯滤波
	Mat gausKernel;
	int kernel_size = 5;
	double sigma = 1;
	Gaussian_kernel(kernel_size, sigma, gausKernel);
	Mat gausImage;
	//Convolves an image with kernel 即利用内核实现对图像的卷积运算
	filter2D(grayImage, gausImage, grayImage.depth(), gausKernel);
	imshow("gaus image", gausImage);
	imwrite("高斯滤波后的图像.bmp", gausImage);
	//计算XY方向梯度
	Mat imageX, imageY, imageXY;
	//theta为梯度方向，theta=arctan(imageY/imageX)
	Mat theta;
	GradDirection(gausImage, imageX, imageY, imageXY, theta);
	imshow("XY grad", imageXY);
	imwrite("梯度模值.bmp", imageXY);
	//对梯度幅值进行非极大值抑制
	Mat localImage;
	NonLocalMaxValue(imageXY, localImage, theta, imageX, imageY);
	imshow("Non local maxinum image", localImage);
	imwrite("非极大性值抑制.bmp", localImage);
	//双阈值算法检测和边缘连接
	DoubleThreshold(localImage, 60, 100);
	DoubleThresholdLink(localImage, 60, 100);
	imshow("canny image", localImage);
	imwrite("双阈值和边缘连接.bmp", localImage);
	Mat temMat;
	/*
	void Canny( InputArray image, OutputArray edges,
                         double threshold1, double threshold2,
                         int apertureSize = 3, bool L2gradient = false );
	输入图像 输出图像 低阈值 高阈值 算子大小 是否采用更精确的方式计算图像梯度
	void Canny( InputArray dx, InputArray dy,
						 OutputArray edges,
						 double threshold1, double threshold2,
						 bool L2gradient = false );
	输入图像x,y方向的导数 输出图像 低阈值 高阈值 算子大小 是否采用更精确的方式计算图像梯度
	*/
	Canny(grayImage, temMat, 60, 100);
	imshow("opencv canny image", temMat);
	imwrite("opencv自带的canny算法得到的效果图.bmp", temMat);
	waitKey(0);
	return 0;
}

C++实现均值滤波器和中值滤波器

C++实现均值滤波器和中值滤波器

代码实现均值滤波器和中值滤波器
在这里插入图片描述
由于中值滤波器是非线性滤波，不是卷积，所以均值和中值滤波分开实现。opencv版本为3.4.5

my_convolution.h

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#ifndef MY_CONVOLUTION
#define MY_CONVOLUTION
#include <opencv2/opencv.hpp>
class My_Convolution {
public:
	My_Convolution();
	~My_Convolution();
	bool load_kernal(const cv::Mat kernal);//加载卷积核
	void convolute(const cv::Mat &image, cv::Mat &dst);//卷积操作
private:
	bool kernal_loaded;//是否已经加载卷积核
	cv::Mat curr_kernal;//当前卷积核
	int bios_x, bios_y;//记录偏移量
	//计算每一个像素的掩模乘积之和
	void compute_sum_of_product(int i, int j, int chan, cv::Mat &complete_image, cv::Mat & dst);
	//将原图像转换成边框补全的图像
	void complete_image_transform(const cv::Mat &image, cv::Mat &dst);
};
#endif // MY_CONVOLUTION

my_convolution.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
#include "my_convolution.h"
using namespace std;
using namespace cv;
My_Convolution::My_Convolution() {
	kernal_loaded = false;
}
My_Convolution::~My_Convolution() {}
//加载卷积核
bool My_Convolution::load_kernal(const Mat kernal) {
	if (kernal.cols % 2 == 1 && kernal.rows % 2 == 1) {
		curr_kernal = kernal.clone();
		bios_x = (kernal.cols - 1) / 2;
		bios_y = (kernal.rows - 1) / 2;
		kernal_loaded = true;
		return true;
	}
	else {
		cout << "The size of kernal is not suitable!" << endl;
		return false;
	}
}
//卷积操作
void My_Convolution::convolute(const Mat &image, Mat &dst) {
	if (!kernal_loaded) {
		cout << "kernal is empty!Please load the kernal first!" << endl;return;
	}
	Mat complete_image;
	complete_image_transform(image, complete_image);
	dst = Mat::zeros(image.rows, image.cols, image.type());
	int channels = image.channels();//获取图像的通道数
	if (channels == 3) {
		for (int chan = 0;chan < channels;chan++) {
			for (int i = 0;i < dst.rows;i++) {
				for (int j = 0;j < dst.cols;j++) {
					compute_sum_of_product(i, j, chan, complete_image, dst);
				}
			}
		}
		return;
	}
	if (channels == 1) {
		for (int i = 0;i < dst.rows;i++) {
			for (int j = 0;j < dst.cols;j++) {
				compute_sum_of_product(i, j, 0, complete_image, dst);
			}
		}
	}
}
//计算掩模乘积之和
void My_Convolution::compute_sum_of_product(int i, int j, int chan, Mat &complete_image, Mat &dst) {
	if (complete_image.channels() == 3) {
		float sum = 0;
		int bios_rows = i;
		int bios_cols = j;
		for (int curr_rows = 0;curr_rows < curr_kernal.rows;curr_rows++) {
			for (int curr_cols = 0;curr_cols < curr_kernal.cols;curr_cols++) {
				float a = curr_kernal.at<float>(curr_rows, curr_cols)*complete_image.at<Vec3b>(curr_rows + bios_rows, curr_cols + bios_cols)[chan];
				sum += a;
			}
		}
		dst.at<Vec3b>(i, j)[chan] = (int)sum;
	}
	else if (complete_image.channels() == 1) {
			float sum = 0;
			int bios_rows = i;
			int bios_cols = j;
			for (int curr_rows = 0;curr_rows < curr_kernal.rows;curr_rows++) {
				for (int curr_cols = 0;curr_cols < curr_kernal.cols;curr_cols++) {
					float a = curr_kernal.at<float>(curr_rows, curr_cols)*complete_image.at<uchar>(curr_rows + bios_rows, curr_cols + bios_cols);
					sum += a;
				}
			}
			dst.at<uchar>(i, j) = (int)sum;
		}
	else {
			cout << "the type of image is not suitable!" << endl;return;
		}
	}
//边框像素补全
void My_Convolution::complete_image_transform(const Mat &image, Mat &dst) {
	if (!kernal_loaded) {
		cout << "kernal is empty!" << endl;
		return;
	}
	//Mat的type（）成员函数生成CV_<位数>(S/U/F)C<通道数>
	//初始化一个补全图像的大小。
	dst = Mat::zeros(2 * bios_y + image.rows, 2 * bios_x + image.cols, image.type());
	Rect real_roi_of_image = Rect(bios_x, bios_y, image.cols, image.rows);
	Mat real_mat_of_image = dst(real_roi_of_image);
	image.copyTo(real_mat_of_image);
}

my_nedianfilter.h

1
2
3
4
5
6
7
8
9
10
11
12
13
#ifndef MY_MEDIANFILTER
#define MY_MEDIANFILTER
#include <vector>
#include <opencv2/opencv.hpp>
using namespace cv;
using namespace std;
void MedianFilter(Mat& src, Mat& dst, int win_size);
void Complete_Image_Transform(Mat &src, Mat &comp, int size);
void Calc_Median(Mat& comp, Mat& dst, int r, int c, int s);
#endif

my_nedianfilter.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
#include "my_medianfilter.h"
void MedianFilter(Mat& src, Mat& dst, int win_size) {
	Mat comp = src.clone();
	Complete_Image_Transform(src, comp, win_size);
	int rows = comp.rows, cols = comp.cols;
	int start = win_size / 2;
	if (comp.channels() > 1)
	{
		//彩色图片通道分离
		vector<Mat> channels_c, channels_d;
		split(comp, channels_c);
		split(comp, channels_d);
		//滤波
		for (int i = 0; i < 3; i++)
			Calc_Median(channels_c[i], channels_d[i], rows, cols, start);
		//合并彩色通道
		merge(channels_d, dst);
	}
	else
		Calc_Median(comp, dst, rows, cols, start);
}
//边框像素补全
void Complete_Image_Transform(Mat &src, Mat &comp, int size) {
	
	//初始化一个补全图像的大小
	comp = Mat::zeros(2 * (size / 2) + src.rows, 2 * (size / 2) + src.cols, src.type());
	Rect real_roi_of_image = Rect((size / 2), (size / 2), src.cols, src.rows);
	Mat real_mat_of_image = comp(real_roi_of_image);
	src.copyTo(real_mat_of_image);
}
void Calc_Median(Mat& comp, Mat& dst, int r, int c, int s)
{
	for (int m = s; m < r - s; m++) {
		for (int n = s; n < c - s; n++) {
			vector<uchar> model;
			for (int i = -s + m; i <= s + m; i++) {
				for (int j = -s + n; j <= s + n; j++) {
					//cout << int(src.at<uchar>(i, j)) << endl;
					model.push_back(comp.at<uchar>(i, j));
				}
			}
			sort(model.begin(), model.end());     //采用快速排序进行
			dst.at<uchar>(m - s, n - s) = model[(s * 2 + 1) * (s * 2 + 1) / 2];
		}
	}
}

C++默认初始化和值初始化、 直接初始化和拷贝初始化

默认初始化和值初始化一一对应

前提知识
声明： 在环境/上下文中指定一个变量的名字。也就是说，声明仅仅是让编译器知道，而没有实际分配空间。

初始化：给一个声明后尚未初始化的变量一个有意义的初始值。

赋值 ： 销毁一个变量原来的值，并赋予一个新值。相当于改变了一个变量的状态。

默认初始化：

对象可能产生未定义的值（是否绝对尚待确定，例如类类型（定义了默认构造函数）的默认初始化是否属于值初始化等）。

值初始化：

对象的值是确定（预设）的。

值初始化出现场景：

1.数组初始化时初始值数量小于维度，剩下的元素会进行值初始化；

1
{ int array[10] = {1, 2, 3}; }

2.当我们不使用初始值定义一个静态变量；

1
2
3
//未经初始化的全局静态变量和局部静态变量都会被自动初始化为0
static int n1;		//n1值初始化为0
{ static int n2; }	//n2值初始化为0

3.形如T()的表达式显示地请求值初始化；

1
2
3
4
{ 
	std::string *pia1 = new int[10]();	//动态分配10个值初始化为0的int
	std::string *pia2 = new int[10];	//动态分配10个未初始化的int
}

4.只提供vector对象容纳的元素数量而略去初始值，此时库会创建一个值初始化的元素初值(c++ primer p88)

1
vector<int> ivec(10);	//10个元素都被初始化为0

默认初始化出现场景：

1.块作用域内不使用任何初始值定义一个非静态变量；

1
{ int a; }//a的值是未定义的，但若在任何函数之外，会被初始化为0

2.类通过默认构造函数来控制默认初始化过程，默认构造函数以如下规则初始化类的数据成员

  • 如果存在类内初始值，用它来初始化成员
  • 否则，默认初始化该成员

直接初始化和拷贝初始化一一对应

直接初始化：

直接调用与实参匹配的构造函数，形式如“T t(u)”。

拷贝初始化：

拷贝初始化首先使用指定构造函数创建一个临时对象，然后用拷贝构造函数将那个临时对象复制到正在创建的对象”，形式如“T t=u”。

分布式和集中式的版本控制的区别




  1. 开始了解Git


  1. 1. 先谈谈版本控制的一些事

  2. 2. Git诞生背后的一些故事

  3. 3. 版本控制：集中式VS分布式

  4. 4. Git的思想和基本工作原理

  5. 5. Git在Windows下的安装




先说集中式版本控制系统，版本库是集中存放在中央服务器的，而干活的时候，用的都是自己的电脑，所以要先从中央服务器取得最新的版本，然后开始干活，干完活了，再把自己的活推送给中央服务器。中央服务器就好比是一个图书馆，你要改一本书，必须先从图书馆借出来，然后回到家自己改，改完了，再放回图书馆。





集中式版本控制系统最大的毛病就是必须联网才能工作，如果在局域网内还好，带宽够大，速度够快，可如果在互联网上，遇到网速慢的话，可能提交一个10M的文件就需要5分钟，这还不得把人给憋死啊。



那分布式版本控制系统与集中式版本控制系统有何不同呢？首先，分布式版本控制系统根本没有“中央服务器”，每个人的电脑上都是一个完整的版本库，这样，你工作的时候，就不需要联网了，因为版本库就在你自己的电脑上。既然每个人电脑上都有一个完整的版本库，那多个人如何协作呢？比方说你在自己电脑上改了文件A，你的同事也在他的电脑上改了文件A，这时，你们俩之间只需把各自的修改推送给对方，就可以互相看到对方的修改了。



和集中式版本控制系统相比，分布式版本控制系统的安全性要高很多，因为每个人电脑里都有完整的版本库，某一个人的电脑坏掉了不要紧，随便从其他人那里复制一个就可以了。而集中式版本控制系统的中央服务器要是出了问题，所有人都没法干活了。



在实际使用分布式版本控制系统的时候，其实很少在两人之间的电脑上推送版本库的修改，因为可能你们俩不在一个局域网内，两台电脑互相访问不了，也可能今天你的同事病了，他的电脑压根没有开机。因此，分布式版本控制系统通常也有一台充当“中央服务器”的电脑，但这个服务器的作用仅仅是用来方便“交换”大家的修改，没有它大家也一样干活，只是交换修改不方便而已。





当然，Git的优势不单是不必联网这么简单，后面我们还会看到Git极其强大的分支管理，把SVN等远远抛在了后面。



CVS作为最早的开源而且免费的集中式版本控制系统，直到现在还有不少人在用。由于CVS自身设计的问题，会造成提交文件不完整，版本库莫名其妙损坏的情况。同样是开源而且免费的SVN修正了CVS的一些稳定性问题，是目前用得最多的集中式版本库控制系统。



除了免费的外，还有收费的集中式版本控制系统，比如IBM的ClearCase（以前是Rational公司的，被IBM收购了），特点是安装比Windows还大，运行比蜗牛还慢，能用ClearCase的一般是世界500强，他们有个共同的特点是财大气粗，或者人傻钱多。



微软自己也有一个集中式版本控制系统叫VSS，集成在Visual Studio中。由于其反人类的设计，连微软自己都不好意思用了。



分布式版本控制系统除了Git以及促使Git诞生的BitKeeper外，还有类似Git的Mercurial和Bazaar等。这些分布式版本控制系统各有特点，但最快、最简单也最流行的依然是Git！