MT6.581/20.482J 生物系统中的算法与计算技术基础 Bruce Tidor教授 Jacob K.White教授 MT6.581/20.482J问题集#6 2006.5.18,星期三，5pm 1.噪声图像数据恢复。图像系统会模糊对象数据并加入噪声，因此我们的 图像采集装置获得的信号如下 8=Hf°+w (1) 其中f·是原始对象数据，H是点扩散函数，w是添加的噪声，g是采集 的图像数据。在这个问题中，我们将探讨使用奇异值分解(SVD)在不过 度放大噪声w的情况下尽可能的恢复原始信号°。 (a)运行脚本load data.m初始化对象数据f,PSFh,以及随机白噪声w。 回忆我们在课堂上讲的，h关于0]对称，并且以长度n为周期重复， 即hn-1]=h1],hn-2]=[2]等。请生成循环行列式矩阵H并计算g。 (b)对H进行奇异值分解绘制U和V的前面几列。对结果进行解释。 (C)绘制H的奇异值可看到它们都是成对出现。所有的循环行列式矩阵 都展现出这种特征，为什么？ (d结果图像数据的盲反卷积为 f=H-g (2) 请绘制出结果。象预料的一样，反卷积放大了噪声淹没了原始对象 数据。计算H的SVD得到 H=U∑VI (3) H=Σ-U (4) 其中Σ是对角矩阵，∑：表示将输入基向量'，映射为输出基向量U 时H的伸展度。如果Σ：非常小，H的逆会放大U中的元素。使用 partial inv.m只会在奇异值∑较大的时候才反转H的方向。请试着 改变阈值大小并绘制结果数据。恢复原始图像的最佳阈值是多少？ (e)生成2到3个其他的信号。象在(a)中一样加入噪声数据来模糊信号。 使用(d)中得到的最优阈值来恢复原始信号。试着改变阈值一从(d) 中得到的答案仍然是最优的吗？ 2.不连续建模。在课堂上，我们选择傅里叶级数来表示图像数据。这个问 题将探讨这种表示方法带来的一个后果：所谓吉布斯现象导致的图像数
MIT 6.581/20.482J 生物系统中的算法与计算技术基础 Bruce Tidor 教授 Jacob K. White 教授 1 MIT 6.581/20.482J 问题集 #6 2006.5.18，星期三，5pm 1. 噪声图像数据恢复。图像系统会模糊对象数据并加入噪声，因此我们的 图像采集装置获得的信号如下 g = Hf + w 0 (1) 其中 f 0 是原始对象数据，H 是点扩散函数，w 是添加的噪声，g 是采集 的图像数据。在这个问题中，我们将探讨使用奇异值分解(SVD)在不过 度放大噪声 w 的情况下尽可能的恢复原始信号 f 0 。 (a) 运行脚本 load_data.m 初始化对象数据 f，PSF h，以及随机白噪声 w。 回忆我们在课堂上讲的，h 关于 h[0]对称，并且以长度 n 为周期重复， 即 h[n-1]=h[1]，h[n-2]=h[2]等。请生成循环行列式矩阵 H 并计算 g。 (b) 对 H 进行奇异值分解绘制 U 和 V 的前面几列。对结果进行解释。 (c) 绘制 H 的奇异值可看到它们都是成对出现。所有的循环行列式矩阵 都展现出这种特征，为什么？ (d) 结果图像数据的盲反卷积为 f H g −1 = (2) 请绘制出结果。象预料的一样，反卷积放大了噪声淹没了原始对象 数据。计算 H 的 SVD 得到 T H = UΣV (3) T H V U −1 −1 = Σ (4) 其中 Σ 是对角矩阵，Σi,i 表示将输入基向量 Vi 映射为输出基向量 Ui 时 H 的伸展度。如果 Σi,i非常小，H 的逆会放大 Ui 中的元素。使用 partial_inv.m 只会在奇异值 Σi,i较大的时候才反转 H 的方向。请试着 改变阈值大小并绘制结果数据。恢复原始图像的最佳阈值是多少？ (e) 生成 2 到 3 个其他的信号。象在(a)中一样加入噪声数据来模糊信号。 使用(d)中得到的最优阈值来恢复原始信号。试着改变阈值——从(d) 中得到的答案仍然是最优的吗？ 2. 不连续建模。在课堂上，我们选择傅里叶级数来表示图像数据。这个问 题将探讨这种表示方法带来的一个后果：所谓吉布斯现象导致的图像数
MT6.581/20.482J 生物系统中的算法与计算技术基础 Bruce Tidor教授 Jacob K.White教授 据的不连续性。 (a)使用load data.m中的方脉冲函数x,并通过f找到X的傅里叶变换。 在仁10,20,40的地方截断傅里叶级数，其中k是非0傅里叶系数 的个数。绘制逆转换图.注意截断操作要小心（它并不是直接通过SVD 建立一个低阶矩阵)。结果信号应该超过原始的不连续性。增加k值 有什么会变化？而什么又不变呢？ (b)现在采用问题1中的对象数据。转换这个平滑的方脉冲函数，并采 用()中相同的k值来截断表示。当将这些截断的表示数据逆转换时 会出现什么现象？ (c)如果采集锐化边缘很重要，那么采用类小波变换来表示数据可能更合 适。从loadWaveletBasis.m载入小波基Q。将信号x和xf加入这个基， 并在=k1,2,时截断表示，则 x)=Q(:,L,:k)*Q(:L,:k)*x (5) 如果我们用截短的小波级数来表示这个函数会有什么变化？绘制结 果并说明在哪些情况下使用小波表示更有优势，而哪些情况下傅里 叶表示更好？ 3.二维图像处理。在这个问题中，我们将探讨几种不同的方法来处理对象 和点扩散函数(PSF)之间的卷积。 ()图1显示了一个简单的对象以及其在PSF卷积下的图像。两个矩阵 中的非0元素（即图中非黑色像素对应的部分）分别为： 11 and 你能推出相应的PSF吗？在MATLAB中采用下述语句证实你的答案 abs(ifft2(fft2(Mobi).*fft2(Ms))) 其中每个矩阵包含对应的8×8个像素值。 object image 图18×8像素的对象和图像 2
MIT 6.581/20.482J 生物系统中的算法与计算技术基础 Bruce Tidor 教授 Jacob K. White 教授 2 据的不连续性。 (a) 使用 load_data.m 中的方脉冲函数 x，并通过 fft 找到 X 的傅里叶变换。 在 k=10，20，40 的地方截断傅里叶级数，其中 k 是非 0 傅里叶系数 的个数。绘制逆转换图。注意截断操作要小心(它并不是直接通过 SVD 建立一个低阶矩阵)。结果信号应该超过原始的不连续性。增加 k 值 有什么会变化？而什么又不变呢？ (b) 现在采用问题 1 中的对象数据 f。转换这个平滑的方脉冲函数，并采 用(a)中相同的 k 值来截断表示。当将这些截断的表示数据逆转换时 会出现什么现象？ (c) 如果采集锐化边缘很重要，那么采用类小波变换来表示数据可能更合 适。从 loadWaveletBasis.m 载入小波基 Q。将信号 x 和 xf 加入这个基， 并在 k=k1，k2，k3 时截断表示，则 x Q k Q k x k T (:,1,: ) * (:,1,: ) * ( ) = (5) 如果我们用截短的小波级数来表示这个函数会有什么变化？绘制结 果并说明在哪些情况下使用小波表示更有优势，而哪些情况下傅里 叶表示更好？ 3. 二维图像处理。在这个问题中，我们将探讨几种不同的方法来处理对象 和点扩散函数(PSF)之间的卷积。 (a) 图 1 显示了一个简单的对象以及其在 PSF 卷积下的图像。两个矩阵 中的非 0 元素(即图中非黑色像素对应的部分)分别为： ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ 4 1 4 1 4 1 4 1 4 1 4 1 4 1 4 1 1 1 1 1 1 1 1 1 and 你能推出相应的PSF吗？在MATLAB中采用下述语句证实你的答案 abs(ifft2(fft2(Mobj)⋅*fft2(Mpsf)))) 其中每个矩阵包含对应的 8×8 个像素值。 图 1 8×8 像素的对象和图像
MT6.581/20.482J 生物系统中的算法与计算技术基础 Bruce Tidor教授 Jacob K.White教授 (b)对象与PSF间的二维卷积可用矩阵的形式表示： H(0,) H(n-1,:)…… H(1,)f(1,1) g(1,1) H1,) H(0,) H(2,) H(L,) f(L.n) g(1,n) H(n-1,:)） H(0,)f°(n,n) g(n,n) 其中∫)和g)是()位置像素对应的对象和图像数据。这个矩阵 是“二重循环行列式”，每个块元素也是一个循环行列式 h(k,0) h(k,n-1)…… h(k,1) h(k,) hk,0) h(k,2) h(k,1) H(k,:)= hk,2) h(k,n-1) hk,0) 其中每一项都来自于二维PSF的第i行。采用(a)中找到的PSF建立 该矩阵，并验证其图像数据与先前矩阵相乘获得的数据相同。 (c)循环行列式乘法能通过一个一维FFT来计算.MATLAB中表示如下： H(k.:)*f=ifft(fft(h(k,:))-*fft(f)) 采用一维F℉T来实现二维卷积运算，并用前面获得的相同图像数据 来进行验证。 (d从课程网站上下载mandrill.bmp,使用一维FFT来执行PSF的卷积。 你应该会得到一副与图2类似的图片。这张图片比原图要显得亮但不 真实。规范化PSF,使其不会放大输入信号的能量。重新执行仿真程 序并说明原始数据和卷积数据的区别。 图2狒狒图像以及PFS的卷积图像
MIT 6.581/20.482J 生物系统中的算法与计算技术基础 Bruce Tidor 教授 Jacob K. White 教授 3 (b) 对象与 PSF 间的二维卷积可用矩阵的形式表示： ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − ( , ) (1, ) (1,1) ( , ) (1, ) (1,1) ( 1,:) (0,:) (2,:) (1,:) (1,:) (0,:) (0,:) ( 1,:) (1,:) 0 0 0 g n n g n g f n n f n f H n H H H H H H H n H M M M M M M M O M O O O O L L 其中 f 0 (i,j)和 g(i,j)是(i,j)位置像素对应的对象和图像数据。这个矩阵 是“二重循环行列式”，每个块元素也是一个循环行列式 ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − − = ( , 1) ( ,0) ( ,2) ( ,2) ( ,1) ( ,1) ( ,0) ( ,0) ( , 1) ( ,1) ( ,:) h k n h k h k h k h k h k h k h k h k n h k H k M O M O O O O L L 其中每一项都来自于二维 PSF 的第 i 行。采用(a)中找到的 PSF 建立 该矩阵，并验证其图像数据与先前矩阵相乘获得的数据相同。 (c) 循环行列式乘法能通过一个一维 FFT 来计算。MATLAB 中表示如下： H(k,:)*f = ifft(fft(h(k,:))⋅*fft(f)) 采用一维 FFT 来实现二维卷积运算，并用前面获得的相同图像数据 来进行验证。 (d) 从课程网站上下载 mandrill.bmp，使用一维 FFT 来执行 PSF 的卷积。 你应该会得到一副与图 2 类似的图片。这张图片比原图要显得亮但不 真实。规范化 PSF，使其不会放大输入信号的能量。重新执行仿真程 序并说明原始数据和卷积数据的区别。 图 2 狒狒图像以及 PFS 的卷积图像