ajn25953
2008-05-12, 10:45
正在做光流法运动估计的设计，这段程序是先把两幅图灰度化处理，然后求运动场，最后有运动场和第一幅图恢复出第二幅图像，程序运行后提示Warning: Matrix is singular to working precision.，请高手务必帮助修改一下，万分感谢。
clear
clc
%读图像的值
x=imread('q009.bmp');
x=double(rgb2gray(x));
y=imread('q011.bmp');
y=double(rgb2gray(y));
%显示原图像
figure(1);
subplot(221)
image(x);
colormap(gray(255));
title('第一幅图像');
figure(1);
subplot(222)
image(y);
colormap(gray(255));
title('第二幅图像');
[m,n]=size(x);
%扩充图像边缘
x(m+1,n+1)=0;
y(m+1,n+1)=0;
%有限差分的梯度估计
for i=1:m
for j=1:n
dx1(i,j)=(x(i+1,j)-x(i,j)+x(i+1,j+1)-x(i,j+1)+y(i+1,j)-y(i,j)+y(i+1,j+1)-y(i,j+1))/4;
dx2(i,j)=(x(i,j+1)-x(i,j)+x(i+1,j+1)-x(i+1,j)+y(i,j+1)-y(i,j)+y(i+1,j+1)-y(i+1,j))/4;
dtt(i,j)=(y(i,j)-x(i,j)+y(i+1,j)-x(i+1,j)+y(i,j+1)-x(i,j+1)+y(i+1,j+1)-x(i+1,j+1))/4;
end
end
%计算的运动矢量的数量
block_size=16;
slide_size=8;
p1=m/slide_size;
p=double(int8(p1));
if p<p1
p=p+1;
end
q1=n/slide_size;
q=double(int8(q1));
if q<q1
q=q+1;
end
v1=zeros(p,q);
v2=zeros(p,q);
%画运动场
for mm=1:p
for nn=1:q
a=zeros(2,2);
b=zeros(2,1);
for i=1:block_size
for j=1:block_size
ii=(mm-1)*slide_size+i;
jj=(nn-1)*slide_size+j;
if ii>m|jj>n
dx1(ii,jj)=0;
dx2(ii,jj)=0;
dtt(ii,jj)=0;
end
a(1,1)=a(1,1)+dx1(ii,jj)*dx1(ii,jj);
a(1,2)=a(1,2)+dx1(ii,jj)*dx2(ii,jj);
a(2,2)=a(2,2)+dx2(ii,jj)*dx2(ii,jj);
b(1)=b(1)-dx1(ii,jj)*dtt(ii,jj);
b(2)=b(2)-dx2(ii,jj)*dtt(ii,jj);
end
end
v=inv(a)*b;
if isnan(v(1))
if mm==1|nn==1
v1(mm,nn)=0;
else
v1(mm,nn)=v1(mm-1,nn-1);
end
else
v1(mm,nn)=v(1);
end
if isnan(v(2))
if mm==1|nn==1
v2(mm,nn)=0;
else
v2(mm,nn)=v2(mm-1,nn-1);
end
else
v2(mm,nn)=v(2);
end
end
end
figure(2);
quiver(v1,v2);
axis ij;
title('运动场');
%帧差分的PSNR
x=x(1:m,1:n);
y=y(1:m,1:n);
diff1=mean(mean((y-x).^2));
PSNR1=10*log10(255*255/diff1)
%L-K方法的PSNR
evaluate=x;
for i=1:p
for j=1:q
for M=slide_size*(i-1)+1:slide_size*i
for N=slide_size*(j-1)+1:slide_size*j
d=v1(i,j);
d1=round(d);
d=v2(i,j);
d2=round(d);
l1=M+d1;
l2=N+d2;
if l1<=0
l1=1;
end
if l2<=0
l2=1;
end
if l1>m
l1=m;
end
if l2>n
l2=n;
end
evaluate(l1,l2)=x(M,N);
end
end
end
end
diff2=mean(mean((y-evaluate).^2));
PSNR2=10*log10(255*255/diff2)
%画恢复后的图像
figure(1);
subplot(224)
colormap(gray);
imagesc(evaluate,[0,255]);
title('恢复后的图象');
clear
clc
%读图像的值
x=imread('q009.bmp');
x=double(rgb2gray(x));
y=imread('q011.bmp');
y=double(rgb2gray(y));
%显示原图像
figure(1);
subplot(221)
image(x);
colormap(gray(255));
title('第一幅图像');
figure(1);
subplot(222)
image(y);
colormap(gray(255));
title('第二幅图像');
[m,n]=size(x);
%扩充图像边缘
x(m+1,n+1)=0;
y(m+1,n+1)=0;
%有限差分的梯度估计
for i=1:m
for j=1:n
dx1(i,j)=(x(i+1,j)-x(i,j)+x(i+1,j+1)-x(i,j+1)+y(i+1,j)-y(i,j)+y(i+1,j+1)-y(i,j+1))/4;
dx2(i,j)=(x(i,j+1)-x(i,j)+x(i+1,j+1)-x(i+1,j)+y(i,j+1)-y(i,j)+y(i+1,j+1)-y(i+1,j))/4;
dtt(i,j)=(y(i,j)-x(i,j)+y(i+1,j)-x(i+1,j)+y(i,j+1)-x(i,j+1)+y(i+1,j+1)-x(i+1,j+1))/4;
end
end
%计算的运动矢量的数量
block_size=16;
slide_size=8;
p1=m/slide_size;
p=double(int8(p1));
if p<p1
p=p+1;
end
q1=n/slide_size;
q=double(int8(q1));
if q<q1
q=q+1;
end
v1=zeros(p,q);
v2=zeros(p,q);
%画运动场
for mm=1:p
for nn=1:q
a=zeros(2,2);
b=zeros(2,1);
for i=1:block_size
for j=1:block_size
ii=(mm-1)*slide_size+i;
jj=(nn-1)*slide_size+j;
if ii>m|jj>n
dx1(ii,jj)=0;
dx2(ii,jj)=0;
dtt(ii,jj)=0;
end
a(1,1)=a(1,1)+dx1(ii,jj)*dx1(ii,jj);
a(1,2)=a(1,2)+dx1(ii,jj)*dx2(ii,jj);
a(2,2)=a(2,2)+dx2(ii,jj)*dx2(ii,jj);
b(1)=b(1)-dx1(ii,jj)*dtt(ii,jj);
b(2)=b(2)-dx2(ii,jj)*dtt(ii,jj);
end
end
v=inv(a)*b;
if isnan(v(1))
if mm==1|nn==1
v1(mm,nn)=0;
else
v1(mm,nn)=v1(mm-1,nn-1);
end
else
v1(mm,nn)=v(1);
end
if isnan(v(2))
if mm==1|nn==1
v2(mm,nn)=0;
else
v2(mm,nn)=v2(mm-1,nn-1);
end
else
v2(mm,nn)=v(2);
end
end
end
figure(2);
quiver(v1,v2);
axis ij;
title('运动场');
%帧差分的PSNR
x=x(1:m,1:n);
y=y(1:m,1:n);
diff1=mean(mean((y-x).^2));
PSNR1=10*log10(255*255/diff1)
%L-K方法的PSNR
evaluate=x;
for i=1:p
for j=1:q
for M=slide_size*(i-1)+1:slide_size*i
for N=slide_size*(j-1)+1:slide_size*j
d=v1(i,j);
d1=round(d);
d=v2(i,j);
d2=round(d);
l1=M+d1;
l2=N+d2;
if l1<=0
l1=1;
end
if l2<=0
l2=1;
end
if l1>m
l1=m;
end
if l2>n
l2=n;
end
evaluate(l1,l2)=x(M,N);
end
end
end
end
diff2=mean(mean((y-evaluate).^2));
PSNR2=10*log10(255*255/diff2)
%画恢复后的图像
figure(1);
subplot(224)
colormap(gray);
imagesc(evaluate,[0,255]);
title('恢复后的图象');