octave_packages/image-1.0.15/iradon.m

   1 ## Copyright (C) 2010 Alex Opie <lx_op@orcon.net.nz>
   2 ##
   3 ##
   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.
   8 ##
   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.
  13 ##
  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/>.
  17
  18 ## -*- texinfo -*-
  19 ## @defun @var{recon} = iradon (@var{proj}, @var{theta}, @var{interp}, @
  20 ##                              @var{filter}, @var{scaling}, @var{output_size})
  21 ##
  22 ## Performs filtered back-projection on the projections in @var{proj}
  23 ## to reconstruct an approximation of the original image.
  24 ##
  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.
  32 ##
  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.
  36 ##
  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.
  39 ##
  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
  43 ## @group
  44 ## 2 * floor (size (proj, 1) / (2 * sqrt (2))).
  45 ## @end group
  46 ##
  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.
  50 ##
  51 ## @defunx [@var{recon}, @var{filt}] = iradon (...)
  52 ##
  53 ## This form also returns the filter frequency response in the vector
  54 ## @var{filt}.
  55 ## @end defun
  56 ##
  57 ## Performs filtered back-projection in order to reconstruct an
  58 ## image based on its projections.
  59 ##
  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.
  67 ##
  68 ## Usage example:
  69 ## @example
  70 ##   P = phantom ();
  71 ##   projections = radon (P, 1:179);
  72 ##   reconstruction = iradon (filtered_projections, 1:179, 'Spline', 'Hann');
  73 ##   figure, imshow (reconstruction, [])
  74 ## @end example
  75
  76 function [recon, filt] = iradon (proj, theta, interp, filter, scaling, output_size)
  77
  78   if (nargin == 0)
  79     error ("No projections provided to iradon");
  80   endif
  81
  82   if (nargin < 6)
  83     output_size = 2 * floor (size (proj, 1) / (2 * sqrt (2)));
  84   endif
  85   if (nargin < 5) || (length (scaling) == 0)
  86     scaling = 1;
  87   endif
  88   if (nargin < 4) || (length (filter) == 0)
  89     filter = "Ram-Lak";
  90   endif
  91   if (nargin < 3) || (length (interp) == 0)
  92     interp = "linear";
  93   endif
  94   if (nargin < 2) || (length (theta) == 0)
  95     theta = 180 * (0:1:size (proj, 2) - 1) / size (proj, 2);
  96   endif
  97
  98   if (isscalar (theta)) && (size (proj, 2) != 1)
  99     theta = (0:size (proj, 2) - 1) * theta;
 100   endif
 101
 102   if (length (theta) != size (proj, 2))
 103     error ("iradon: Number of projections does not match number of angles")
 104   endif
 105   if (!isscalar (scaling))
 106     error ("iradon: Frequency scaling value must be a scalar");
 107   endif
 108   if (!length (find (strcmpi (interp, {'nearest', 'linear', 'spline', \
 109                                        'pchip', 'cubic'}))))
 110     error ("iradon: Invalid interpolation method specified");
 111   endif
 112
 113   ## Convert angles to radians
 114   theta *= pi / 180;
 115
 116   ## First, filter the projections
 117   [filtered, filt] = rho_filter (proj, filter, scaling);
 118
 119   ## Next, back-project
 120   recon = back_project (filtered, theta, interp, output_size);
 121
 122 endfunction
 123
 124
 125 function recon = back_project (proj, theta, interpolation, dim)
 126   ## Make an empty image
 127   recon = zeros (dim, dim);
 128
 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);
 135   endif
 136
 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);
 141
 142   y = (dim - 1: -1 : 0)' - centre;
 143   y = repmat (y, 1, dim);
 144
 145   ## s axis for projections, needed by interp1
 146   s = (0:size (proj, 1) - 1) - floor (size (proj, 1) / 2);
 147
 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);
 153   endfor
 154
 155   ## Scale the reconstructed values to their original size
 156   recon *= pi / (2 * length (theta));
 157
 158 endfunction
 159
 160 %!demo
 161 %! P = phantom ();
 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")
 166
 167 % $Log$
 168 % 2010-30-03 lxop
 169 % First submitted to Octave-Forge