123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144 |
- %% Author: Rodrigo de Barros Vimieiro
- % Date: April, 2018
- % rodrigo.vimieiro@gmail.com
- % =========================================================================
- %{
- % -------------------------------------------------------------------------
- % filterProj(proj,param,cutoff,filterOrder)
- % -------------------------------------------------------------------------
- % DESCRIPTION:
- % This function filters the projection based on Feldkamp algorithms. It
- % also applies a window filter to discourage noise and avoid ringing
- % effect.
- %
- % The geometry is for DBT with half cone-beam. All parameters are set
- % in "ParameterSettings" code.
- %
- % INPUT:
- %
- % - proj = 2D projection images
- % - param = Parameter of all geometry
- % - cutoff = Cut off frequency up to Nyquist
- % - filterOrder = Butterworth filter order
- %
- % OUTPUT:
- %
- % - filteredProj = Filtered projections.
- %
- % Reference: Jiang Hsieh's book (second edition)
- % Reference: Fessler, J. A.: Fundamentals of CT Reconstruction in 2D
- % and 3D. In: Brahme, A. (eds.) Comprehensive Biomedical Physics,
- % 1st ed., vol. 2, pp 263-295. Elsevier, Netherlands (2014).
- % doi:10.1016/B978-0-444-53632-7.00212-4.
- % Reference: Fessler Book -> (http://web.eecs.umich.edu/~fessler/book/c-tomo.pdf)
- %
- % ---------------------------------------------------------------------
- % Copyright (C) <2018> <Rodrigo de Barros Vimieiro>
- %
- % This program is free software: you can redistribute it and/or modify
- % it under the terms of the GNU General Public License as published by
- % the Free Software Foundation, either version 3 of the License, or
- % (at your option) any later version.
- %
- % This program is distributed in the hope that it will be useful,
- % but WITHOUT ANY WARRANTY; without even the implied warranty of
- % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- % GNU General Public License for more details.
- %
- % You should have received a copy of the GNU General Public License
- % along with this program. If not, see <http://www.gnu.org/licenses/>.
- %}
- % =========================================================================
- %% Filtering Code
- function filteredProj = filterProj(proj,param,cutoff,filterOrder)
- filteredProj = zeros(size(proj),'single');
- % Detector Coordinate sytem in (mm)
- [uCoord,vCoord] = meshgrid(param.us,param.vs);
- % Compute weighted projections (Fessler Book Eq. (3.10.6))
- weightFunction = param.DSO ./ sqrt((param.DSD.^2)+(vCoord.^2)+(uCoord.^2));
-
- % Apply weighting function on each proj
- for i=1:param.nProj
- filteredProj(:,:,i) = proj(:,:,i).* weightFunction;
- end
- % Increase filter length to two times nv to decrease freq step
- h_Length = 2^nextpow2(2*param.nv);
- % Builds ramp filter in space domain
- ramp_kernel = ramp_builder(h_Length);
- % Window filter in freq domain
- H_filter = filter_window(ramp_kernel, h_Length, cutoff, filterOrder);
- % Replicate to all columns to build a 2D filter kernel
- H_filter = repmat(H_filter',1,param.nu);
- % Proj in freq domain
- H_proj = zeros(h_Length,param.nu,'single');
- fprintf('%2d/%2d', 1, param.nProj)
- % Filter each projection
- for i=1:param.nProj
- fprintf('\b\b\b\b\b%2d/%2d', i, param.nProj)
- H_proj(1:param.nv,:) = filteredProj(:,:,i);
- % Fourier transfor in projections
- H_proj = fftshift(fft(H_proj));
-
- % Multiplication in frequency = convolution in space
- H_proj = H_proj .* H_filter;
-
- % Inverse Fourier transfor
- H_proj = (real(ifft(ifftshift(H_proj))));
- filteredProj(:,:,i) = H_proj(1:param.nv,:);
- end
- end
- %% Function Ramp Filter
- %{
- The function builds Ramp Filter in space domain
- Reference: Jiang Hsieh's book (second edition,page 73) Eq. 3.29
- Reference: Fessler Book Eq.(3.4.14)
- %}
- function h = ramp_builder(h_Length)
- n = (-(h_Length/2):(h_Length/2-1));
- h = zeros(size(n),'single');
- h(h_Length/2+1) = 1/4; % Eq. 3.29
- odd = mod(n,2) == 1; % Eq. 3.29
- h(odd) = -1 ./ (pi * n(odd)).^2; % Eq. 3.29
- end
- %% Function Window Ramp Filter
- %{
- The function builds Ramp Filter apodizided in frequency domain
- Reference: Fessler Book and MIRT
- %}
- function H_filter = filter_window(ramp_kernel, h_Length, cutoff, filterOrder)
- H_ramp = abs(fftshift(fft(ramp_kernel))); % Bring filter to freq domain
- w = round(h_Length*cutoff); % Cut off freq
- n = (-(h_Length/2):(h_Length/2-1));
- % Hanning filter + applying cut off frequency
- H_window = 0.5 * (1 + cos(2*pi*n/w)); % Normal Hanning
- H_window = H_window .* (abs(n) < w/2); % Apply cut off freq
- %H_window = 1 ./ (1 + (n/w).^ (2*filterOrder)); %Butterworth
- H_filter = H_ramp .* H_window; % Apply window filter
- end
|