图像分割配准与融合
程序员文章站
2024-01-16 09:08:40
...
内容：
- 图像目标边界描述
- 图像分割
- 图像配准
- 图像融合
GithubPage:https://zhangquan1995.github.io
图像目标边界描述
图像边界链码
给定一个二值图像，需要快速获取二值图像中连通体边界线的链码。
function out=chaincode4(image)
%功能：实现4连通链码
%输入： 二值图像
%输出：链码的结果
%4连通的扫描方向
n=[0 1;-1 0;0 -1;1 0];
%设置标志
flag=1;
%初始输出的链码串为空
cc=[];
%找到起始点
[x y]=find(image==1);
x=min(x);
imx=image(x,:);
y=min(find(imx==1));
first=[x y];
dir=3;
while flag==1
tt=zeros(1,4);
newdir=mod(dir+3,4);
for i=0:3
j=mod(newdir+i,4)+1;
tt(i+1)=image(x+n(j,1),y+n(j,2));
end
d=min(find(tt==1));
dir=mod(newdir+d-1,4);
%找到下一个像素点的方向码后补充在链码的后面
cc=[cc,dir];
x=x+n(dir+1,1);y=y+n(dir+1,2);
%判别链码的结束标志
if x==first(1)&&y==first(2)
flag=0;
end
end
out=cc;
function out=chaincode8(image)
%功能：实现8连通链码
%输入： 二值图像
%输出：链码的结果
n=[0 1;-1 1;-1 0;-1 -1;0 -1;1 -1;1 0;1 1];
%设置标志
flag=1;
%初始输出的链码串为空
cc=[];
%找到起始点
[x y]=find(image==1);
x=min(x);
imx=image(x,:);
y=min(find(imx==1));
first=[x y];
dir=7;
while flag==1
tt=zeros(1,8);
newdir=mod(dir+7-mod(dir,2),8);
for i=0:7
j=mod(newdir+i,8)+1;
tt(i+1)=image(x+n(j,1),y+n(j,2));
end
d=min(find(tt==1));
dir=mod(newdir+d-1,8);
%找到下一个像素点的方向码后补充在链码的后面
cc=[cc,dir];
x=x+n(dir+1,1);y=y+n(dir+1,2);
%判别链码的结束标志
if x==first(1)&&y==first(2)
flag=0;
end
end
out=cc;clear all;clc;
L=zeros(7,6);L(2:6,2:3)=1;L(5:6,4:5)=1;
c4=chaincode4(L) %4连通边界链码
c8=chaincode8(L) %8连通边界链码 目标平移，链码不会变化；目标旋转，链码发生变化。可对链码进行旋转归一化。
图像分割
在对图像的研究和应用中，人们往往仅仅对图像中的某些部分感兴趣，这些部分被称为前景或目标，其余部分则称为背景。目标一般对应于图像中特定的、具有独特性质的区域。独特性质可以是像素的灰度值、物体轮廓曲线、颜色和纹理等。为了识别和分析图像中的目标，需要将它们从图像中分离提取出来，在此基础上才有可能进一步对目标进行测量和对图像进行利用。
基于阈值的图像分割
阈值分割法是图像分割中最具代表性的一种分割算法。由于图像阈值处理的直观性和易于实现的特点，以及阈值分割总能用封闭而连通的边界定义不交叠的区域，使得阈值化分割算法成为图像分割中最为常见的分割法之一。它特别适合于目标和背景占据不同灰度级范围的图像。分为全局阈值和局部阈值。
全局阈值法
- 基于灰度直方图的阈值分割
%读入图像
A=imread('hehua.jpg');
B=rgb2gray(A);
B=double(B);
%求图像的灰度直方图
hist(B)
[m,n]=size(B);
%根据直方图进行阈值分割
for i=1:m
for j=1:n
%阈值
if B(i,j)>70&B(i,j)<130
B(i,j)=1;
else
B(i,j)=0;
end
end
end
%显示分割结果
subplot(121),imshow(A)
subplot(122),imshow(B)- 迭代阈值分割
%读入图像,并进行灰度转换
A=imread('baihe.jpg');
B=rgb2gray(A);
%初始化阈值
T=0.5*(double(min(B(:)))+double(max(B(:))));
d=false;
%通过迭代求最佳阈值
while~d
g=B>=T;
Tn=0.5*(mean(B(g))+mean(B(~g)));
d=abs(T-Tn)<0.5;
T=Tn;
end
% 根据最佳阈值进行图像分割
level=Tn/255;
BW=im2bw(B,level);
% 显示分割结果
subplot(121),imshow(A)
subplot(122),imshow(BW)- 最大类间方差阈值分割法
I=imread('moli.jpg');
%rgb转灰度
if isrgb(I)==1
I_gray=rgb2gray(I);
else
I_gray=I;
end
subplot(121),imshow(I_gray);
%转化为双精度
I_double=double(I_gray);
[wid,len]=size(I_gray);
%灰度级
colorlevel=256;
%直方图
hist=zeros(colorlevel,1);
%计算直方图
for i=1:wid
for j=1:len
m=I_gray(i,j)+1;
hist(m)=hist(m)+1;
end
end
%直方图归一化
hist=hist/(wid*len);
miuT=0;
for m=1:colorlevel
miuT=miuT+(m-1)*hist(m);
end
xigmaB2=0;
for mindex=1:colorlevel
threshold=mindex-1;
omega1=0;
omega2=0;
for m=1:threshold-1
omega1=omega1+hist(m);
end
omega2=1-omega1;
miu1=0;
miu2=0;
for m=1:colorlevel
if m<threshold
miu1=miu1+(m-1)*hist(m);
else
miu2=miu2+(m-1)*hist(m);
end
end
miu1=miu1/omega1;
miu2=miu2/omega2;
xigmaB21=omega1*(miu1-miuT)^2+omega2*(miu2-miuT)^2;
xigma(mindex)=xigmaB21;
if xigmaB21>xigmaB2
finalT=threshold;
xigmaB2=xigmaB21;
end
end
%阈值归一化
fT=finalT/255
for i=1:wid
for j=1:len
if I_double(i,j)>finalT
bin(i,j)=1;
else
bin(i,j)=0;
end
end
end
subplot(122),imshow(bin);MATLAB提供了全局阈值函数graythresh()，采用最大类间方差阈值计算图像的全局阈值，与im2bw()结合可将灰度图转化为二值图像。
I=imread('shuixian.jpg');
Level=graythresh(I);
BW=im2bw(I,Level);
subplot(121),imshow(I);
subplot(122),imshow(BW);局部阈值法
当图像中有阴影、光照不均匀、各处的对比度不同、有突发噪声、背景灰度变化等，如果只用一个固定的全局阈值对整幅图像进行分割，则由于不能兼顾图像各处的情况而使分割效果受到影响。有一种解决办法就是用与像素值位置相关的一组阈值来对图像各部分分别进行分割。这种与坐标相关的阈值称为动态阈值、局部阈值或自适应阈值。优点是抗噪声能力强，对一些用全局阈值不能分割的图像有较好的效果，缺点是时间复杂度和空间复杂度比较大。
clear all;close all;clc;
im1 = double(imread('tshape.png'));
bwim1 = adaptivethreshold(im1,15,0.02,0);
subplot(1,2,1);
imshow(uint8(im1));
subplot(1,2,2);
imshow(bwim1);function bw=adaptivethreshold(IM,ws,C,tm)
% 功能：自适应图像分割
% 输入：
% IM-待分割的原始图像
% ws-平均滤波时的窗口大小，可参考fspecial的用法
% C-常量，需要根据经验选取合适的参数
% tm-开关变量，tm=1进行中值滤波，tm=0则进行均值滤波
% bw-图像分割后输出的二值图像
% 输入参数处理
if (nargin<3)
error('You must provide the image IM, the window size ws, and C.');
elseif (nargin==3)
tm=0;
elseif (tm~=0 && tm~=1)
error('tm must be 0 or 1.');
end
IM=mat2gray(IM);
if tm==0
% 图像均值滤波
mIM=imfilter(IM,fspecial('average',ws),'replicate');
else
% 图像中值滤波
mIM=medfilt2(IM,[ws ws]);
end
sIM=mIM-IM-C;
bw=im2bw(sIM,0);
bw=imcomplement(bw);基于区域生长法的图像分割
区域生长法是根据同一物体区域内像素的相似性质来聚集像素点的方法，从初始区域开始，将相邻的具有同样性质的像素或其他区域并到目前的区域中，从而逐步增长区域，直至没有可以归并的点或其他小区域为止。区域内像素的相似性度量可以包括平均灰度值、纹理、颜色等信息。
- 选取待分割的区域内的一点作为种子点
- 以种子点为中心，考虑其四个领域像素，如果满足生长规则，则并入种子点，将此像素压入堆栈
- 从堆栈中选取一个像素，当做种子点，回到上一步
- 当堆栈为空时，生长结束
clear all;close all;clc;
I = im2double(imread('medtest.png'));
x=198;
y=359;
J = regiongrowing(I,x,y,0.2);
figure(1), imshow(I);
figure(2), imshow(J);function J=regiongrowing(I,x,y,reg_maxdist)
% 功能：基于区域生长法的图像分割
% 输入：I-待分割的输入图像 x,y-种子点的坐标
% t –最大密度距离(默认值为0.2)
% 输出：J-分割后的图像
% 示例:
% I = im2double(imread('medtest.png'));
% x=198; y=359;
% J = regiongrowing(I,x,y,0.2);
% figure, imshow(I+J);
if(exist('reg_maxdist','var')==0), reg_maxdist=0.2; end
if(exist('y','var')==0), figure, imshow(I,[]); [y,x]=getpts; y=round(y(1)); x=round(x(1)); end
J = zeros(size(I));
Isizes = size(I);
% 分割区域的均值
reg_mean = I(x,y);
% 区域中的像素数
reg_size = 1;
neg_free = 10000; neg_pos=0;
neg_list = zeros(neg_free,3);
% 区域中新的像素距区域的距离
pixdist=0;
neigb=[-1 0; 1 0; 0 -1;0 1];
% 基于区域生长的图像分割
while(pixdist<reg_maxdist&®_size<numel(I))
% 添加新的邻域像素
for j=1:4,
% 计算邻域坐标
xn = x +neigb(j,1); yn = y +neigb(j,2);
% 检查邻域是否在图像内部
ins=(xn>=1)&&(yn>=1)&&(xn<=Isizes(1))&&(yn<=Isizes(2));
if(ins&&(J(xn,yn)==0))
neg_pos = neg_pos+1;
neg_list(neg_pos,:) = [xn yn I(xn,yn)]; J(xn,yn)=1;
end
end
if(neg_pos+10>neg_free), neg_free=neg_free+10000; neg_list((neg_pos+1):neg_free,:)=0; end
dist = abs(neg_list(1:neg_pos,3)-reg_mean);
[pixdist, index] = min(dist);
J(x,y)=2; reg_size=reg_size+1;
% 计算区域新的均值
reg_mean= (reg_mean*reg_size + neg_list(index,3))/(reg_size+1);
x = neg_list(index,1); y = neg_list(index,2);
neg_list(index,:)=neg_list(neg_pos,:); neg_pos=neg_pos-1;
end
J=J>1;基于最大方差法灰度门限的图像分割
分割门限选择的准确性直接影响分割精度及图像描述分析的正确性。
clear all;close all;clc;
a=imread('shuixian.jpg');
a=rgb2gray(a);
th=thresh_md(a)%计算灰度阈值function th=thresh_md(a)
%功能：实现最大方差法计算分割门限
%输入：a-为灰度图像
%输出: th-灰度分割门限
%返回图像矩阵a各个灰度等级像素个数
count=imhist(a);
[m,n]=size(a);
N=m*n-sum(sum(find(a==0),1));
%指定图像灰度等级为256级
L=256;
%计算出各灰度出现的概率
count=count/N;
%找出出现概率不为0的最小灰度
for i=2:L
if count(i)~=0
st=i-1;
break;
end
end
%找出出现概率不为0的最大灰度
for i=L:-1:1
if count(i)~=0
nd=i-1;
break;
end
end
f=count(st+1:nd+1);
%p和q分别为灰度起始和结束值
p=st;
q=nd-st;
%计算图像的平均灰度
u=0;
for i=1:q
u=u+f(i)*(p+i-1);
ua(i)=u;
end
%计算出选择不同k值时，A区域的概率
for i=1:q
w(i)=sum(f(1:i));
end
%求出不同k值时类间的方差
d=(u*w-ua).^2./(w.*(1-w));
%求出最大方差对应的灰度级
[y,tp]=max(d);
th=tp+p; 基于k-means算法的图像分割
k-means算法，也被称为k-平均或k-均值，是一种得到广泛使用的聚类算法。它是将各个聚类子集内的数据样本的均值作为该聚类的代表点，算法的主要思想是通过迭代过程把数据集划分为不同的类别，使得评价聚类性能的准则函数达到最优，从而使生成的每个聚类紧凑、类间独立。
- 选定某种距离作为数据样本间的相似性度量
- 选择评价聚类性能的准则函数
- 相似度的计算根据一个簇中对象的平均值来进行
算法步骤：
- 为每个聚类确定一个初始聚类中心，即k个聚类中心
- 将样本集中的样本按照最小距离原则分配到最邻近聚类
- 使用每个聚类中的样本均值作为新的聚类中心
- 重复以上两步骤，直到聚类中心不再变换
图像配准
图像配准可定义成两相邻图像之间的空间变换和灰度变换，即先将图像像素的坐标映射到一个新坐标系中的某一坐标，再对其像素进行重采样。图像配准要求相邻图像之间有一部分在逻辑上是相同的，即相邻的图像有一部分反映了同一目标区域，这一点是实现图像配准的基本条件。如果确定了相邻图像代表同一场景的所有像素之间的坐标关系，采用相应的处理算法即可实现图像配准。一般地，空间变换和灰度变换是非线性变换。
相似性测度：均方根误差、差绝对值和误差、马氏距离和相似度。
基于灰度的图像配准
clear all;close all;clc;
I=imread('lena.jpg');
I=rgb2gray(I);
T=I(50:100,75:125);
% 输入： T-模板 I-输入的原始图像
% I_SSD-采用像素差平方和法（SSD）的匹配结果
% I_NCC-采用标准化互相关匹配法（NCC）的匹配结果
[I_SSD,I_NCC]=template_matching(T,I);
subplot(1,2,1);imshow(I_SSD);
subplot(1,2,2);imshow(I_NCC);function [I_SSD,I_NCC]=template_matching(T,I)
% 功能：图像配准
% [I_SSD,I_NCC]=template_matching(T,I)
% 输入： T-模板 I-输入的原始图像
% I_SSD-采用像素差平方和法（SSD）的匹配结果
% I_NCC-采用标准化互相关匹配法（NCC）的匹配结果
% 将图像转换成双精度型
T=double(T); I=double(I);
if(size(T,3)==3)
% 如果是彩色图像，则按彩色图像匹配方法进行匹配
[I_SSD,I_NCC]=template_matching_color(T,I);
else
% 如果是灰度图像，则按灰度图像匹配方法进行匹配
[I_SSD,I_NCC]=template_matching_gray(T,I);
end
function [I_SSD,I_NCC]=template_matching_color(T,I)
% 功能：对彩色图像进行匹配子函数，其核心原理是从R、G、B三个子色调进行匹配
[I_SSD_R,I_NCC_R]=template_matching_gray(T(:,:,1),I(:,:,1));
[I_SSD_G,I_NCC_G]=template_matching_gray(T(:,:,2),I(:,:,2));
[I_SSD_B,I_NCC_B]=template_matching_gray(T(:,:,3),I(:,:,3));
% 融合三次匹配结果
I_SSD=(I_SSD_R+I_SSD_G+I_SSD_B)/3;
I_NCC=(I_NCC_R+I_NCC_G+I_NCC_B)/3;
function [I_SSD,I_NCC]=template_matching_gray(T,I)
% 功能：对灰度图像进行匹配子函数
T_size = size(T); I_size = size(I);
outsize = I_size + T_size-1;
% 在频域内进行相关运算
if(length(T_size)==2)
FT = fft2(rot90(T,2),outsize(1),outsize(2));
FI = fft2(I,outsize(1),outsize(2));
Icorr = real(ifft2(FI.* FT));
else
FT = fftn(rot90_3D(T),outsize);
FI = fftn(I,outsize);
Icorr = real(ifftn(FI.* FT));
end
LocalQSumI= local_sum(I.*I,T_size);
QSumT = sum(T(:).^2);
% 计算模板和图像的像素差平方和
I_SSD=LocalQSumI+QSumT-2*Icorr;
% 将其皈依化到0和1之间
I_SSD=I_SSD-min(I_SSD(:));
I_SSD=1-(I_SSD./max(I_SSD(:)));
I_SSD=unpadarray(I_SSD,size(I));
if (nargout>1)
LocalSumI= local_sum(I,T_size);
stdI=sqrt(max(LocalQSumI-(LocalSumI.^2)/numel(T),0) );
stdT=sqrt(numel(T)-1)*std(T(:));
meanIT=LocalSumI*sum(T(:))/numel(T);
I_NCC= 0.5+(Icorr-meanIT)./ (2*stdT*max(stdI,stdT/1e5));
I_NCC=unpadarray(I_NCC,size(I));
end
function T=rot90_3D(T)
T=flipdim(flipdim(flipdim(T,1),2),3);
function B=unpadarray(A,Bsize)
Bstart=ceil((size(A)-Bsize)/2)+1;
Bend=Bstart+Bsize-1;
if(ndims(A)==2)
B=A(Bstart(1):Bend(1),Bstart(2):Bend(2));
elseif(ndims(A)==3)
B=A(Bstart(1):Bend(1),Bstart(2):Bend(2),Bstart(3):Bend(3));
end
function local_sum_I= local_sum(I,T_size)
B = padarray(I,T_size);
if(length(T_size)==2)
s = cumsum(B,1);
c = s(1+T_size(1):end-1,:)-s(1:end-T_size(1)-1,:);
s = cumsum(c,2);
local_sum_I= s(:,1+T_size(2):end-1)-s(:,1:end-T_size(2)-1);
else
s = cumsum(B,1);
c = s(1+T_size(1):end-1,:,:)-s(1:end-T_size(1)-1,:,:);
s = cumsum(c,2);
c = s(:,1+T_size(2):end-1,:)-s(:,1:end-T_size(2)-1,:);
s = cumsum(c,3);
local_sum_I = s(:,:,1+T_size(3):end-1)-s(:,:,1:end-T_size(3)-1);
end用以上相关算法进行模板匹配的计算量很大，其中，模板要在（N-M+1）^2个参考位置上做相关计算。
序贯相似性检测算法（SSDA）
- 定义绝对误差值
- 取一个不变阈值
- 在子图中随机选取像点，计算它同T中对应点的误差值，然后把这个差值和其他点对的差值累加起来，当累加r次误差超过Tk,则停止累加，并记下次数，作为检测曲面。
- 把取值最大的检测曲面对应的点作为匹配点。
基于特征点的图像配准
基于Harris角点的图像配准
- 采用Harris角点检测算法分别检测两幅输入图像的Harris角点
- 角点领域八个像素值作为匹配特征点向量，则第一幅图像的特征匹配点向量为N1 X 8维，第二幅图像特征匹配点为N2 X 8维。
- 将特征点向量进行匹配，取相似度小于阈值的作为两幅图像的配准点。
- 需要剔除匹配错误点
function test()
%直接运行
%读入第一幅图像并进行Harris角点检测;
img11 = imread('scence1.jpg');
img1 = rgb2gray(img11);
img1 = double(img1(:,:));
pt1 = kp_harris(img1);
%读入第二幅图像并进行Harris角点检测;
img21 = imread('scence2.jpg');
img2 = rgb2gray(img21);
img2 = double(img2(:,:));
pt2 = kp_harris(img2);
% 进行匹配
result = match(img1,pt1,img2,pt2);
result(1,intersect(find(result(1,:) > 0),find(result(2,:) == 0))) = 0;
while(length(find(result(1,:)>0)) > 3)
result
draw2(img11,img21,pt1,pt2,result);
pause;
[index index] = max(result(2,:));
result(1,index(1)) = 0;
result(2,index(1)) = 0;
end
draw2(img11,img21,pt1,pt2,result);
endfunction draw2(img1,img2,pt1,pt2,result)
% 功能：功能显示匹配特征点子程序
h = figure;
subplot(1,2,1);
imshow(img1);
subplot(1,2,2);
imshow(img2);
s = size(pt1,2);
subplot(1,2,1);
for i=1:size(pt1,1)
rectangle('Position',[pt1(i,2)-s,pt1(i,1)-s,2*s,2*s],'Curvature',[00],'EdgeColor','b','LineWidth',2);
end
subplot(1,2,2);
for i=1:size(pt2,1)
rectangle('Position',[pt2(i,2)-s,pt2(i,1)-s,2*s,2*s],'Curvature',[00],'EdgeColor','b','LineWidth',2);
end
for i=1:size(result,2)
if(result(1,i)~=0)
subplot(1,2,1);
text(pt1(result(1,i),2)+3,pt1(result(1,i),1)+3,num2str(i),'BackgroundColor',[1 1 1]);
subplot(1,2,2);
text(pt2(i,2)+3,pt2(i,1)+3,num2str(i),'BackgroundColor',[1 1 1]);
end
end
endfunction result = match(img1,pt1,img2,pt2)
%功能：进行匹配子程序
regionValue1 = getRegionValue(img1,pt1);
len1 = size(regionValue1,2);
regionValue2 = getRegionValue(img2,pt2);
len2 = size(regionValue2,2);
%找出最佳匹配点
result = zeros(2,len2);
for i=1:len1
B = regionValue1(:,i);
[value,index] = sort(sum(abs(regionValue2-B(:,ones(1,size(regionValue2,2))))));
if((result(1,index(1))==0)||(result(2,index(1))>value(1)))
result(1,index(1))=i;
result(2,index(1))=value(1);
end
end
endfunction regionValue = getRegionValue(img,pt)
% 功能：分别取角点附近的8个像素值，作为匹配特征点向量
len = size(pt,1);
regionValue = zeros(8,len);
maxX = size(img,1);
maxY = size(img,2);
for i=1:len
x = pt(i,1);
y = pt(i,2);
%1
if(x-1<1||y-1<1)
regionValue(1,i)=100;
else
regionValue(1,i)=img(x,y)-img(x-1,y-1);
end
%2
if(x-1<1)
regionValue(2,i)=200;
else
regionValue(2,i)=img(x,y)-img(x-1,y);
end
%3
if(x-1<1||y+1>maxY)
regionValue(3,i)=300;
else
regionValue(3,i)=img(x,y)-img(x-1,y+1);
end
%4
if(y+1>maxY)
regionValue(4,i)=400;
else
regionValue(4,i)=img(x,y)-img(x,y+1);
end
%5
if(x+1>maxX||y+1>maxY)
regionValue(5,i)=500;
else
regionValue(5,i)=img(x,y)-img(x+1,y+1);
end
%6
if(x+1>maxX)
regionValue(6,i)=600;
else
regionValue(6,i)=img(x,y)-img(x+1,y);
end
%7
if(x+1>maxX||y-1<1)
regionValue(7,i)=700;
else
regionValue(7,i)=img(x,y)-img(x+1,y-1);
end
%8
if(y-1<1)
regionValue(8,i)=800;
else
regionValue(8,i)=img(x,y)-img(x,y-1);
end
end
endfunction points = kp_harris(im)
% 功能：检测图像的Harris角点
im = double(im(:,:,1));
sigma = 1.5;
s_D = 0.7*sigma;
x = -round(3*s_D):round(3*s_D);
dx = x .* exp(-x.*x/(2*s_D*s_D)) ./ (s_D*s_D*s_D*sqrt(2*pi));
dy = dx';
Ix = conv2(im, dx, 'same');
Iy = conv2(im, dy, 'same');
s_I = sigma;
g = fspecial('gaussian',max(1,fix(6*s_I+1)), s_I);
Ix2 = conv2(Ix.^2, g, 'same');
Iy2 = conv2(Iy.^2, g, 'same');
Ixy = conv2(Ix.*Iy, g, 'same');
cim = (Ix2.*Iy2 - Ixy.^2)./(Ix2 + Iy2 + eps);
[r,c,max_local] = findLocalMaximum(cim,3*s_I);
t = 0.6*max(max_local(:));
[r,c] = find(max_local>=t);
points = [r,c];
end
function [row,col,max_local] = findLocalMaximum(val,radius)
mask = fspecial('disk',radius)>0;
nb = sum(mask(:));
highest = ordfilt2(val, nb, mask);
second_highest = ordfilt2(val, nb-1, mask);
index = highest==val & highest~=second_highest;
max_local = zeros(size(val));
max_local(index) = val(index);
[row,col] = find(index==1);
end剔除匹配错误的点
clear all;close all;clc;
im1=imread('scence1.jpg');
p1=kp_harris(im1);
im1=rgb2gray(im1);
im2=imread('scence2.jpg');
p2=kp_harris(im2);
im2=rgb2gray(im2);
w=5;
p1=p1';
p2=p2';
[m1,m2,cormat]=matchbycorrelation(im1,p1,im2,p2,w);
% m1、m2为匹配点坐标，cormat为相关矩阵function[m1,m2,cormat] = matchbycorrelation(im1, p1, im2, p2, w, dmax)
% 功能:基于特征点的图像配准
% 输入：im1-灰度图像im1(需要转换成double型) im2-灰度图像im1(需要转换成double型)
% p1-图像im1上的特征点（2×N矩阵） p2-图像im2上的特征点（2×N矩阵）
% w-窗口大小 damx-最大匹配半径
% 输出：m1-灰度图像im1中匹配点的坐标 m2-灰度图像im1中匹配点的坐标
% cormat-相关矩阵
if nargin == 5
dmax = Inf;
end
im1 = double(im1);
im2 = double(im2);
im1 = im1 - filter2(fspecial('average',w),im1);
im2 = im2 - filter2(fspecial('average',w),im2);
% 产生相关矩阵
cormat = correlatiomatrix(im1, p1, im2, p2, w, dmax);
[corrows,corcols] = size(cormat);
% 查找相关矩阵中的最大值
[mp2forp1, colp2forp1] = max(cormat,[],2);
[mp1forp2, rowp1forp2] = max(cormat,[],1);
p1ind = zeros(1,length(p1));
p2ind = zeros(1,length(p2));
indcount = 0;
for n = 1:corrows
if rowp1forp2(colp2forp1(n)) == n
indcount = indcount + 1;
p1ind(indcount) = n;
p2ind(indcount) = colp2forp1(n);
end
end
% 整理匹配点坐标
p1ind = p1ind(1:indcount);
p2ind = p2ind(1:indcount);
% 从原始矩阵中提取匹配点坐标
m1 = p1(:,p1ind);
m2 = p2(:,p2ind); function cormat = correlatiomatrix(im1, p1, im2, p2, w, dmax)
% 功能:计算相关矩阵
% 输入：im1-灰度图像im1(需要转换成double型) im2-灰度图像im1(需要转换成double型)
% p1-图像im1上的特征点（2×N矩阵） p2-图像im2上的特征点（2×N矩阵）
% w-窗口大小 damx-最大匹配半径
% 输出：cormat-相关矩阵
if mod(w, 2) == 0
error('Window size should be odd');
end
[rows1, npts1] = size(p1);
[rows2, npts2] = size(p2);
cormat = -ones(npts1,npts2)*Inf;
if rows1 ~= 2 | rows2 ~= 2
error('Feature points must be specified in 2xN arrays');
end
[im1rows, im1cols] = size(im1);
[im2rows, im2cols] = size(im2);
% 相关运算窗口半径
r = (w-1)/2;
n1ind = find(p1(1,:)>r & p1(1,:)<im1rows+1-r & p1(2,:)>r & p1(2,:)<im1cols+1-r);
n2ind = find(p2(1,:)>r & p2(1,:)<im2rows+1-r & p2(2,:)>r & p2(2,:)<im2cols+1-r);
for n1 = n1ind
w1 = im1(p1(1,n1)-r:p1(1,n1)+r, p1(2,n1)-r:p1(2,n1)+r);
w1 = w1./sqrt(sum(sum(w1.*w1)));
if dmax == inf
n2indmod = n2ind;
else
p1pad = repmat(p1(:,n1),1,length(n2ind));
dists2 = sum((p1pad-p2(:,n2ind)).^2);
n2indmod = n2ind(find(dists2 < dmax^2));
end
for n2 = n2indmod
w2 = im2(p2(1,n2)-r:p2(1,n2)+r, p2(2,n2)-r:p2(2,n2)+r);
cormat(n1,n2) = sum(sum(w1.*w2))/sqrt(sum(sum(w2.*w2)));
end
end尝试通过小波提升算法建立感兴趣区域，提高图像灰度匹配的速度。有图像处理、目标粗定位、目标精确定位和模板更新组成。
图像融合
- 基于小波变换的图像融合
- 基于拉普拉斯金字塔变换的图像融合
图像融合就是将不同来源的同一对象的图像数据进行空间配准，然后采用一定的算法将各图像数据中所含的信息优势互补，产生新图像数据的信息。
融合三个层次：
- 像素级融合：低级，精度高，但实时性差
- 特征级融合：信息得到压缩，实时性好，精度差
基于空域的图像融合
- 图像像素灰度极值的融合法
- 图像像素灰度值加权融合法
- TOET图像融合方法
基于变换域的图像融合
- 基于多分辨率金字塔融合法
- 基于傅里叶变换的融合法
- 基于小波变换的融合法
%导入待融合图像1
load bust
X1=X;
map1=map;
subplot(131);image(X1);
colormap(map1);title('原始图像1');
axis square
%导入待融合图像2
load mask
X2=X;
map2=map;
%对灰度值大于100的像素进行增强，小于100的像素进行减弱
for i=1:256
for j=1:256
if(X2(i,j)>100)
X2(i,j)=1.2*X2(i,j);
else
X2(i,j)=0.5*X2(i,j);
end
end
end
subplot(132)
image(X2);colormap(map2);title('原始图像2');
axis square
%对原始图像1进行小波分解
[c1,s1]=wavedec2(X1,2,'sym4');
%对分解后的低频部分进行增强
sizec1=size(c1);
for I=1:sizec1(2)
c1(I)=1.2*c1(I);
end
%对原始图像2进行分解
[c2,s2]=wavedec2(X2,2,'sym4');
%将分解后的低频分量和高频分量进行相加，并乘以权重系数0.5
c=c1+c2;
c=0.5*c;
s=s1+s2;
s=0.5*s;
%进行小波重构
xx=waverec2(c,s,'sym4');
subplot(133);image(xx);title('融合图像');
axis square基于小波变换的图像融合
clear all;close all;clc;
I1=imread('clock1.bmp');
I2=imread('clock2.bmp');
dim=2;
% I1，I2 为两幅原始图像
% 对输入图像进行小波分解
y1=mywavedec2(I1,dim);
y2=mywavedec2(I2,dim);
% 根据低频融合算法进行图像融合
[r,c]=size(y1);
% 首先取两幅源图像相应的小波分解系数绝对值最大者的值作为融合图像的分解系数
for i=1:r
for j=1:c
if ( abs(y1(i,j)) >= abs(y2(i,j)) )
y3(i,j)=y1(i,j);
elseif ( abs(y1(i,j)) < abs(y2(i,j)) )
y3(i,j)=y2(i,j);
end
end
end
figure();
subplot(1,2,1);imshow(I1);
subplot(1,2,2);imshow(I2);
LLa = y1(1:r/(2^dim),1:c/(2^dim));
LLb = y2(1:r/(2^dim),1:c/(2^dim));
%调用lowfrefus函数对低频部分的小波分解系数进行融合
y4(1:r/(2^dim),1:c/(2^dim)) = lowfrefus(LLa,LLb);
figure();imshow(uint8(y4));%融合后的图像function y=mywavedec2(x,dim)
% 函数 MYWAVEDEC2() 对输入矩阵 x 进行 dim 层分解，得到相应的分解系数矩阵 y
% 输入参数：x —— 输入矩阵；
% dim —— 分解层数。
% 输出参数：y —— 分解系数矩阵。
x=modmat(x,dim); % 首先规范化输入矩阵，使其行列数均能被 2^dim 整除
subplot(121);imshow(x);title('原始图像'); % 画出规范化后的源图像
[m,n]=size(x); % 求出规范化矩阵x的行列数
xd=double(x); % 将矩阵x的数据格式转换为适合数值处理的double格式
for i=1:dim
xd=modmat(xd,1);
[dLL,dHL,dLH,dHH]=mydwt2(xd); % 矩阵小波分解
tmp=[dLL,dHL;dLH,dHH]; % 将分解系数存入缓存矩阵
xd=dLL; % 将缓存矩阵左上角部分的子矩阵作为下一层分解的源矩阵
[row,col]=size(tmp); % 求出缓存矩阵的行列数
y(1:row,1:col)=tmp; % 将缓存矩阵存入输出矩阵的相应行列
end
yd=uint8(y); % 将输出矩阵的数据格式转换为适合显示图像的uint8格式
for i=1:dim % 对矩阵 yd 进行分界线处理，画出分解图像的分界线
m=m-mod(m,2);
n=n-mod(n,2);
yd(m/2,1:n)=255;
yd(1:m,n/2)=255;
m=m/2;n=n/2;
end
subplot(122);imshow(yd);title([ num2str(dim) ' 层小波分解图像']);
%% 子函数 %%
function y=modmat(x,dim)
% 函数 MODMAT() 对输入矩阵x进行规范化，使其行列数均能被 2^dim 整除
% 输入参数：x —— r*c 维矩阵；
% dim —— 矩阵重构的维数
% 输出参数：y —— rt*ct 维矩阵，mod(rt,2^dim)=0，mod(ct,2^dim)=0
[row,col]=size(x); % 求出输入矩阵的行列数row,col
rt=row - mod(row,2^dim); % 将row,col分别减去本身模 2^dim 得到的数
ct=col - mod(col,2^dim); % 所得的差为rt、ct，均能被 2^dim 整除
y=x(1:rt,1:ct); % 输出矩阵 y 为输入矩阵 x 的 rt*ct 维子矩阵
%% 子函数 %%
function [cA,cD] = mydwt(x,lpd,hpd,dim);
% 函数 [cA,cD]=MYDWT(X,LPD,HPD,DIM) 对输入序列x进行一维离散小波分解，输出分解序列[cA,cD]
% 输入参数：x——输入序列；
% lpd——低通滤波器；
% hpd——高通滤波器；
% dim——小波分解级数。
% 输出参数：cA——平均部分的小波分解系数；
% cD——细节部分的小波分解系数。
cA=x; % 初始化cA，cD
cD=[];
for i=1:dim
cvl=conv(cA,lpd); % 低通滤波，为了提高运行速度，调用MATLAB提供的卷积函数conv()
dnl=downspl(cvl); % 通过下抽样求出平均部分的分解系数
cvh=conv(cA,hpd); % 高通滤波
dnh=downspl(cvh); % 通过下抽样求出本层分解后的细节部分系数
cA=dnl; % 下抽样后的平均部分系数进入下一层分解
cD=[cD,dnh]; % 将本层分解所得的细节部分系数存入序列cD
end
%% 子函数 %%
function [LL,HL,LH,HH]=mydwt2(x);
% 函数 MYDWT2() 对输入的r*c维矩阵 x 进行二维小波分解，输出四个分解系数子矩阵[LL,HL,LH,HH]
% 输入参数：x —— 输入矩阵，为r*c维矩阵。
% 输出参数：LL,HL,LH,HH —— 是分解系数矩阵的四个相等大小的子矩阵,大小均为 r/2 * c/2 维
% LL：低频部分分解系数； HL：垂直方向分解系数；
% LH：水平方向分解系数； HH：对角线方向分解系数。
lpd=[1/2 1/2];hpd=[-1/2 1/2]; % 默认的低通、高通滤波器
[row,col]=size(x); % 读取输入矩阵的大小
for j=1:row % 首先对输入矩阵的每一行序列进行一维离散小波分解
tmp1=x(j,:);
[ca1,cd1]=mydwt(tmp1,lpd,hpd,1);
x(j,:)=[ca1,cd1]; % 将分解系数序列再存入矩阵x中，得到[L|H]
end
for k=1:col % 再对输入矩阵的每一列序列进行一维离散小波分解
tmp2=x(:,k);
[ca2,cd2]=mydwt(tmp2,lpd,hpd,1);
x(:,k)=[ca2,cd2]; % 将分解所得系数存入矩阵x中，得到[LL,Hl;LH,HH]
end
LL=x(1:row/2,1:col/2); % LL是矩阵x的左上角部分
LH=x(row/2+1:row,1:col/2); % LH是矩阵x的左下角部分
HL=x(1:row/2,col/2+1:col); % HL是矩阵x的右上角部分
HH=x(row/2+1:row,col/2+1:col); % HH是矩阵x的右下角部分
%% 子函数 %%
function y=downspl(x);
%%%%%% 子函数 %%%%%%
% 函数 Y=DOWMSPL(X) 对输入序列进行下抽样，输出序列 Y。
% 下抽样是对输入序列取其偶数位，舍弃奇数位。例如 x=[x1,x2,x3,x4,x5]，则 y=[x2,x4].
N=length(x); % 读取输入序列长度
M=floor(N/2); % 输出序列的长度是输入序列长度的一半（带小数时取整数部分）
i=1:M;
y(i)=x(2*i); function y = lowfrefus(A,B);
% 功能：对输入的小波分解系数矩阵，根据融合算法，得出融合图像的低频小波分解系数
%求出分解系数矩阵的行列数
[row,col]=size(A);
% alpha是方差匹配度比较的阈值
alpha=0.5;
% 根据低频融合算法，先求出矩阵A,B中以点P为中心的区域方差和方差匹配度
for i=1:row
for j=1:col
% 再根据方差匹配度与阈值的比较确定融合图像的小波分解系数
[m2p(i,j),Ga(i,j),Gb(i,j)] = area_var_match(A,B,[i,j]);
Wmin=0.2-0.5*((1-m2p(i,j))/(1-alpha));
Wmax=1-Wmin;
% m2p表示方差匹配度
if m2p(i,j)<alpha
if Ga(i,j)>=Gb(i,j)
% 若匹配度小于阈值，则取区域方差大的相应点的分解系数作为融合图像的分解系数
y(i,j)=A(i,j);
else
y(i,j)=B(i,j);
end
% 若匹配度大于阈值，则采取加权平均方法得出相应的分解系数
else
if Ga(i,j)>=Gb(i,j)
y(i,j)=Wmax*A(i,j)+Wmin*B(i,j);
else
y(i,j)=Wmin*A(i,j)+Wmax*B(i,j);
end
end
end
end
%% 子函数 %%
function w = weivec(x,p)
%功能：对输入的r*c矩阵，计算出以点p为中心时矩阵各点的对应权值
% 距离点p越近，权值就越大。权值是通过行和列的高斯分布加权相加得到的
[r,c]=size(x);
p1=p(1); p2=p(2);
sig=1;
for i=1:r
for j=1:c
w(i,j)=0.5*(gaussmf(i,[sig p1])+gaussmf(j,[sig p2]));
end
end
%% 子函数 %%
function [m2p,Ga,Gb] = area_var_match(A,B,p)
% 功能：计算两个输入矩阵以点p为中心的区域方差以及区域方差匹配度
% 设置区域的大小
level=1;
[subA,mpa,npa]=submat(A,p,level);
% submat 函数取输入矩阵中以点P为中心、阶数为（2*level+1）的方阵作为子矩阵
[subB,mpb,npb]=submat(B,p,level);
[r,c]=size(subA);
% 获取子矩阵的权值分布
w=weivec(subA,[mpa npa]);
% 计算子矩阵的平均值
averA=sum(sum(subA))/(r*c);
averB=sum(sum(subB))/(r*c);
% 计算子矩阵的区域方差
Ga=sum(sum(w.*(subA-averA).^2));
Gb=sum(sum(w.*(subB-averB).^2));
% 计算两个子矩阵的区域方差匹配度
if (Ga==0)&(Gb==0)
m2p=0;
else
m2p=2*sum(sum(w.*abs(subA-averA).*abs(subB-averB)))/(Ga+Gb);
end
%% 子函数 %%
function [smat,mp,np] = submat(x,p,level);
% 功能：取输入矩阵中以点P为中心、阶数为（2*level+1）的方阵作为输出的子矩阵
[row,col]=size(x);
m=p(1); n=p(2);
if (m>row)||(n>col)
error('Point p is out of matrix X !');
return;
end
if ((2*level+1)>row)||((2*level+1)>col)
error('Too large sample area level !');
return;
end
% 设置子矩阵的边界值
up=m-level;
down=m+level;
left=n-level;
right=n+level;
% 若子矩阵的某一边界值超出输入矩阵的相应边界，就进行边界处理，
% 即超出边界后往相反方向平移，使其恰好与边界重合
if left<1
right=right+1-left;
left=1;
end
if right>col
left=left+col-right;
right=col;
end
if up<1
down=down+1-up;
up=1;
end
if down>row
up=up+row-down;
down=row;
end
% 获取作为输出的子矩阵，并计算点p在输出的子矩阵中的位置
smat = x(up:down,left:right);
mp=m-up+1;np=n-left+1;原始图像
融合后
基于拉普拉斯金字塔变换的图像融合
clear all;close all;clc;
I1=imread('clock1.bmp');
I2=imread('clock2.bmp');
I1=double(I1);
I2=double(I2);
Y= fuse_lap(I1, I2, 2, 1, 1);
figure();
subplot(1,2,1);imshow(uint8(I1));
subplot(1,2,2);imshow(uint8(I2));
figure();imshow(uint8(Y));function Y = fuse_lap(M1, M2, zt, ap, mp)
% 功能：基于拉普拉斯金字塔对输入的两幅灰度图像进行融合
% 输入：M1 – 输入的灰度图像A M2 – 输入的灰度图像B
% zt – 最大分解层数 ap – 高通滤波器系数选择 (见 selc.m)
% mp – 基本图像系数选择 (见 selb.m)
%输出：% Y - 融合后的图像
% 检验输入的图像大小是否相同
[z1 s1] = size(M1);
[z2 s2] = size(M2);
if (z1 ~= z2) | (s1 ~= s2)
error('Input images are not of same size');
end;
% 定义滤波器
w = [1 4 6 4 1] / 16;
E = cell(1,zt);
for i1 = 1:zt
% 计算并存储实际图像尺寸
[z s] = size(M1);
zl(i1) = z; sl(i1) = s;
% 检验图像是否需要扩展
if (floor(z/2) ~= z/2), ew(1) = 1; else, ew(1) = 0; end;
if (floor(s/2) ~= s/2), ew(2) = 1; else, ew(2) = 0; end;
% 如果需要扩展的话对图像进行扩展
if (any(ew))
M1 = adb(M1,ew);
M2 = adb(M2,ew);
end;
% 进行滤波
G1 = conv2(conv2(es2(M1,2), w, 'valid'),w', 'valid');
G2 = conv2(conv2(es2(M2,2), w, 'valid'),w', 'valid');
M1T = conv2(conv2(es2(undec2(dec2(G1)), 2), 2*w, 'valid'),2*w', 'valid');
M2T = conv2(conv2(es2(undec2(dec2(G2)), 2), 2*w, 'valid'),2*w', 'valid');
E(i1) = {selc(M1-M1T, M2-M2T, ap)};
M1 = dec2(G1);
M2 = dec2(G2);
end;
M1 = selb(M1,M2,mp);
for i1 = zt:-1:1
M1T = conv2(conv2(es2(undec2(M1), 2), 2*w, 'valid'), 2*w', 'valid');
M1 = M1T + E{i1};
M1 = M1(1:zl(i1),1:sl(i1));
end;
Y = M1;
function Y = selb(M1, M2, mp)
% 功能：对基本图像选择系数
% 输入：
% M1 - 系数A
% M2 - 系数B
% mp - 类型选择
% mp == 1: 选择A
% mp == 2: 选择B
% mp == 3: 选择A和B的均值
%输出：
% Y - 融合系数
switch (mp)
case 1, Y = M1;
case 2, Y = M2;
case 3, Y = (M1 + M2) / 2;
otherwise, error('unknown option');
end;
function Y = selc(M1, M2, ap)
%功能：选择高通滤波器系数
%输入：M1 - 系数A M2 - 系数B mp - 类型选择（输入1，2，3，4）
%输出：Y - 融合系数
% 检验输入的图像大小是否相同
[z1 s1] = size(M1);
[z2 s2] = size(M2);
if (z1 ~= z2) | (s1 ~= s2)
error('Input images are not of same size');
end;
% 方法选择
switch(ap(1))
case 1,
mm = (abs(M1)) > (abs(M2));
Y = (mm.*M1) + ((~mm).*M2);
case 2,
um = ap(2); th = .75;
S1 = conv2(es2(M1.*M1, floor(um/2)), ones(um), 'valid');
S2 = conv2(es2(M2.*M2, floor(um/2)), ones(um), 'valid');
MA = conv2(es2(M1.*M2, floor(um/2)), ones(um), 'valid');
MA = 2 * MA ./ (S1 + S2 + eps);
m1 = MA > th; m2 = S1 > S2;
w1 = (0.5 - 0.5*(1-MA) / (1-th));
Y = (~m1) .* ((m2.*M1) + ((~m2).*M2));
Y = Y + (m1 .* ((m2.*M1.*(1-w1))+((m2).*M2.*w1) + ((~m2).*M2.*(1-w1))+((~m2).*M1.*w1)));
case 3,
um = ap(2);
A1 = ordfilt2(abs(es2(M1, floor(um/2))), um*um, ones(um));
A2 = ordfilt2(abs(es2(M2, floor(um/2))), um*um, ones(um));
mm = (conv2((A1 > A2), ones(um), 'valid')) > floor(um*um/2);
Y = (mm.*M1) + ((~mm).*M2);
case 4,
mm = M1 > M2;
Y = (mm.*M1) + ((~mm).*M2);
otherwise,
error('unkown option');
end;
function Y = dec2(X);
% 功能： 以步长2降采样
% 输入： X – 输入的灰度图像
% 输出： Y – 采样后的图像
[a b] = size(X);
Y = X(1:2:a, 1:2:b);
function Y = undec2(X)
% 功能：以步长2升采样
% 输入：X – 需要升采样的二维图像
% 输出：Y – 升采样后的图像
[z s] = size(X);
Y = zeros(2*z, 2*s);
Y(1:2:2*z,1:2:2*s) = X;
function Y = es2(X, n)
% 功能：将图像矩阵进行扩展
% 输入：X - 输入的二维图像矩阵 n - 要扩展的行数和列数
% 输出：Y - 扩展后的矩阵
[z s] = size(X);
Y = zeros(z+2*n, s+2*n);
Y(n+1:n+z,n:-1:1) = X(:,2:1:n+1);
Y(n+1:n+z,n+1:1:n+s) = X(:,:);
Y(n+1:n+z,n+s+1:1:s+2*n) = X(:,s-1:-1:s-n);
Y(n:-1:1,n+1:s+n) = X(2:1:n+1,:);
Y(n+z+1:1:z+2*n,n+1:s+n) = X(z-1:-1:z-n,:); 原始图像：
融合后：
此博客均属原创或译文，欢迎转载但请注明出处
GithubPage:https://zhangquan1995.github.io
上一篇： HIT 软件构造Lab1 思考
下一篇： 【图像处理与分析】图像分割方法总结