ramp滤波函数matlab,matlab - MATLAB如何在频域中实现Ram-Lak滤波器(斜坡滤波器)? - 堆栈内存溢出...

如果您想在傅立叶域中进行无需过滤的反Radon变换,则列出的公式是中间结果。 另一种方法是使用空间域中的卷积来完成整个滤波反投影算法,这可能在40年前更快; 你最终会重新发布你发布的公式。 但是,我现在不建议,特别是不适合你的第一次重建; 你应该先了解希尔伯特变换。

无论如何,这里有一些Matlab代码,它执行强制性的Shepp-Logan模拟滤波反投影重建。 我将展示如何使用Ram-Lak过滤器进行自己的过滤。 如果我真的很有动力,我会用一些interp2命令和总结来代替radon / iradon。

phantomData=phantom();

N=size(phantomData,1); theta = 0:179; N_theta = length(theta); [R,xp] = radon(phantomData,theta); % make a Ram-Lak filter. it's just abs(f). N1 = length(xp); freqs=linspace(-1, 1, N1).'; myFilter = abs( freqs ); myFilter = repmat(myFilter, [1 N_theta]); % do my own FT domain filtering ft_R = fftshift(fft(R,[],1),1); filteredProj = ft_R .* myFilter; filteredProj = ifftshift(filteredProj,1); ift_R = real(ifft(filteredProj,[],1)); % tell matlab to do inverse FBP without a filter I1 = iradon(ift_R, theta, 'linear', 'none', 1.0, N); subplot(1,3,1);imagesc( real(I1) ); title('Manual filtering') colormap(gray(256)); axis image; axis off % for comparison, ask matlab to use their Ram-Lak filter implementation I2 = iradon(R, theta, 'linear', 'Ram-Lak', 1.0, N); subplot(1,3,2);imagesc( real(I2) ); title('Matlab filtering') colormap(gray(256)); axis image; axis off % for fun, redo the filtering wrong on purpose % exclude high frequencies to create a low-resolution reconstruction myFilter( myFilter > 0.1 ) = 0; ift_R = real(ifft(ifftshift(ft_R .* myFilter,1),[],1)); I3 = iradon(ift_R, theta, 'linear', 'none', 1.0, N); subplot(1,3,3);imagesc( real(I3) ); title('Low resolution filtering') colormap(gray(256)); axis image; axis off

aHR0cHM6Ly9pLnN0YWNrLmltZ3VyLmNvbS9PMnFTZy5naWY=

  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
根据需要进行修改: 在MATLAB,反投影算法是基于滤波投影数据的非常有用的方法之一。Ramp滤波器的主要思想是消除投影数据的高频噪声,以获得更清晰和更准确的重建图像。以下是Ramp滤波器反投影算法的步骤: 1. 读取投影数据,并将其转换为极坐标形式。 2. 在极坐标对投影数据进行Ramp滤波,并将其转换回直角坐标。 3. 对过滤后的投影数据执行反投影。 4. 将所有反投影图像叠加在一起,以获得最终重建图像。 下面是MATLAB代码的示例,演示如何使用Ramp滤波器进行反投影: % 读取投影数据 proj_data = read_projection_data(filename); % 将投影数据转换为极坐标形式 [theta, r] = meshgrid(linspace(0, 2*pi, size(proj_data, 1)), ... linspace(-1, 1, size(proj_data, 2))); x = r .* cos(theta); y = r .* sin(theta); % 进行Ramp滤波 ramp_filter = 2 .* abs(linspace(0, 1, size(proj_data, 1)) - 0.5); filt_proj_data = fft(proj_data) .* repmat(ramp_filter', 1, size(proj_data, 2)); filt_proj_data = real(ifft(filt_proj_data)); % 将过滤后的投影数据转换回直角坐标 [x_grid, y_grid] = meshgrid(linspace(-1, 1, size(proj_data, 2))); filt_proj_interp = interp2(x, y, filt_proj_data, x_grid, y_grid); % 执行反投影 back_proj = backprojection(filt_proj_interp, theta); % 叠加所有反投影图像,获得最终重建图像 recon = sum(back_proj, 3); 以上代码可以根据具体的数据和需求进行修改。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值