How can I obtain filtered back projection using 1D Fourier Transform
조회 수: 8 (최근 30일)
이전 댓글 표시
I have an assignment which consist of two part. First part is forward projection of an image and the second part is filtered back projection from the projection data. I managed the first part. In the second part, I can obtain the reconstructed image without filter; however, I cannot apply the filter. The built-in functions such as iradon is forbidden. I am storing the projection data in a matrix MAT, which is MXN. M = length(angles); N = length(beam_number); I have tried the following:
% Triangular Filter
triang_filter = zeros(1,length(tValues));
for ind8 = 1:1:((length(tValues)+1)/2) %positive slope part of the triangular wave
triang_filter(ind8) = 2 * ((ind8 - 1) / (length(tValues) - 1));
end
for ind9 = ((length(tValues)+1)/2):1:(length(tValues)) %negative slope part of the triangular wave
triang_filter(ind9) = 2 * ((length(tValues) - ind9) / (length(tValues) - 1));
end
ft_p_theta_of_t = fft(p_theta_of_t');
ft_p_theta_of_t = ft_p_theta_of_t';
filtered_p_theta_of_t = [];
for ind11 = 1:1:length(thetaValues) %multiply the fft of projection with filter
filtered_p_theta_of_t(ind11,:) = triang_filter .* ft_p_theta_of_t(ind11,:);
end
inv_filtered_p_theta_of_t = ifft(filtered_p_theta_of_t');
inv_filtered_p_theta_of_t = inv_filtered_p_theta_of_t';
Then, I am using back projection algorithm; however, unfiltered image has less blurring than filtered image.
Thank you in advance :)
댓글 수: 0
답변 (2개)
Ross Hanna
2020년 6월 15일
For future reference for people seeing this, two suggestions.
For the triangle filter, either use the triang() function;
or in one line
filt = (-abs(linspace(-1,1,pixNum))+1); % Creates a triangle (Ram-Lak) filter
% tmp = linspace(-1,1,pixNum) % creates a single line from -1 to +1 intercepting x-axis at (0,0) with length "pixNum"
% tmp = abs(tmp) % finds the absolute (all positive) of the previous tmp variable (triangle, with point at bottom)
% tmp = -tmp % translates filter around y = 0, (triangle with point at top, but in negative domain)
% filt = tmp + 1 % shifts filter up by 1
this creates the Ram-Lak (triangle filter) in one line as opposed to two loops.
For applying the filter, I'm reconstructing using *.tif images for the x-ray projections, hence the variable name "tf".
*.tif images are imported to matlab, converted to double() and stored in "A" (300 x 370 Matrix)
pixNum = 300; % number of pixels in reconstruction
Nrot = 370; % number of rotation steps (360 degrees rotation in 370 steps)
filt = (-abs(linspace(-1,1,pixNum))+1);
for i = 1:Nrot
tf = A(:,rot); % pick the 1D values to backproject
tf = fft(tf); % move into frequency domain
tf = tf .* filt; % multiply frequency domain values by the Ram-Lak filter. (multiplication in frequency domain is the same as convolution in spatial domain)
tf = ifft(tf); % bring back into spatial domain
tf = real(tf); % pick the real part of the signal (this could be debated, something I'm currently investigating)
Output = backProject(tf); %Not an actual function, just an example that you then back project the filtered data for each rotation step
end
Hopefully this helps
EDIT - just tried running your example. it throws an error if using an even number for the tValues variable. as it then gives non-integer indices.
댓글 수: 0
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!