18 SAR图像和光学图像的配准算法（matlab程序）-编程知识

18 SAR图像和光学图像的配准算法（matlab程序）

1.简述

合成孔径雷达(synthetic aperture radar,SAR)图像配准的主要目标是对同一或不同传感器在不同时间、不同视点捕获的SAR图像进行配准。SAR因具有全天候成像能力和地物穿透能力，因此具有非常广泛的应用，如变化检测[1]、图像融合[2]、目标检测与识别[3]。

图像配准方法可分为两类:基于区域的配准方法和基于特征的配准方法[4]。基于区域的配准方法对噪声敏感，当目标发生形变、旋转、尺寸变换时容易导致算法失败且计算量过大。SAR图像中通常包含许多不同的特征，因此基于特征的配准方法具有较高的计算效率[5]。对SAR图像匹配而言，更具研究意义和实用价值。

基本方案：
(1)基于边缘提取的方法，即提取 SAR 图像和光学图像边缘特征的基础上按照一定的相似性测度进行匹配。缺点是依赖图像边缘的准确提取，如果边缘提取不够准确则配准效果就受到影响。
(2)基于封闭区域的方法，采用各种分割的方法提取封闭区域或子区，在闭合区域边缘进行匹配，精度和基于边缘类似，但是必须要求图像上存在封闭的边缘，如果没有封闭的边缘则很难实现配准。
(3)基于图像互相关的方法，但是由于 SAR 和光学机理差异大，互相关匹配的方法并不能完全有效。无论是提取图像的边缘还是封闭区域，都是根据原始图像提取边缘，而边缘的提取又会受到局部阈值的影响，阈值的不同会得到不同的边缘，因此基于边缘的方法十分依赖于边缘的提取效果。

2.代码

clear;clc;
%% 基于各向异性热扩散方程的SAR图像分割
Img = imread('6.bmp');
Img = rgb2gray(Img);
figure;
imshow(Img);
Img = double(Img);
Img = exp(Img/17); %11 6

SIZE = size(Img);
SIZEX = SIZE(1,1);
SIZEY = SIZE(1,2);

Total = 0;time = 40;
%% 设定初始阈值
Vector_1 = Img(Img<3);Vector_1 = Vector_1.^2;%2,180 7 3
Vector_2 = Img(Img>=3 & Img<=254);Vector_2 = Vector_2.^2;
Vector_3 = Img(Img>254);Vector_3 = Vector_3.^2;
Mark_1 = size(Vector_1);Mark_2 = size(Vector_2);Mark_3 = size(Vector_3);
P_Class(:,:,1) = Mark_1(1,2)/(SIZEX*SIZEY)*ones(SIZEX,SIZEY);
P_Class(:,:,2) = Mark_2(1,2)/(SIZEX*SIZEY)*ones(SIZEX,SIZEY);
P_Class(:,:,3) = Mark_3(1,2)/(SIZEX*SIZEY)*ones(SIZEX,SIZEY);

Sigma(1) = mean(Vector_1);
Sigma(2) = mean(Vector_2);
Sigma(3) = mean(Vector_3);
Mark_1 = 1;Mark_2 = 1;Mark_3 = 1;
disp('准备开始初始分割......');
for m = 1:time
Vector_1 = 0;Vector_2 = 0;Vector_3 = 0;k = 1:3;
str = '初始分割已经完成';
str = strcat(str,num2str(m/time*100.0),'%');
disp(str);
for i =1:SIZEX
for j = 1:SIZEY
Total = 0;
for k = 1:3
Total = Total+1/(Sigma(k)+0.00000001)*exp(-1*Img(i,j)^2/(Sigma(k)+0.00000001))*P_Class(i,j,k);
end

for k = 1:3

P_Class_tem(i,j,k) = 1/(Sigma(k)+0.00000001)*exp(-1*Img(i,j)^2/(Sigma(k)+0.00000001))*P_Class(i,j,k)/(Total+1/inf);
end
end
end
P_Class = P_Class_tem;
[Value,Flag] = max(P_Class,[],3);
Img1 = Img.^2;
Vector_1 = Img1(Flag==1);
Vector_2 = Img1(Flag==2);
Vector_3 = Img1(Flag==3);

Sigma(1) = std(Vector_1,1);
Sigma(2) = std(Vector_2,1);
Sigma(3) = std(Vector_3,1);
Mark_1 = 1;Mark_2 = 1;Mark_3 = 1;
end
figure, imshow(Flag,[]);
P_Class = abs(1000*P_Class);
disp('初始分割已经完成');
disp('开始平滑后验概率矩阵......');
for i = 1:3
if i == 3
P_Class3(:,:,i) = RSRADFilter(P_Class(:,:,i),1.06);
else
P_Class3(:,:,i) = RSRADFilter(P_Class(:,:,i),2.36);
end
end
for i = 1:SIZEX
for j = 1:SIZEY
[Value,Flag(i,j)] = max([P_Class3(i,j,1),P_Class3(i,j,2),P_Class3(i,j,3)]);
end
end
Flag = WipeSpeckles(Flag);
figure;
imshow(Flag,[]);