图像分割配准与融合

内容：

• 图像目标边界描述
• 图像分割
• 图像配准
• 图像融合

Github个人博客：https://joeyos.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 = ze