clear;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%读取待校正图像和标准图像
[FileName,PathName,FilterIndex] = uigetfile({'*.tif;*.bmp;*.jpg;*.png','All image files'},'Select the warp image');%获取待校正图像的信息(文件名、路径等)
[FileName1,PathName1,FilterIndex1] = uigetfile({'*.tif;*.bmp;*.jpg;*.png','All image files'},'Select the base image');%获取标准图像的信息(文件名、路径等)
if FilterIndex
pic_orig=imread([PathName,FileName]); %读取待校正图像
else
return
end
if FilterIndex1
base_orig=imread([PathName1,FileName1]); %读取标准图像
else
return
end
figure(1);
% subplot(2,4,1);
imshow(pic_orig);%显示待校正图像
title('待校正图像');
figure(2);
% subplot(2,4,5);
imshow(base_orig);%显示标准图像
title('标准图像');
%%%%%%%%%%%%%%%%%%%%%%%%%% 判断图像是彩色图像还是灰度图像,如果是彩色图像把它转换为灰度图像
[x,y,z] = size(pic_orig);
[x1,y1,z1] = size(base_orig);
if(z~=1)
pic_gray = rgb2gray(pic_orig);
end
if(z1~=1)
base_gray2 = rgb2gray(base_orig);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%对标准图像做灰度拉伸,增强其灰度对比
base_gray=imadjust(base_gray2,[0.4 0.7],[0 1]);
%%%%%%%%%%%%%%%%%%%%%%%%%%把灰度图像转换成二值图像
pic_bw=im2bw(pic_gray,0.9);%待校正图像可以直接二值化,效果较好,阈值自己调整
base_gray=double(base_gray); %标准图像先要进行低通滤波去除背景条纹,数据格式转变为double类型
base_gray_fft=(fftshift(fft2(base_gray)));%快速傅里叶变换
% base_gray_magfft=abs(base_gray_fft);
% base_gray_magfft=(base_gray_magfft-min(min(base_gray_magfft)))/(max(max(base_gray_magfft))-min(min(base_gray_magfft)))*255;%调整幅度
r = 120; %低通的中心半径,自己观察
[m,n]=size(base_gray_fft);
m1=fix(m/2);
n1=fix(n/2);
for i=1:m
for j=1:n
if (((i-m1)^2+(j-n1)^2) < r^2 )
mask(i,j) = 1;
else
mask(i,j) = 0;
end
end
end
base_gray2_fft=base_gray_fft.*mask;
base_gray2=ifft2(ifftshift(base_gray2_fft));
base_gray2=uint8(real(base_gray2));
base_bw=im2bw(base_gray2,0.9);%标准图像二值化,阈值自己调整
base_bw=~base_bw;%二值图像取反
figure(3);
% subplot(2,4,2);
imshow(pic_bw);
title('待校正图像二值化');
figure(4);
% subplot(2,4,6);
imshow(base_bw);
title('标准图像二值化');
[L,num] = bwlabel(base_bw,8);
figure(5);
imshow(uint8(L));
%%%%%%%%%%%%%%%%%%%%%%%%%%forstner提取特征点
%第一步:先用4个方向的差分算子提取初选点
result_pic=zeros(x,y); %特征提取结果
for i=2:x-1
for j=2:y-1
dg1_pic=abs(pic_bw(i,j)-pic_bw(i+1,j));
dg2_pic=abs(pic_bw(i,j)-pic_bw(i,j+1));
dg3_pic=abs(pic_bw(i,j)-pic_bw(i-1,j));
dg4_pic=abs(pic_bw(i,j)-pic_bw(i,j-1));
dg_pic=[dg1_pic,dg2_pic,dg3_pic,dg4_pic];
temp_pic=sort(dg_pic); %对4个差分算子进行升序排列
if temp_pic(3)==1 %有任意两个大于阈值,则该像素有可能是一特征点
result_pic(i,j)=255; %是初选点
else
result_pic(i,j)=0; %不是初选点
end
end
end
result_base=zeros(x1,y1);
for i=2:x1-1
for j=2:y1-1
dg1_base=abs(base_bw(i,j)-base_bw(i+1,j));
dg2_base=abs(base_bw(i,j)-base_bw(i,j+1));
dg3_base=abs(base_bw(i,j)-base_bw(i-1,j));
dg4_base=abs(base_bw(i,j)-base_bw(i,j-1));
dg_base=[dg1_base,dg2_base,dg3_base,dg4_base];
temp_base=sort(dg_base); %对4个差分算子进行升序排列
if temp_base(3)==1 %有任意两个大于阈值,则该像素有可能是一特征点
result_base(i,j)=255; %是初选点
else
result_base(i,j)=0; %不是初选点
end
end
end
figure(5);
% subplot(2,4,3);
imshow(result_pic);
title('待校正图像的初选点');
figure(6);
% subplot(2,4,7);
imshow(result_base);
title('标准图像的初选点');
%第二步:在以初选点为中心的3*3的窗口中计算协方差矩阵与圆度
wMatrix_pic=zeros(x,y); %权重矩阵
Tq_pic=0.9; %阈值,可设置
for i=2:x-1
for j=2:y-1
if result_pic(i,j)==255 %如果是初选点
gu2_pic=0.0;
gv2_pic=0.0;
guv_pic=0.0;
for ii=i-1:i
for jj=j-1:j
gu2_pic=gu2_pic+(pic_bw(ii+1,jj+1)-pic_bw(ii,jj))^2;
gv2_pic=gv2_pic+(pic_bw(ii,jj+1)-pic_bw(ii+1,jj))^2;
guv_pic=guv_pic+(pic_bw(ii+1,jj+1)-pic_bw(ii,jj))*(pic_bw(ii,jj+1)-pic_bw(ii+1,jj));
end
end
DetN_pic=gu2_pic*gv2_pic-guv_pic^2;
trN_pic=gu2_pic+gv2_pic;
q_pic=4*DetN_pic/(trN_pic*trN_pic);
%第三步:设定阈值Tq,若满足则计算权值
if q_pic > Tq_pic
wMatrix_pic(i,j)=DetN_pic/trN_pic;
else
result_pic(i,j)=0;
end
end
end
end
wMatrix_base=zeros(x1,y1); %下面是针对标准图像的第二步和第三步处理
Tq_base=0.8; %阈值
for i=2:x1-1
for j=2:y1-1
if result_base(i,j)==255
gu2_base=0.0;
gv2_base=0.0;
guv_base=0.0;
for ii=i-1:i
for jj=j-1:j
gu2_base=gu2_base+(base_bw(ii+1,jj+1)-base_bw(ii,jj))^2;
gv2_base=gv2_base+(base_bw(ii,jj+1)-base_bw(ii+1,jj))^2;
guv_base=guv_base+(base_bw(ii+1,jj+1)-base_bw(ii,jj))*(base_bw(ii,jj+1)-base_bw(ii+1,jj));
end
end
DetN_base=gu2_base*gv2_base-guv_base^2;
trN_base=gu2_base+gv2_base;
q_base=4*DetN_base/(trN_base*trN_base);
if q_base > Tq_base
wMatrix_base(i,j)=DetN_base/trN_base;
else
result_base(i,j)=0;
end
end
end
end
%第四步:以权值为基础,在一定窗口内抑制局部非最大值候选点,取出局部极大值点
vwsize=5; %选择5*5的窗口
wradius=floor(vwsize/2);
for i=wradius+1:x-wradius
for j=wradius+1:y-wradius
if result_pic(i,j)==255
tempiv_pic=wMatrix_pic(i-wradius:i+wradius,j-wradius:j+wradius);
tempsort_pic=sort(tempiv_pic(:),'descend');%将区域内像素按从大至小排列
if wMatrix_pic(i,j)==tempsort_pic(1) && wMatrix_pic(i,j)~=tempsort_pic(2)%排除整个区域权值相等的情况
;
else
result_pic(i,j)=0;
end
end
end
end
for i=wradius+1:x1-wradius
for j=wradius+1:y1-wradius
if result_base(i,j)==255
tempiv_base=wMatrix_base(i-wradius:i+wradius,j-wradius:j+wradius);
tempsort_base=sort(tempiv_base(:),'descend');%将区域内像素按从大至小排列
if wMatrix_base(i,j)==tempsort_base(1) && wMatrix_base(i,j)~=tempsort_base(2)%排除整个区域权值相等的情况
;
else
result_base(i,j)=0;
end
end
end
end
figure(7);
% subplot(2,4,4);
imshow(result_pic);
title('待校正图像的最终匹配点');
figure(8);
% subplot(2,4,8);
imshow(result_base);
title('标准图像的最终匹配点');
%%%%%%%%%%%%%%%%%%%%%%%%%%特征点匹配
cnt_pic=0;
cnt_base=0;
for i=1:x
for j=1:y
if result_pic(i,j)==255 %计算特征点个数
cnt_pic=cnt_pic+1;
end
end
end
for i=1:x1
for j=1:y1
if result_base(i,j)==255
cnt_base=cnt_base+1;
end
end
end
ssd=zeros(cnt_base,cnt_pic);%以邻域像素差的平方和最小为基准寻找同名点
ssd=zeros();
n=0;
for i=1:x
for j=1:y
if result_pic(i,j)==255
PIC=pic_bw(i-2:i+2,j-2:j+2);
n=n+1
m=0;
for ii=1:x1
for jj=1:y1
if result_base(ii,jj)==255
BASE=base_bw(i-2:i+2,j-2:j+2);
m=m+1
ssd(m,n)=sum(sum((PIC-BASE).*(PIC-BASE)));
[ssd_max,max_index]=max(ssd,[],1);
end
end
end
end
end
end