登录论坛

查看完整版本 : [MATLAB基础] 有关求pearson相关系数的问题


实况高手
2009-07-16, 16:29
我用matlab程序运行了一个求pearson相关系数的程序,但是得出的结果值却很大,按道理说相关系数值应落在-1到1之间,可能是程序哪里出了点问题,哪位达人帮忙看看解决解决,我是初学者,很多东西都还不懂,谢谢大家。
这是程序:
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);
v=fix(m*n/2);

h=m/2; %%定义小窗口的bandwidth
W=zeros(m,n); %权系数矩阵
Q=zeros(Width1,Height1); %方差矩阵
sum1=0.0;
sum2=0.0;
sum3=0.0;
sum4=0.0;
sum5=0.0;
sum6=0.0;
%计算权系数矩阵
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-1;
for j=n1+1:Height1-1;

%地理加权局部统计
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;
sum2=0.0;
sum4=0.0;
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;
sum3=0.0;
sum5=0.0;


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
p=sum6/(sqrt(Q1(i,j))*sqrt(Q2(i,j)));
end
end
end
end
end
res=p;