电子科技大学图像与视频处理（图像部分）课后作业

                                    电子科技大学图像与视频处理（图像部分）课后作业

作业2-1

1.已知一幅图像有8个灰度级，各个灰度级出现的概率如下表所示，画出该图像直方图，然后对该图像的直方图进行均衡化处理，并画出均衡化后的灰度直方图。

2. 编程练习：用Matlab绘制出一幅M×N 的灰度图像的灰度直方图；对该图像进行均衡化处理，绘制出均衡后图像的直方图。

提示：
(1) 读图；
(2) 调用灰度直方图显示函数；
(3) 调用灰度直方图均衡化函数。
有兴趣的同学尝试不调用后两个函数进行编程实现

I=imread('pout.tif');       %读取原图像
I=im2double(I);
 
%% 显示原图像
subplot(2,2,1);
imshow(I);
subplot(2,2,2);
imhist(I);
 
%% 图像均衡化
subplot(2,2,3);
imshow(histeq(I));
subplot(2,2,4);
imhist(histeq(I));

作业2-2

试分别计算I与A和B两个模板的滤波结果，输出结果要求与I的维数相同，需要考虑边界处理问题（边界填充0值）。并说明A与B两个模板对处理结果有什么不同响应？

作业2-3

用Harris算子和归一化互相关实现下面两幅图的匹配。原图见群文件。

主程序：

I1=imread('C:\Users\matalab2018a\ch03\boat1.bmp'); % 注意修改自己的路径
I2=imread('C:\Users\matalab2018a\ch03\boat3.bmp');
I1=rgb2gray(I1);
I2=rgb2gray(I2);
coords1=Harris(I1);
coords2=Harris(I2);
%===================================
% 计算互相关矩阵R
%===================================
[~,m]=size(coords1); %获取图像1的角度个数
[~,n]=size(coords2); %获取图像2的角度个数
R=zeros(m,n);      %互相关矩阵
[r1,c1]=size(I1);
[r2,c2]=size(I2);
r=2;               %角点匹配领域半径
I1_tmp=zeros(r1+2*r,c1+2*r);
I2_tmp=zeros(r2+2*r,c2+2*r);
I1_tmp(r+1:r1+r,r+1:c1+r)=I1;
I2_tmp(r+1:r2+r,r+1:c2+r)=I2;
for i=1:m
x1=coords1(1,i)+r;
y1=coords1(2,i)+r;
for j=1:n
x2=coords2(1,j)+r;
y2=coords2(2,j)+r;
I1_mean=mean(mean(I1_tmp(x1-r:x1+r,y1-r:y1+r)));
I2_mean=mean(mean(I2_tmp(x2-r:x2+r,y2-r:y2+r)));
a=I1_tmp(x1-r:x1+r,y1-r:y1+r)-I1_mean;
b=I2_tmp(x2-r:x2+r,y2-r:y2+r)-I2_mean;
p=sum(sum(a.*b));
q=sqrt(sum(sum(a.^2))*sum(sum(b.^2)));
R(i,j)=p/q;
end
end
%===================================
% 寻找相关性最大的角点对并绘图
%===================================
% 构造画布
canvas=ones(max(r1,r2),c1+c2+100)*255;
canvas(1:r1,1:c1)=I1;
canvas(1:r2,c1+101:c1+c2+100)=I2;
imshow(uint8(canvas));
hold on
 
[~,match_points]=max(R,[],2);       % 选取相关系数最大的匹配点
d=10;
for i=1:m
    x1=coords1(2,i);
    x2=coords2(2,match_points(i));
    y1=coords1(1,i);
    y2=coords2(1,match_points(i));
    plot([x1 x2+c1+100],[y1 y2],'ro');% 绘制角点
    % 绘制匹配角点连线
    % 去除相关系数小、在图像边缘的角点对
    if R(i)>0.32 &&...
        x1>d && x1< c1-d && y1>d && y2< r1-d &&...
        x2>d && x2<c2-d && y2>d && y2<r2-d
        line([x1,x2+c1+100],[y1,y2]);
        % plot([x1 x2+c1+100],[y1 y2],'ro');
    end
end

子程序

function coords=Harris(I)
[m,n]=size(I);
%===================================
% Step 1:计算相关矩阵M
%===================================
tmp=zeros(m+2,n+2);
tmp(2:m+1,2:n+1)=I;     % 补零
Ix=zeros(m+2,n+2);      % x方向梯度
Iy=zeros(m+2,n+2);      % y方向梯度
Ix(:,2:n)=tmp(:,3:n+1)-tmp(:,1:n-1);
Iy(2:m,:)=tmp(3:m+1,:)-tmp(1:m-1,:);
Ix2=Ix(2:m+1,2:n+1).^2;
Iy2=Iy(2:m+1,2:n+1).^2;
Ixy=Ix(2:m+1,2:n+1).*Iy(2:m+1,2:n+1);
h=fspecial('gaussian',[7 7],2);
Ix2=filter2(h,Ix2);
Iy2=filter2(h,Iy2);
Ixy=filter2(h,Ixy);
%===================================
% Step 2:计算Harris角点响应
%===================================
k=0.05;               % k常取0.04-0.06
Rmax=0;
R=zeros(m,n);
for i=1:m
for j=1:n
M=[Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];
R(i,j)=det(M)-k*(trace(M))^2;
if R(i,j)>Rmax
Rmax=R(i,j);
end
end
end
%===================================
% Step 3:寻找Harris角点
%===================================
tmp(2:m+1,2:n+1)=R;
%res=zeros(m+2,n+2);
%res(2:m+1,2:n+1)=I;
coords=[];
for i=2:m+1
for j=2:n+1
% 在3*3邻域内非极大值抑制
if tmp(i,j)>0.01*Rmax &&...
tmp(i,j)>tmp(i-1,j-1) && tmp(i,j)>tmp(i-1,j) && tmp(i,j)>tmp(i-1,j+1) &&...
tmp(i,j)>tmp(i,j-1) && tmp(i,j)>tmp(i,j+1) &&...
tmp(i,j)>tmp(i+1,j-1) && tmp(i,j)>tmp(i+1,j) && tmp(i,j)>tmp(i+1,j+1)
coords=[coords,[i-1;j-1]];
end
end
end
end

效果图