2020-1-24
dudoNet中的反投影
表达式
在dudo-Net这篇文章中,用到了对反投影的一个表达,这个表达令我眼前一亮。
是这样的:
其中,
整个过程实际上就是滤波反投影(FBP)的过程。
通过这个表达,我们可以得到CT的图像域与sino图域之间的梯度的一个关系:
从此,我们得到了一个关于两个神经网络:sinogram domain enhancement net(SE-NET)和image domain enhancement net(IE-NET)的一个关系
matlab实现反投影
我们通过该表达式来实现反投影
整个过程是这样的:
1.建立图形
2.投影
3.傅里叶变换,滤波,傅里叶逆变换
4.反投影
%% RadonInversionLayer
% 2020-1-24
close all; clear; clc;
%% 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;
%% Generate model
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');
%% Generate projection
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');
%% Filter
% filter = [0: 2/N: 1-2/N, 1: -2/N: 2/N];
filter = 2*[0: (N/2-1), N/2: -1: 1]/N;
proj_filter = zeros(views, bins);
for i = 1: views
proj_filter(i, :) = fft(proj(i, :));
proj_filter(i, :) = proj_filter(i, :) .* filter;
proj_filter(i, :) = ifft(proj_filter(i, :));
end
figure; imshow(proj_filter, []); title('Filter Projection image');
%%BackProjection2
rec = zeros(N, N);
for i = 1: views
theta = alpha(i);
pos = xx .* cos(theta) + yy .* sin(theta);
temp1 = interp1(vec, proj_filter(i, :), floor(pos), 'linear');
temp2 = interp1(vec, proj_filter(i, :), ceil(pos), 'linear');
temp = (ceil(pos)-pos).*temp1+(pos-floor(pos)).*temp2;
rec = rec + temp/N;
rec(rec<0) = 0;
end
figure; imshow(rec, []); title('Reconstruction2');
在这里,投影和反投影过程中运用了插值,第三步我运用了快速傅里叶变换(fft)
效果如下