1 ## Copyright (C) 2010 Alex Opie <lx_op@orcon.net.nz>
4 ## This program is free software; you can redistribute it and/or modify it
5 ## under the terms of the GNU General Public License as published by
6 ## the Free Software Foundation; either version 3 of the License, or (at
7 ## your option) any later version.
9 ## This program is distributed in the hope that it will be useful, but
10 ## WITHOUT ANY WARRANTY; without even the implied warranty of
11 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 ## General Public License for more details.
14 ## You should have received a copy of the GNU General Public License
15 ## along with this program; see the file COPYING. If not, see
16 ## <http://www.gnu.org/licenses/>.
19 ## @defun {@var{filtered} =} rho_filter (@var{proj}, @var{type}, @var{scaling})
21 ## Filters the parallel ray projections in the columns of @var{proj},
22 ## according to the filter type chosen by @var{type}. @var{type}
26 ## @item 'Ram-Lak' (default)
27 ## @item 'Shepp-Logan'
33 ## If given, @var{scaling} determines the proportion of frequencies
34 ## below the nyquist frequency that should be passed by the filter.
35 ## The window function is compressed accordingly, to avoid an abrupt
36 ## truncation of the frequency response.
38 ## @defunx {[@var{filtered}, @var{filter}] =} rho_filter (...)
40 ## This form also returns the frequency response of the filter in
41 ## the vector @var{filter}.
45 ## Performs rho filtering on the parallel ray projections provided.
47 ## Rho filtering is performed as part of the filtered back-projection
48 ## method of CT image reconstruction. It is the filtered part of
50 ## The simplest rho filter is the Ramachadran-Lakshminarayanan (Ram-Lak),
51 ## which is simply |rho|, where rho is the radial component of spatial
52 ## frequency. However, this can cause unwanted amplification of noise,
53 ## which is what the other types attempt to minimise, by introducing
54 ## roll-off into the response. The Hann and Hamming filters multiply
55 ## the standard response by a Hann or Hamming window, respectively.
56 ## The cosine filter is the standard response multiplied by a cosine
57 ## shape, and the Shepp-Logan filter multiplies the response with
58 ## a sinc shape. The 'none' filter performs no filtering, and is
59 ## included for completeness and to enable incorporating this function
60 ## easily into scripts or functions that may offer the ability to choose
61 ## to apply no filtering.
63 ## This function is designed to be used by the function @command{iradon},
64 ## but has been exposed to facilitate custom inverse radon transforms
65 ## and to more clearly break down the process for educational purposes.
68 ## filtered = rho_filter (proj);
69 ## reconstruction = iradon (filtered, 1, 'linear', 'none');
71 ## are exactly equivalent to
73 ## reconstruction = iradon (proj, 1, 'linear', 'Ram-Lak');
79 ## projections = radon (P);
80 ## filtered_projections = rho_filter (projections, 'Hamming');
81 ## reconstruction = iradon (filtered_projections, 1, 'linear', 'none');
82 ## figure, imshow (reconstruction, [])
85 function [filtered_proj, filt] = rho_filter (proj, type, scaling)
92 if (nargin < 2) || (size (type) == 0)
96 if (strcmpi (type, 'none'))
101 if (scaling > 1) || (scaling < 0)
102 error ('Scaling factor must be in [0,1]');
105 ## Extend the projections to a power of 2
106 new_len = 2 * 2^nextpow2 (size (filtered_proj, 1));
107 filtered_proj (new_len, 1) = 0;
109 ## Scale the frequency response
110 int_len = (new_len * scaling);
111 if (mod (floor (int_len), 2))
112 int_len = ceil (int_len);
114 int_len = floor (int_len);
117 ## Create the basic filter response
118 rho = scaling * (0:1 / (int_len / 2):1);
119 rho = [rho'; rho(end - 1:-1:2)'];
121 ## Create the window to apply to the filter response
122 if (strcmpi (type, 'ram-lak'))
124 elseif (strcmpi (type, 'hamming'))
125 filt = fftshift (hamming (length (rho)));
126 elseif (strcmpi (type, 'hann'))
127 filt = fftshift (hann (length (rho)));
128 elseif (strcmpi (type, 'cosine'))
129 f = 0.5 * (0:length (rho) - 1)' / length (rho);
130 filt = fftshift (sin (2 * pi * f));
131 elseif (strcmpi (type, 'shepp-logan'))
132 f = (0:length (rho) / 2)' / length (rho);
133 filt = sin (pi * f) ./ (pi * f);
135 filt = [filt; filt(end - 1:-1:2)];
137 error ('rho_filter: Unknown window type');
143 ## Pad the response to the correct length
144 len_diff = new_len - int_len;
147 filt = padarray (fftshift (filt), pad);
148 filt = fftshift (filt);
151 filtered_proj = fft (filtered_proj);
153 ## Perform the filtering
154 for i = 1:size (filtered_proj, 2)
155 filtered_proj (:, i) = filtered_proj (:, i) .* filt;
158 ## Finally bring the projections back to the spatial domain
159 filtered_proj = real (ifft (filtered_proj));
161 ## Chop the projections back to their original size
162 filtered_proj (size (proj, 1) + 1:end, :) = [];
168 %! projections = radon (P);
169 %! filtered_projections = rho_filter (projections, 'Hamming');
170 %! reconstruction = iradon (filtered_projections, 1, 'linear', 'none');
171 %! figure, imshow (reconstruction, [])
175 % Function now returns a value in `filt' when `none' is specified
176 % as the filter type. Fixes the warning that was occuring with
180 % Fixed default argument settings and checking.
183 % Function completed to Matlab compatible level.