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{recon} = iradon (@var{proj}, @var{theta}, @var{interp}, @
20 ## @var{filter}, @var{scaling}, @var{output_size})
22 ## Performs filtered back-projection on the projections in @var{proj}
23 ## to reconstruct an approximation of the original image.
25 ## @var{proj} should be a matrix whose columns are projections of an
26 ## image (or slice). Each element of @var{theta} is used as the angle
27 ## (in degrees) that the corresponding column of @var{proj} was
28 ## projected at. If @var{theta} is omitted, it is assumed that
29 ## projections were taken at evenly spaced angles between 0 and 180 degrees.
30 ## @var{theta} can also be a scalar, in which case it is taken as the
31 ## angle between projections if more than one projection is provided.
33 ## @var{interp} determines the type of interpolation that is used
34 ## in the back-projection. It must be one of the types accepted by
35 ## @command{interp1}, and defaults to 'Linear' if it is omitted.
37 ## @var{filter} and @var{scaling} determine the type of rho filter
38 ## to apply. See the help for @command{rho_filter} for their use.
40 ## @var{output_size} sets the edge length of the output image (it
41 ## is always square). This argument does not scale the image. If it
42 ## is omitted, the length is taken to be
44 ## 2 * floor (size (proj, 1) / (2 * sqrt (2))).
47 ## If @var{proj} was obtained using @command{radon}, there is no
48 ## guarantee that the reconstructed image will be exactly the same
49 ## size as the original.
51 ## @defunx [@var{recon}, @var{filt}] = iradon (...)
53 ## This form also returns the filter frequency response in the vector
57 ## Performs filtered back-projection in order to reconstruct an
58 ## image based on its projections.
60 ## Filtered back-projection is the most common means of reconstructing
61 ## images from CT scans. It is a two step process: First, each of
62 ## the projections is filtered with a `rho filter', so named due
63 ## to its frequency domain definition, which is simply |rho|, where
64 ## rho is the radial axis in a polar coordinate system. Second,
65 ## the filtered projections are each `smeared' across the image
66 ## space. This is the back-projection part.
71 ## projections = radon (P, 1:179);
72 ## reconstruction = iradon (filtered_projections, 1:179, 'Spline', 'Hann');
73 ## figure, imshow (reconstruction, [])
76 function [recon, filt] = iradon (proj, theta, interp, filter, scaling, output_size)
79 error ("No projections provided to iradon");
83 output_size = 2 * floor (size (proj, 1) / (2 * sqrt (2)));
85 if (nargin < 5) || (length (scaling) == 0)
88 if (nargin < 4) || (length (filter) == 0)
91 if (nargin < 3) || (length (interp) == 0)
94 if (nargin < 2) || (length (theta) == 0)
95 theta = 180 * (0:1:size (proj, 2) - 1) / size (proj, 2);
98 if (isscalar (theta)) && (size (proj, 2) != 1)
99 theta = (0:size (proj, 2) - 1) * theta;
102 if (length (theta) != size (proj, 2))
103 error ("iradon: Number of projections does not match number of angles")
105 if (!isscalar (scaling))
106 error ("iradon: Frequency scaling value must be a scalar");
108 if (!length (find (strcmpi (interp, {'nearest', 'linear', 'spline', \
109 'pchip', 'cubic'}))))
110 error ("iradon: Invalid interpolation method specified");
113 ## Convert angles to radians
116 ## First, filter the projections
117 [filtered, filt] = rho_filter (proj, filter, scaling);
119 ## Next, back-project
120 recon = back_project (filtered, theta, interp, output_size);
125 function recon = back_project (proj, theta, interpolation, dim)
126 ## Make an empty image
127 recon = zeros (dim, dim);
129 ## Zero pad the projections if the requested image
130 ## has a diagonal longer than the projections
131 diagonal = ceil (dim * sqrt (2)) + 1;
132 if (size (proj, 1) < diagonal)
133 diff = 2 * ceil ((diagonal - size (proj, 1)) / 2);
134 proj = padarray (proj, diff / 2);
137 ## Create the x & y values for each pixel
138 centre = floor ((dim + 1) / 2);
139 x = (0:dim - 1) - centre + 1;
140 x = repmat (x, dim, 1);
142 y = (dim - 1: -1 : 0)' - centre;
143 y = repmat (y, 1, dim);
145 ## s axis for projections, needed by interp1
146 s = (0:size (proj, 1) - 1) - floor (size (proj, 1) / 2);
148 ## Sum each projection's contribution
149 for i = 1:length (theta)
150 s_dash = (x * cos (theta (i)) + y * sin (theta (i)));
151 interpolated = interp1 (s, proj (:, i), s_dash (:), ["*", interpolation]);
152 recon += reshape (interpolated, dim, dim);
155 ## Scale the reconstructed values to their original size
156 recon *= pi / (2 * length (theta));
162 %! figure, imshow (P, []), title ("Original image")
163 %! projections = radon (P, 0:179);
164 %! reconstruction = iradon (projections, 0:179, 'Spline', 'Hann');
165 %! figure, imshow (reconstruction, []), title ("Reconstructed image")
169 % First submitted to Octave-Forge