声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2868|回复: 3

[图像处理] 关于Divide by zero

[复制链接]
发表于 2008-5-27 20:41 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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

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 编辑 ]
回复
分享到:

使用道具 举报

发表于 2008-5-27 20:49 | 显示全部楼层

回复 楼主 的帖子

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

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

评分

1

查看全部评分

 楼主| 发表于 2008-5-28 18:37 | 显示全部楼层
谢谢楼上指点。
分母加eps后就不会有被0除的问题,但是imshow(rC/max(max(rC)),'InitialMagnification','fit');title('recoveC');图像会出现“伪影”,当c小于3是不会,c大于3之后就会。

一开始c=3出现“伪影”时(c=3还不至于分母为0),我以为是舍入误差造成的——现在是不是我也不明白了。请问你认为会是什么原因呢?
 楼主| 发表于 2008-5-28 18:51 | 显示全部楼层

恢复2楼

我一开始的代码是这样的,能否帮我运行看看?
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 编辑 ]
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-12 05:29 , Processed in 0.057944 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表