登录论坛

查看完整版本 : 求解地理加权相关系数


实况高手
2009-09-14, 17:02
小弟最近在编写一个求解两幅影像之间地理加权相关系数的算法,
程序如下:function res=forth(name1,name2)
%res 返回方差数组
imgData1=imread(name1);
imgData1=double(imgData1)+1; %把图像数据阵转换为双精度类型矩阵

imgData2=imread(name2);
imgData2=double(imgData2)+1; %把图像数据阵转换为双精度类型矩阵

if ndims(imgData1)>2
imgData1=imgData1(:,:,1); %如果图像为6层tiff格式文件,则提取第一层数据
end

if ndims(imgData2)>2
imgData2=imgData2(:,:,1); %如果图像为6层tiff格式文件,则提取第一层数据
end

s1=size(imgData1);
Width1=s1(1); %图像的像素高
Height1=s1(2); %图像的像素宽

s2=size(imgData2);
Width2=s2(1); %图像的像素高
Height2=s2(2); %图像的像素宽

if Width1~=Width2 && Height1~=Height2
mexErrMsgTxt('the two image must be the same size!');
end

m=3; %定义小窗口的宽
n=3; %定义小窗口的高
m1=fix(m/2); %fix为取整函数
n1=fix(n/2);


h=m/2.0; %%定义小窗口的bandwidth
W=zeros(m,n); %权系数矩阵
Q1=zeros(Width1,Height1); %方差矩阵
Q2=zeros(Width1,Height1);

%计算权系数矩阵
for a=1:m;
for b=1:n;
if a==m1+1 && b==n1+1
D=h;
else
D=(m1+1-a)*(m1+1-a)+(n1+1-b)*(n1+1-b);
D=sqrt(D);
end
if D<h
x=D/h;
x=x*x;
W(a,b)=(1-x)*(1-x);
else
W(a,b)=0;
end
D=0;
end
end
%用小窗口扫描图像,并求Q
for i=m1+1:Width1-m1;
for j=n1+1:Height1-n1;
sum1=0.0;
sum2=0.0;
sum3=0.0;
sum4=0.0;
sum5=0.0;
sum6=0.0;
%地理加权局部统计
for t1=-m1:m1;
for t2=-n1:n1;
sum1=sum1+W(t1+m1+1,t2+n1+1);
sum2=sum2+imgData1(i+t1,j+t2)*W(t1+m1+1,t2+n1+1);
sum4=sum4+imgData2(i+t1,j+t2)*W(t1+m1+1,t2+n1+1);
end
end
u1=sum2/sum1;
u2=sum4/sum1;

for t1=-m1:m1;
for t2=-n1:n1;
y1=imgData1(i+t1,j+t2)-u1;
y2=imgData2(i+t1,j+t2)-u2;
y1=y1*y1;
y2=y2*y2;
sum3=sum3+W(t1+m1+1,t2+n1+1)*y1;
sum5=sum5+W(t1+m1+1,t2+n1+1)*y2;
end
end
%方差
Q1(i,j)=sum3/sum1;
Q2(i,j)=sum5/sum1;



for t1=-m1:m1;
for t2=-n1:n1;
y1=imgData1(i+t1,j+t2)-u1;
y2=imgData2(i+t1,j+t2)-u2;
y=y1*y2;

sum6=sum6+W(t1+m1+1,t2+n1+1)*y;
if Q1(i,j)~=0 && Q2(i,j)~=0
res(i-1,j-1)=sum6/(sqrt(Q1(i,j))*sqrt(Q2(i,j)));
end
end
end
end
end
它是用不同大小的窗口去扫描影像上所有的像素点,运行后出了问题,求出来的结果是一个矩阵,但是系数大于1,我反复检查了很多遍,觉得公式没错,也找不出我的算法哪里有错,请高手帮忙调试一下,可以自己找两幅影像计算一下,在下不胜感激。