关于Divide by zero

yuyizhen2004

大虾帮忙看看结果为什么那样，谢谢了。下面是一个先做卷积模拟成像，然后去卷积恢复图像的程序

A=zeros(128,128);
A(49:72,36:59)=7;
A(57:80,69:92)=9; %A相当于被成像的物体。

[x,y]=meshgrid(1:128,1:128);
a=1;
c=5;
b_x=64;
b_y=64;
gauss=zeros(128,128);
gauss=a*exp(-((x-b_x).^2+(y-b_y).^2)/c^2); %建立高斯函数，模拟成像系统的点扩展函数。

C=conv2(A,gauss,'same'); %卷积模拟成像，C为A的像，像由物A与高斯卷积，说明像只有几何因素造成的模糊，没有噪声。

rC=fftshift(ifft2(fft2(C)./fft2(gauss))); %去卷积试图恢复C为A。这个式子的数学根据是：如果a=b卷积c，则a的傅立叶变换=b的傅立叶变换*c的傅立叶比那换。
Warning: Divide by zero. %当高斯的半高宽c=2时，没有该警告。高斯的傅立叶变换fft2(gauss）也是一个高斯分布，gauss平扁，那么fft2(gauss）就尖锐，所以c增大，fft2(gauss）绝大多数元素会因为太小而存为零。

rC=fftshift(ifft2(vpa(fft2(C),20)./vpa(fft2(gauss),20))); %所以这里试图提高精度（gauss的傅立叶变换在c=5时，其元素最小数量级-17）。
a=abs(rC);
imshow(a/max(max(a))); %a不显示，一查是非数。
max(max(a))
ans =
   NaN

请问该怎么办？我想了、弄了两三个星期也没搞明白，多谢指点啊。

[ 本帖最后由 sigma665 于 2008-5-27 20:48 编辑 ]

sigma665

vpa不能提高计算精度，只能提高显示精度

另外，分母为0，可以加个eps

yuyizhen2004

谢谢楼上指点。
分母加eps后就不会有被0除的问题，但是imshow(rC/max(max(rC)),'InitialMagnification','fit');title('recoveC');图像会出现“伪影”，当c小于3是不会，c大于3之后就会。

一开始c=3出现“伪影”时（c=3还不至于分母为0），我以为是舍入误差造成的——现在是不是我也不明白了。请问你认为会是什么原因呢？

yuyizhen2004

我一开始的代码是这样的，能否帮我运行看看？
A=zeros(128,128);
A(49:72,36:59)=7;
A(57:80,69:92)=9;
[x,y]=meshgrid(1:128,1:128);
a=1;
c=2;
b_x=64;
b_y=64;
gauss=zeros(128,128);
gauss=a*exp(-((x-b_x).^2+(y-b_y).^2)/c^2);
C=conv2(A,gauss,'same');
rC=fftshift(ifft2(fft2(C)./fft2(gauss)));  %当c>3.9后会devide by zero，该句改为 rC=fftshift(ifft2(fft2(C)./（fft2(gauss)+eps))); ??
figure
subplot(131);imshow(A/max(max(A)),'InitialMagnification','fit');title('A');
subplot(132);imshow(C/max(max(C)),'InitialMagnification','fit');title('C=conv2(A,gauss)');
subplot(133);imshow(rC/max(max(rC)),'InitialMagnification','fit');title('recoverC');

[ 本帖最后由 yuyizhen2004 于 2008-5-28 18:56 编辑 ]