空域图像增强 | Euan's Blog

空域图像增强

Euan 9月 17, 2019

目的

掌握OPENCV中空域图像增强的一些基本方法
理解图像基于空域像素的处理：图像像素灰度级、点运算、逻辑运算、几何运算
理解图像直方图均衡化空域增强的基本原理
理解图像的空间滤波增强的基本原理和方法：平滑空间滤波、锐化空间滤波处理；比较不同滤波器处理的效果和优缺点

原理

图像的点运算、代数运算、几何运算、直方图增强、空域平滑滤波、空域锐化滤波

实验内容或步骤:

对“lena.jpg”变换后的灰度图像进行如下操作:

1）统计每个灰度级（0-255）的像素点数，并打印显示；
2）实现幂律变换,并显示；（,将其封装为一个函数，公式中的参数(c, )作为函数参数传入）
3）用双线性插值将图像放到1.5倍和缩小0.8倍，并显示
4）将图像进行直方图均衡化处理，显示图像及其直方图；
5）分别以3x3、7x7模板进行均值滤波(或中值滤波)处理，并显示；
6）分别使用拉普拉斯算子和Sobel算子对图像进行锐化处理，并显示。

问题讨论:

图像像素访问一般有三种方式，选择一种使用即可
参考
https://blog.csdn.net/u013921430/article/details/81114644
改变图像大小（resize）
void resize(InputArray src, OutputArray dst, Size dsize, double fx=0, double fy=0, int interpolation=INTER_LINEAR );
src：输入，原图像，即待改变大小的图像；
dst：输出，改变大小之后的图像；
dsize：输出图像的大小。
fx和fy是图像width方向和height方向的缩放比例
注意：dsize和fx/fy不能同时起作用，当dsize不为空的时候，fx，fy被忽略
例如：resize(img, imgDst, Size(100,100));将图像缩放到指定的像素大小
resize(img, imgDst, Size(), 0.5, 0.5); 将图像缩小一倍
interpolation：指定插值的方式
INTER_NEAREST - 最邻近插值
NTER_LINEAR - 双线性插值（默认参数）
INTER_CUBIC - 4x4像素邻域内的双立方插值
INTER_LANCZOS4 - 8x8像素邻域内的Lanczos插值
直方图均衡（equalizeHist）
https://docs.opencv.org/3.4.5/d4/d1b/tutorial_histogram_equalization.html
https://docs.opencv.org/3.4.5/d8/dbc/tutorial_histogram_calculation.html
均值滤波（blur）中值滤波(medianBlur)
https://docs.opencv.org/3.4.5/dc/dd3/tutorial_gausian_median_blur_bilateral_filter.html
拉普拉斯算子（Laplacian）
https://docs.opencv.org/3.4.5/d4/d86/group__imgproc__filter.html#gad78703e4c8fe703d479c1860d76429e6
http://www.opencv.org.cn/opencvdoc/2.3.2/html/modules/imgproc/doc/filtering.html#laplacian
Sobel算子（Sobel）
https://docs.opencv.org/3.0.0/d4/d86/group__imgproc__filter.html#gacea54f142e81b6758cb6f375ce782c8d
http://www.opencv.org.cn/opencvdoc/2.3.2/html/modules/imgproc/doc/filtering.html#sobel

实验报告要求：
（1）描述实验的基本原理和步骤。
（2）用数据和图片给出各个步骤中取得的实验结果，包括原始图像及其处理后的图像，并进行必要的讨论。
（3）完整的源代码。

步骤

统计像素点数，打印显示

实验方式1 傻瓜式操作

灰度化原理：R=G=B

灰度值0-255，则灰度级为256
思路：使用指针访问记录各点像素并且输入到txt中

由于该图像的尺寸为500*500，则有250000个像素点

再通过 python（比较熟悉）进行统计每一集的像素数目

如图：一个键为灰度级，值为该灰度值的数目，并按键值从小到大排序

实现方式2 正确操作

指针访问记录各点像素，通过遍历每个像素值，并通过c++内嵌统计数值得出，之前做得太傻了，有轮子不用自己造……

反转变换并显示

原理: 图像颜色的反转,比较简单的思路就是使用255减去当前值,从而得到反转后的图像

作用: 强图像的暗区中白色或灰色的细节，特别是黑色面积在尺寸上占主导地位时

幂律变换并显示

公式：s=cr^γ

其中c、γ 为常数。考虑偏移量上式可写为 s=c(ε+r)^γ

伽马变换可以很好地拉伸图像的对比度，扩展灰度级。

和s取值范围[0,1]，使用归一化处理参数

使用幂律变换进行对比度增强；灰度级压缩

函数封装：

显示

双线性插值

原理：假设源图像大小为mxn，目标图像为axb。那么两幅图像的边长比分别为：m/a和n/b。注意，通常这个比例不是整数，编程存储的时候要用浮点型。目标图像的第（i,j）个像素点（i行j列）可以通过边长比对应回源图像。其对应坐标为（im/a,jn/b）

显然，这个对应坐标一般来说不是整数，而非整数的坐标是无法在图像这种离散数据上使用的。双线性插值通过寻找距离这个对应坐标最近的四个像素点，来计算该点的值（灰度值或者RGB值）
双线性插值加速以及优化策略

这个过程很折腾，一开始还没报错，后来疯狂内存问题报错，后来花了2个小时解决了……后来重新做时，不小心改了cv的源代码……后来脑子转过来改回来正常运行

resize可以直接跳用函数进行放缩，通过原理的算法经过多次缩放比例实验后，发现一旦有整数的话很容易成功，如果放大1.5杯，算的过程中例如fx等参数最后时int类型，通过cout输入int原本是float类型可能类型强制转换导致最后像素值的值冲突

直方图均衡化

原理：直方图变换其实是一种灰度变换，灰度变换的变换函数决定了输入随机变量与输出随机变量之间的关系，也就是两个随机变量的关系；一副图像是二维离散的数据，不利于使用数学的工具进行处理，在数字图像处理中，我们通常是采用连续的变量进行推导，最后在推广到离散的情况。

作用：通过直方图均衡技术，将图像的灰度分布变得较为均匀，从而使得图像对比度增大，视觉效果更佳

遍历全图，先统计每个灰度级下的像素点个数（为此我们开辟了256大小的数组）；
计算每个灰度级的像素点占总像素的点的比例；
按照第二步求出的比例重新计算每个灰度级下的新的灰度值，即均衡化；
依照新的灰度值表遍历更新图像的灰度值。

直方图绘制

直方图calcHist函数参数介绍

void calcHist(
const Mat* images, //源(图像)数组。应具有相同的深度(CV_8U或CV_32F)和大小。每个数组可以有任意数量的通道。
int nimages, //源图像的数量
const int* channels, //用于计算直方图的dims通道列表
InputArray mask, //可选掩码。如果矩阵不为空，它必须是与images[i]大小相同的8位数组。
OutputArray hist, //输出直方图
int dims, //直方图的维度(必须是正数且不能大于CV_MAX_DIMS)
const int* histSize, //每个维度中的直方图大小的数组。
const float** ranges, //每个维度中直方图边界边界的dims数组的数组。
bool uniform = true, //指示直方图是否均匀的标志
bool accumulate = false //累积标志。若设置为true，分配的直方图不会被清除，用于及时更新直方图。
);

均值滤波及中值滤波

均值滤波：线性平均滤波器，它通过求窗口内所有像素的平均值来得到中心像素点的像素值。这样的好处是可以有效的平滑图像，降低图像的尖锐程度，降低噪声。但缺点是不能消除噪声，而且将周围的景物的像素点平均了，变得模糊了。

中值滤波：如上课所说，椒盐噪声有很好的抑制作用，对模板里的像素来个简单的冒泡排序，然后求个中值代替中心点像素就ok了

图片中一个方块区域内，中心点的像素为全部点像素值的平均值。

缺陷：均值滤波本身存在着固有的缺陷，即它不能很好地保护图像细节，在图像去噪的同时也破坏了图像的细节部分，从而使图像变得模糊，不能很好地去除噪声点。特别是椒盐噪声
在图像中取3*3的矩阵，里面有9个像素点，我们将9个像素进行排序，最后将这个矩阵的中心点赋值为这九个像素的中值。

模糊化太会过于严重。

Laplace&Sobel

一开始看到这两个算法以为是边缘检测，后来感觉不对重新看了下题目发现是锐化。

边缘检测:字面意思就是在局部范围内与周围环境明显不同

Laplace步骤

高斯模糊-去噪声
转为灰度图像
拉普拉斯-二阶导数计算
取绝对值convertScaleAbs()
显示结果

Sobel步骤

高斯平滑
转灰度
求梯度xy
振幅图像（合并近似梯度）

锐化

Laplace原理：根据图像某个像素的周围像素到此像素的突变程度有关，也就是说它的依据是图像像素的变化程度。我们知道，一个函数的一阶微分描述了函数图像是朝哪里变化的，即增长或者降低；而二阶微分描述的则是图像变化的速度，急剧增长下降还是平缓的增长下降。那么据此我们可以猜测出依据二阶微分能够找到图像的色素的过渡程度，例如白色到黑色的过渡就是比较急剧的。
步骤:得到图像矩阵后，用二阶分然后再和愿图像相加

Sobel锐化也可以通过边缘检测然后和原图像相加得出
可以看出在边缘检测方面，sobel 产生的边缘抗噪性好，导致很多原图的像素被去除，而Laplace 对边缘敏感，可能有些是噪声的边缘，也被算进来了

源代码

#include<iostream>
#include<opencv2/core/core.hpp>
#include<opencv2/highgui/highgui.hpp>
#include<fstream>
#include <opencv2/opencv.hpp>

using namespace cv;
using namespace std;
int i;
int j;

// 归一化
//data 进行处理的像素集合
//grayscale 目标灰度级
//rows cols type 目标图像的行，列，以及类型

int gray[256] = { 0 };  //记录每个灰度级别下的像素个数
double gray_prob[256] = { 0 };  //记录灰度分布密度
double gray_distribution[256] = { 0 }; //记录累计密度
int gray_equal[256] = { 0 }; //均衡化后的灰度值

int gray_sum = 0; //像素总数

Mat Normalize(vector<double> data, int grayscale, int rows, int cols, int type)
{
double max = 0.0;
double min = 0.0;
for (int i = 0; i < data.size(); i++)
{
if (data[i] > max)
max = data[i];
if (data[i] < min)
min = data[i];
}
Mat dst;
dst.create(rows, cols, type);
int index = 0;
for (int r = 0; r < dst.rows; r++)
{
uchar* dstRowData = dst.ptr<uchar>(r);
for (int c = 0; c < dst.cols; c++)
{
dstRowData[c] = (uchar)(grayscale * ((data[index++] - min) * 1.0 / (max - min)));
}
}
return dst;
}

// 幂律变化
Mat PowerTranseform(Mat src, double gamma, int parameter)
{
Mat dst;
dst.create(src.size(), src.type());
vector<double> value;
for (int r = 0; r < src.rows; r++)
{
uchar* srcRowData = src.ptr<uchar>(r);
for (int c = 0; c < src.cols; c++)
{
//幂次变换的公式为 s = c * r ^ v r为输入图像像素
value.push_back(parameter * pow((double)srcRowData[c], gamma));
}
}
return Normalize(value, 255, src.rows, src.cols, src.type());
}


// 双线性插值
Mat bilinear(Mat src, int row, int col) {
int rows = src.rows, cols = src.cols;
cv::Mat dst(row, col, src.type());
for (int i = 0; i < row; ++i) {
//以ptr的方式访问dst的数据
uchar* p = dst.ptr<uchar>(i);
//使两个图像的几何中心重合，采样更合理
float x = (i + 0.5) * rows / row - 0.5;
int fx = (int)x;
//x为坐标的小数部分
x -= fx;
//以整数计算速度更快
short x1 = (1.f - x) * 2048;
short x2 = 2048 - x1;
for (int j = 0; j < col; ++j) {
//trick
float y = (j + 0.5) * cols / col - 0.5;
int fy = (int)y;
y -= fy;
//trick
short y1 = (1.f - y) * 2048;
short y2 = 2048 - y1;
//结果右移22位抵消2048的平方
p[j] = (src.at<uchar>(fx, fy) * x1 * y1 + src.at<uchar>(fx + 1, fy) * x2 * y1
+ src.at<uchar>(fx, fy + 1) * x1 * y2 + src.at<uchar>(fx + 1, fy + 1) * x2 * y2) >> 22;
}
}
return dst;
}

Mat equalize_hist(Mat& input)
{
Mat output = input.clone();
gray_sum = input.cols * input.rows;

//统计每个灰度下的像素个数
for (int i = 0; i < input.rows; i++)
{
uchar* p = input.ptr<uchar>(i);
for (int j = 0; j < input.cols; j++)
{
int vaule = p[j];
//cout << input.cols*i+j << ':'<< vaule << endl;
gray[vaule]++;
}
}
//统计灰度频率
for (int i = 0; i < 256; i++)
{
gray_prob[i] = ((double)gray[i] / gray_sum);
cout << i + 1 << ':' << gray_prob[i] * gray_sum << endl;
}
return output;
}

//均值化
void histBar(Mat src)
{
int histSize = 30;
MatND hist;
float range[] = { 0, 255 };
const float* ranges = { range };
//计算直方图
calcHist(&src, 1, 0, Mat(), hist, 1, &histSize, &ranges, true, false);
//将直方图bin的数值归一化到0-255
normalize(hist, hist, 0, 255, NORM_MINMAX, -1, Mat());

//显示直方图
int binsW = cvRound((double)500 / histSize);
Mat histImg = Mat::zeros(500, 500, CV_8UC3);
RNG rng(123);
for (int i = 0; i < histSize; i++)
{
Scalar color = Scalar(rng.uniform(0, 255), rng.uniform(0, 255), rng.uniform(0, 255));
rectangle(histImg, Point(i * binsW, 500), Point((i + 1) * binsW, 500 - cvRound(hist.at<float>(i) * 500 / 255.0)), color, -1);
}
imshow("histogram", histImg);
}

void sharpen(const Mat image, Mat& result)
{
result.create(image.size(), image.type());
for (int j = 1; j < image.rows - 1; j++)
{
const uchar* previous = image.ptr<const uchar>(j - 1);
const uchar* current = image.ptr<const uchar>(j);
const uchar* next = image.ptr<const uchar>(j + 1);

uchar* output = result.ptr<uchar>(j);
for (int i = 1; i < image.cols - 1; i++)
{
//sharpened_pixel=5*current-left-right-up-down;
//cv::saturate_cast用以对计算结果进行截断(0-255)
*output++ = saturate_cast<uchar>(
5 * current[i] - current[i - 1]
- current[i + 1] - previous[i] - next[i]);
}
}
result.row(0).setTo(Scalar(0));
result.row(result.rows - 1).setTo(Scalar(0));
result.col(0).setTo(Scalar(0));
result.col(result.cols - 1).setTo(Scalar(0));
}

int main()
{
Mat img = cv::imread("lena_gray.jpg", 0);
int row = img.rows;//获取行数
int col = img.cols * img.channels();//注意是列数*通道数
cout << row << "," << col << endl;
equalize_hist(img);

/*ofstream ofs("pixel.txt");
for (int i = 0; i < row; ++i)
{
uchar* data = img.ptr<uchar>(i);//获取第i行首地址
for (int j = 0; j < col; ++j)
{
//指针方式
ofs << int(*data) << endl;
data++;
}
}
cout << row << "," << col << endl;
cout << "遍历完毕" << endl;
cout << endl;
ofs.close();*/

//显示原图
imshow("img", img);
Mat dst;
bitwise_not(img, dst, noArray());

//show image 
imshow("dst", dst);

//幂律变换
Mat dstImg;
dstImg = PowerTranseform(img, 0.4, 1);
imshow("变换后", dstImg);

// 双线性插值
Mat scale1, scale2;
/*scale1 = bilinear(img, row * 1.5, col * 1.5);
imshow("1.5", scale1);
cv::imwrite("1.5.jpg", scale1);*/
Mat dst_img1,dst_img2;
// 双线性插值(resize和源算法)
resize(img, dst_img1, Size(), 1.5, 1.5);
imshow("resize_1.5", dst_img1);
imwrite("resize_1.5.jpg", dst_img1);
resize(img, dst_img2, Size(), 0.8, 0.8);
imshow("resize_0.8", dst_img2);
imwrite("resize_0.8.jpg", dst_img2);
scale2 = bilinear(img, row * 0.8, col * 0.8);
imshow("0.8", scale2);
imwrite("0.8.jpg", scale2);

histBar(img);

//进行均值滤波操作
Mat out;
blur(img, out, Size(3, 3));
imshow("均值滤波3*3", out);
blur(img, out, Size(5, 5));
imshow("均值滤波5*5", out);
blur(img, out, Size(7, 7));
imshow("均值滤波7*7", out);

//进行均值滤波操作
Mat g_mDstImage;
medianBlur(img, g_mDstImage, 3);
imshow("中值滤波3*3", out);
medianBlur(img, g_mDstImage, 5);
imshow("中值滤波5*5", out);
medianBlur(img, g_mDstImage, 7);
imshow("中值滤波7*7", out);

//边缘检测
//sobel算子
Mat gao_src, gray_src, xdst, ydst, dst_sobel;
/*GaussianBlur(img, gao_src, Size(3, 3), 0, 0);//高斯滤波
/imshow("gaosioutput", gao_src);
cvtColor(gao_src, gray_src, CV_RGB2GRAY);
//imshow("grayoutput", gray_src);//转换为灰度图*/
Sobel(img, xdst, -1, 1, 0);
imshow("xdst", xdst);
Sobel(img, ydst, -1, 0, 1);
imshow("ydst", ydst);
addWeighted(xdst, 0.5, ydst, 0.5, 1, dst_sobel);
imshow("sobel", dst_sobel);
//拉普拉斯算子
Mat gaosrc, graysrc, ladst, dst_la;
/*GaussianBlur(img, gaosrc, Size(3, 3), 0);
//imshow("gaosrc", gaosrc);
cvtColor(gaosrc, graysrc, CV_RGB2GRAY);
//imshow("graysrc", graysrc);*/
Laplacian(img, ladst, -1);
//imshow("ladst", ladst);
convertScaleAbs(ladst, dst_la);
imshow("la", dst_la);

//锐化
//拉普拉斯算子
Mat la_dst, sobel_dst;
la_dst.create(img.size(), img.type());
sharpen(img, la_dst);
imshow("sharpen_la", la_dst);
//sobel算子
//Generate Gradient X
//imshow("sharpen_la", sobel_dst);

waitKey(0); // 无限制的等待用户的按键事件
return 0;
}

本文使用 CC BY-NC-SA 3.0 中国大陆协议许可
具体请参见知识共享协议
本文链接：https://zyhang8.github.io/2019/09/17/cv-exp2/