clear
m=1;
n=1;
I1=imread（‘before.tif‘）;
I1=double（I1）;
I2=imread（‘after.tif‘）;
I2=double（I2）;

[M?N]=size（I1）;
f1ave=zeros（MN）;
f2ave=zeros（MN）;

for?i=1:M
????for?j=1:N
??????????????f1sum=0.0;
??????????????f2sum=0.0;
??????????????nn=0;
????????
????????for?k=-m:m
????????????for?l=-n:n
????????????????if?（（i+k）<1?||?（j+l）<1?||?（i+k）>M?||?（j+l）>N）
????????????????????f1sum=f1sum;
????????????????????f2sum=f2sum;
????????????????else
????????????????????f1sum=f1sum+I1（i+kj+l）;
????????????????????f2sum=f2sum+I2（i+kj+l）;
????????????????????nn=nn+1;
????????????????end
????????????end
????????end
????????
????????f1ave（ij）=f1sum/nn;
????????f2ave（ij）=f2sum/nn;
????end
end

f1dif=I1-f1ave;
f2dif=I2-f2ave;

for?i=1:M
????for?j=1:N
????????cmult=0.0;
????????c1=0.0;
????????c2=0.0;
????????
?????????for?k=-m:m
????????????for?l=-n:n
????????????????if?（（i+k）<1?||?（j+l）<1?||?（i+k）>M?||?（j+l）>N）
????????????????????cmult=cmult;
????????????????????c1=c1;
????????????????????c2=c2;
????????????????else
????????????????????cmult=cmult+f1dif（i+kj+l）*f2dif（i+kj+l）;
????????????????????c1=c1+f1dif（i+kj+l）*f1dif（i+kj+l）;
????????????????????c2=c2+f2dif（i+kj+l）*f2dif（i+kj+l）;
????????????????end
????????????end
?????????end
????????cmultave（ij）=cmult/（2*m+1）/（2*n+1）;
????????c1ave（ij）=c1/（2*m+1）/（2*n+1）;
????????c2ave（ij）=c2/（2*m+1）/（2*n+1）;
????end
end

cxie=cmultave./sqrt（c1ave.*c2ave）;
I=255.*（1-abs（cxie））;
I=uint8（I）;
imshow（I）;

imwrite（I‘Isub.tif‘）;
????????????????

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????目錄???????????0??2007-02-28?09:00??XieFangXiangGuan\
?????文件?????1536528??2008-01-03?15:02??XieFangXiangGuan\after.tif
?????文件?????1536528??2008-01-03?15:02??XieFangXiangGuan\before.tif
?????文件?????1536364??2008-01-03?15:10??XieFangXiangGuan\Isub.tif
?????文件????????9728??2007-02-28?09:00??XieFangXiangGuan\Thumbs.db
?????文件????????1572??2007-02-28?09:02??XieFangXiangGuan\XieFangXiangGuan.m