CBP

最新推荐文章于 2024-01-08 21:33:02 发布

老李今天学习了吗

最新推荐文章于 2024-01-08 21:33:02 发布

阅读量1k

点赞数

本文链接：https://blog.youkuaiyun.com/weixin_29732003/article/details/103688488

版权

本文详细介绍了如何使用Matlab实现卷积反投影算法，包括选定参数、生成图像、进行投影、构造卷积核、执行卷积及反投影等步骤，并提供了完整的代码示例。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

实验报告

目标

使用matlab实现卷积反投影

实现步骤

1.选定参数
2.做出图像
3.投影
4.构造卷积核（RL滤波器）
5.卷积
7.反投影

实现结果

原图：
在这里插入图片描述
效果图

代码如下

%% Parameters
N = 512;
views = 360;
bins = 512;
dv = 2 * pi / views;
x1 = 0; y1 = 0; r1 = 200; mu1 = 0.5;
x2 = 80; y2 = 80; r2 = 30; mu2 = 1;
%% generateImage
vec = -N/2 + 0.5: N/2 - 0.5;
[xx, yy] = meshgrid(vec, -vec);
img = zeros(N, N);
img((xx - x1).^2 + (yy - y1).^2 <= r1^2) = mu1;
img((xx - x2).^2 + (yy - y2).^2 <= r2^2) = mu2;
figure('name','cirle in circle'); imshow(img, []); title('Original image');
%% projectionAndReconstructionByExistFunction
R = radon(img);
figure('name','projection');imshow(R, []);title('radon projection');
Ir = iradon(R,0:179);
figure('name','Iprojection');imshow(Ir, []);title('Iprojection');
%% projectionByInterproation
dn = 0.5;
inte_vec = -N/2 + dn/2: dn: N/2 - dn/2;
[ox, oy] = meshgrid(vec, -inte_vec);
alpha = dv/2: dv: 2*pi - dv/2;
%开始投影
proj = zeros(views, bins);
for i = 1: views
    theta = alpha(i);
    X = ox .* cos(theta) - oy .* sin(theta);
    Y = ox .* sin(theta) + oy .* cos(theta);
    temp = interp2(xx, yy, img, X, Y, 'linear',0);
    proj(i, :) = sum(temp);
end
figure; imshow(proj, []); title('Projection image');
%% RLkernel
[n,m] = size(proj);
h=zeros(1,(m*2-1));
for i=0:m-1
    if i==0
        h(m-i)=1/4;
    elseif rem(i,2)==0
        h(m-i)=0;
        h(m+i)=0;
    else
        h(m-i)=-1/(i*pi)^2;
        h(m+i)=-1/(i*pi)^2;
    end
end
%% convolution
Rcir = zeros(n,2*m-1);
Rt = proj;
for fi = 1:n
    for r = 1:m
        Rcir(fi,r) = sum(Rt(fi,(1:r)).*h(m-r+1:m));
    end
end
for fi = 1:n
    for r = m+1:2*m-1
        Rcir(fi,r) = sum(Rt(fi,(r-m:m)).*h(1:2*m-r+1));
    end
end
figure('name','convolution');imshow(Rcir, []);title('convolution');
Rcir2 = Rcir(:,1:m);%清除不需要的部分
figure('name','convolution2');imshow(Rcir2, []);title('convolution2');
%% BackProjection
rec = zeros(N, N);
for i = 1: views
    theta = alpha(i);
    pos =  xx .* cos(theta) + yy .* sin(theta);
    temp = interp1(vec, Rcir2(i, :), pos, 'linear');
    rec = rec + temp/N;
    rec(rec<0) = 0;
end
figure; imshow(rec, []); title('Reconstruction');