这是2020年春季学期我选的外专业课“数字图像处理”的课程设计，觉得比较有趣，故拿来记录留念。
在这门课程中，我第一次对图像处理进行实践，老师讲得非常好，推荐大家选这门课，能学到不少知识。
课设题目
天文摄影师常常用照相机对夜晚的天空进行拍摄。v1.avi和v2.avi是两段持续几秒的夜晚天空景象的视频。可以注意到，夜晚较暗的光照使得视频中有明显的噪声。设计一个算法用于去除这两组图像数据中的噪声。
根据学过的图像处理知识可知，为了实现通过多帧图像数据获取一幅去除噪声后的图像，一种常用的方法是平均去噪法，即计算多帧图像的平均值。
分别计算这两组多帧图像数据的平均值，获得两幅去噪图像。
用上述方法可以发现去噪结果并不令人满意，原因是各帧图像间存在一定的偏移。请设计一个算法克服这种偏移同时获得较好的平均法去噪结果。
实验原理
以视频中的第一帧为基准，逐帧遍历视频图片，由于图片位移间隔时间小，位移小，可以通过特殊方法获得当前帧图片与第一帧基准图片的相对运动大小，求出变换矩阵，通过变换的逆矩阵反向平移当前帧，即可认为观察者未运动，再使用平均值法去噪即可。
由于运动时间短，我们忽略运动的旋转部分，仅考虑平移运动。平移运动为二自由度，因此我们只需要找到两张图中一个对应的特征点，即可求得变换矩阵。
涉及的 MATLAB 函数
VideoReader;
用于读取视频中的内容，将其每一帧都读出来。
im2double；
将读入的图像转换为浮点形式，便于后期处理。
rgb2gray；
将rgb图像转化为灰度图，便于后期处理。
strel；
创建形态学操作的构造元素。
imerode；
图像形态学腐蚀处理。
edge
边缘检测。
find、max；
求得数组最大值并确定最大值的位置。
实验内容与方法
首先观察两个视频，可以看出，视频中图像含有许多噪声：
这些噪声分布均匀且均值基本为0，可以通过均值滤波的方法消除噪声。但是，由于视频拍摄的抖动导致直接进行多帧均值滤波会使运动物体产生模糊，导致滤波效果变差：
因此，想到了使用特征点的方法求出前后两帧的相对运动关系，将后一帧移动到前一帧对应的位置，以消除运动产生的影响。
首先使用im2double与rgb2gray将图像处理为浮点型的单值图，之后对两段视频分别选择寻找特征点的方法。
对于第一个视频，使用缩小范围的寻找特征点的方法，在原图中找到一个sector区域，作为特征点运动的范围。选取这个位置要保证范围内有且仅有一个比较明亮的点，并且在运动过程中这个点不会跑出这个范围。
sector = frame(421:472,681:728);
处理前的sector
对其进行处理，首先进行形态学的腐蚀处理
形态学操作可以改变物体形状
形态学腐蚀就是求局部最小值的操作
se = strel(‘disk’,1);
sector = imerode(sector,se);
之后使用find函数找到图片中幅值最大的像素，以这个点的位置作为特征点。
[x_1,y_1] = find(sector == max(max(sector)));
处理后的sector
选择这个点作为特征点进行唯一判断。
对于第二个视频，我们使用图中月亮的圆心作为特征点。
首先对图像进行canny边缘检测：
result_frame = edge(frame,’canny’,[0.3,0.5]);
Canny边缘检测包含以下四个步骤：
1.高斯滤波
滤波的主要目的是降噪，一般的图像处理算法都需要先进行降噪。而高斯滤波主要使图像变得平滑（模糊），同时也有可能增大了边缘的宽度。
2.计算图片梯度
边缘一定就是图像中像素变化较快的区域，因此使用梯度就可以检测出边缘。我们使用sobel算子对图片进行滤波，就可以计算出它的梯度。
3.过滤非最大值
在高斯滤波过程中，边缘有可能被放大了。这个步骤使用一个规则来过滤不是边缘的点，使边缘的宽度尽可能为1个像素点：如果一个像素点属于边缘，那么这个像素点在梯度方向上的梯度值是最大的。否则不是边缘，将灰度值设为0。这样就可以使检测到的边缘始终是一个像素。
4.使用上下阈值确定边缘
通过给定的上下阈值确定最终哪些地方属于边缘，哪些地方不是。梯度过小的地方不应该是边缘，而梯度过大的地方也不应被视为边缘 。
想要检测图片中的圆首先想到，可以使用Hough圆变换来求解圆的圆心：
首先理解霍夫线变换。霍夫线变换的基本理论，就是二进制图片中每一个点都可能属于某些线。如果将过(x0,y0)的这些条线参数化为一个斜率为a，截距为b的线，由于这条线过这一点，可以将a，b作为变量，画出这个点对应的所有直线参数，在a，b平面中的形式：
也是一条直线，并且斜率与截距都是确定的。对于每一个图像中的非零点，都进行这样的操作，并且在a，b平面表示出来。当所有的非零点都表示出来，并进行累加后，a，b平面上数值比较大的点就是原图中可能的直线的参数。
同理，霍夫圆变换与直线变换的原理相同，先将圆参数化，只不过这时需要三个参数，得到的参数平面将会是一个三位体中的坐标。这样也是可行的，只不过运算缓慢，通常使用霍夫梯度法来检测圆。圆心是圆所有法线的交汇点，通过梯度，确定这些法线，并求出他们的交点中最有可能的那一个，就能检测到圆心。随后再判断圆心到圆周上点的距离，相同或相近距离数量多，就能确定圆的半径，并且也能检验圆心选择的正确性。
在思考过程中，又发现进行完canny边缘检测后，这张图像十分完美，故使用运算更为简便的方法进行圆检测。
使用扫描算法检测这个圆的圆心，具体方法为：
逐列遍历整个图像，若在每一列检测到两次白点，计算两个白点之间的距离，所有距离中最长的距离就是圆的直径，这时两个点中间的的那个点就是圆的圆心。
实验结果
实验所使用代码（Matlab）如下：
对于第一个视频：
clear,clc;
vidobj = VideoReader(['v1' '.avi']);
numFrames = vidobj.NumFrames;
sum_frame = zeros(533,800,3);
origin_frame = im2double(read(vidobj,1));
frame = rgb2gray(origin_frame);
sector = frame(421:472,681:728);
se = strel('disk',1);
sector = imerode(sector,se);
[x_1,y_1] = find(sector == max(max(sector)));
sum_frame = origin_frame;
for i = 2:numFrames
origin_frame = im2double(read(vidobj,i));
frame = rgb2gray(origin_frame);
%find offset
sector = frame(421:472,681:728);
se = strel('disk',1);
sector = imerode(sector,se);
[x,y] = find(sector == max(max(sector)));
x_offset = x - x_1;
y_offset = y - y_1;
if x_offset > 0 && y_offset > 0
origin_frame(1:533-x_offset,1:800-y_offset,:) = origin_frame(x_offset+1:533,y_offset+1:800,:);
elseif x_offset < 0 && y_offset > 0
origin_frame(-x_offset+1:533,1:800-y_offset,:) = origin_frame(1:533+x_offset,y_offset+1:800,:);
elseif x_offset > 0 && y_offset < 0
origin_frame(1:533-x_offset,-y_offset+1:800,:) = origin_frame(x_offset+1:533,1:800+y_offset,:);
elseif x_offset < 0 && y_offset < 0
origin_frame(-x_offset+1:533,-y_offset+1:800,:) = origin_frame(1:533+x_offset,1:800+y_offset,:);
end
sum_frame = sum_frame + origin_frame;
end
sum_frame = sum_frame/60;
imshow(sum_frame);对于第二个视频：
clear,clc;
vidobj = VideoReader(['v2' '.avi']);
numFrames = vidobj.NumFrames;
sum_frame = zeros(600,800,3);
origin_frame = im2double(read(vidobj,1));
frame = rgb2gray(origin_frame);
result_frame = edge(frame,'canny',[0.3,0.5]);
%find x,y
for j = 1:600
n = 0;
for k = 1:800
if result_frame(j,k) == 1
n = n+1;
if n == 1
startpoint = k;
else
if k - startpoint > 10
maxVal(j,1) = k - startpoint;
maxVal(j,2) = startpoint;
else
n = n-1;
end
end
end
if k == 800 && (n == 1 || n == 0)
maxVal(j,1) = 0;
maxVal(j,2) = 0;
end
end
end
r = max(maxVal(:,1));
xpos = find(maxVal(:,1)==r);
xpos = int16(median(xpos));
ypos = maxVal(xpos,2) + int16(r/2);
baseX = xpos;
baseY = ypos;
sum_frame = origin_frame;
for i = 2:numFrames
origin_frame = im2double(read(vidobj,i));
frame = rgb2gray(origin_frame);
result_frame = edge(frame,'canny',[0.3,0.5]);
%find x,y
for j = 1:600
n = 0;
for k = 1:800
if result_frame(j,k) == 1
n = n+1;
if n == 1
startpoint = k;
else
if k - startpoint > 10
maxVal(j,1) = k - startpoint;
maxVal(j,2) = startpoint;
else
n = n-1;
end
end
end
if k == 800 && (n == 1 || n == 0)
maxVal(j,1) = 0;
maxVal(j,2) = 0;
end
end
end
r = max(maxVal(:,1));
xpos = find(maxVal(:,1)==r);
xpos = int16(median(xpos));
ypos = maxVal(xpos,2) + int16(r/2);
x_offset = xpos - baseX;
y_offset = ypos - baseY;
if x_offset > 0 && y_offset > 0
origin_frame(1:600-x_offset,1:800-y_offset,:) = origin_frame(x_offset+1:600,y_offset+1:800,:);
elseif x_offset < 0 && y_offset > 0
origin_frame(-x_offset+1:600,1:800-y_offset,:) = origin_frame(1:600+x_offset,y_offset+1:800,:);
elseif x_offset > 0 && y_offset < 0
origin_frame(1:600-x_offset,-y_offset+1:800,:) = origin_frame(x_offset+1:600,1:800+y_offset,:);
elseif x_offset < 0 && y_offset < 0
origin_frame(-x_offset+1:600,-y_offset+1:800,:) = origin_frame(1:600+x_offset,1:800+y_offset,:);
end
sum_frame = sum_frame + origin_frame;
end
sum_frame = sum_frame/60;
imshow(sum_frame);运行结果：
可以看到，经过处理的图像噪声明显减少，并且并未受到观察者视角移动而出现动态模糊的影响。