--- a/main/image/INDEX	Sun Mar 28 08:55:37 2010 +0000
+++ b/main/image/INDEX	Tue Mar 30 07:37:36 2010 +0000
@@ -56,6 +56,8 @@
  makelut applylut
  deriche
  radon
+ rho_filter
+ iradon
  nonmax_supress
 Black and white image functions
  bwarea

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/image/inst/iradon.m	Tue Mar 30 07:37:36 2010 +0000
@@ -0,0 +1,169 @@
+## Copyright (C) 2010 Alex Opie <lx_op@orcon.net.nz>
+##
+##
+## 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; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @defun @var{recon} = iradon (@var{proj}, @var{theta}, @var{interp}, @
+##                              @var{filter}, @var{scaling}, @var{output_size})
+##
+## Performs filtered back-projection on the projections in @var{proj}
+## to reconstruct an approximation of the original image.
+##
+## @var{proj} should be a matrix whose columns are projections of an
+## image (or slice).  Each element of @var{theta} is used as the angle
+## (in degrees) that the corresponding column of @var{proj} was
+## projected at.  If @var{theta} is omitted, it is assumed that
+## projections were taken at evenly spaced angles between 0 and 180 degrees.
+## @var{theta} can also be a scalar, in which case it is taken as the
+## angle between projections if more than one projection is provided.
+## 
+## @var{interp} determines the type of interpolation that is used
+## in the back-projection.  It must be one of the types accepted by
+## @command{interp1}, and defaults to 'Linear' if it is omitted.
+##
+## @var{filter} and @var{scaling} determine the type of rho filter 
+## to apply.  See the help for @command{rho_filter} for their use.
+##
+## @var{output_size} sets the edge length of the output image (it
+## is always square).  This argument does not scale the image.  If it
+## is omitted, the length is taken to be
+## @group
+## 2 * floor (size (proj, 1) / (2 * sqrt (2))).
+## @end group
+## 
+## If @var{proj} was obtained using @command{radon}, there is no
+## guarantee that the reconstructed image will be exactly the same
+## size as the original.
+## 
+## @defunx [@var{recon}, @var{filt}] = iradon (...)
+##
+## This form also returns the filter frequency response in the vector
+## @var{filt}.
+## @end defun
+##
+## Performs filtered back-projection in order to reconstruct an
+## image based on its projections.
+##
+## Filtered back-projection is the most common means of reconstructing
+## images from CT scans.  It is a two step process: First, each of 
+## the projections is filtered with a `rho filter', so named due
+## to its frequency domain definition, which is simply |rho|, where
+## rho is the radial axis in a polar coordinate system.  Second, 
+## the filtered projections are each `smeared' across the image
+## space.  This is the back-projection part.
+##
+## Usage example:
+## @example
+##   P = phantom ();
+##   projections = radon (P, 1:179);
+##   reconstruction = iradon (filtered_projections, 1:179, 'Spline', 'Hann');
+##   figure, imshow (reconstruction, [])
+## @end example
+
+function [recon, filt] = iradon (proj, theta, interp, filter, scaling, output_size)
+  
+  if (nargin == 0)
+    error ("No projections provided to iradon");
+  endif
+  
+  if (nargin < 6)
+    output_size = 2 * floor (size (proj, 1) / (2 * sqrt (2)));
+  endif
+  if (nargin < 5) || (length (scaling) == 0)
+    scaling = 1;
+  endif
+  if (nargin < 4) || (length (filter) == 0)
+    filter = "Ram-Lak";
+  endif
+  if (nargin < 3) || (length (interp) == 0)
+    interp = "linear";
+  endif
+  if (nargin < 2) || (length (theta) == 0)
+    theta = 180 * (0:1:size (proj, 2) - 1) / size (proj, 2);
+  endif
+  
+  if (isscalar (theta)) && (size (proj, 2) != 1)
+    theta = (0:size (proj, 2) - 1) * theta;
+  endif
+  
+  if (length (theta) != size (proj, 2))
+    error ("iradon: Number of projections does not match number of angles")
+  endif
+  if (!isscalar (scaling))
+    error ("iradon: Frequency scaling value must be a scalar");
+  endif
+  if (!length (find (strcmpi (interp, {'nearest', 'linear', 'spline', \
+                                       'pchip', 'cubic'}))))
+    error ("iradon: Invalid interpolation method specified");
+  endif
+  
+  ## Convert angles to radians
+  theta *= pi / 180;
+  
+  ## First, filter the projections
+  [filtered, filt] = rho_filter (proj, filter, scaling);
+  
+  ## Next, back-project
+  recon = back_project (filtered, theta, interp, output_size);
+  
+endfunction
+
+
+function recon = back_project (proj, theta, interpolation, dim)
+  ## Make an empty image
+  recon = zeros (dim, dim);
+  
+  ## Zero pad the projections if the requested image
+  ## has a diagonal longer than the projections
+  diagonal = ceil (dim * sqrt (2)) + 1;
+  if (size (proj, 1) < diagonal)
+    diff = 2 * ceil ((diagonal - size (proj, 1)) / 2);
+    proj = padarray (proj, diff / 2);
+  endif
+  
+  ## Create the x & y values for each pixel
+  centre = floor ((dim + 1) / 2);
+  x = (0:dim - 1) - centre + 1;
+  x = repmat (x, dim, 1);
+   
+  y = (dim - 1: -1 : 0)' - centre;
+  y = repmat (y, 1, dim);
+  
+  ## s axis for projections, needed by interp1
+  s = (0:size (proj, 1) - 1) - floor (size (proj, 1) / 2);
+  
+  ## Sum each projection's contribution
+  for i = 1:length (theta)
+    s_dash = (x * cos (theta (i)) + y * sin (theta (i)));
+    interpolated = interp1 (s, proj (:, i), s_dash (:), ["*", interpolation]);
+    recon += reshape (interpolated, dim, dim);
+  endfor
+  
+  ## Scale the reconstructed values to their original size
+  recon *= pi / (2 * length (theta));
+  
+endfunction
+
+%!demo
+%! P = phantom ();
+%! figure, imshow (P, []), title ("Original image")
+%! projections = radon (P, 0:179);
+%! reconstruction = iradon (projections, 0:179, 'Spline', 'Hann');
+%! figure, imshow (reconstruction, []), title ("Reconstructed image")
+
+% $Log$
+% 2010-30-03 lxop
+% First submitted to Octave-Forge

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/image/inst/rho_filter.m	Tue Mar 30 07:37:36 2010 +0000
@@ -0,0 +1,175 @@
+## Copyright (C) 2010 Alex Opie <lx_op@orcon.net.nz>
+##
+##
+## 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; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @defun {@var{filtered} =} rho_filter (@var{proj}, @var{type}, @var{scaling})
+##
+## Filters the parallel ray projections in the columns of @var{proj},
+## according to the filter type chosen by @var{type}.  @var{type}
+## can be chosen from
+## @itemize
+## @item 'none'
+## @item 'Ram-Lak' (default)
+## @item 'Shepp-Logan'
+## @item 'Cosine'
+## @item 'Hann'
+## @item 'Hamming'
+## @end itemize
+## 
+## If given, @var{scaling} determines the proportion of frequencies
+## below the nyquist frequency that should be passed by the filter.
+## The window function is compressed accordingly, to avoid an abrupt
+## truncation of the frequency response.
+##
+## @defunx {[@var{filtered}, @var{filter}] =} rho_filter (...)
+##
+## This form also returns the frequency response of the filter in
+## the vector @var{filter}.
+##
+## @end defun
+## 
+## Performs rho filtering on the parallel ray projections provided.
+##
+## Rho filtering is performed as part of the filtered back-projection
+## method of CT image reconstruction.  It is the filtered part of
+## the name.
+## The simplest rho filter is the Ramachadran-Lakshminarayanan (Ram-Lak),
+## which is simply |rho|, where rho is the radial component of spatial
+## frequency.  However, this can cause unwanted amplification of noise, 
+## which is what the other types attempt to minimise, by introducing
+## roll-off into the response.  The Hann and Hamming filters multiply
+## the standard response by a Hann or Hamming window, respectively.
+## The cosine filter is the standard response multiplied by a cosine
+## shape, and the Shepp-Logan filter multiplies the response with
+## a sinc shape.  The 'none' filter performs no filtering, and is
+## included for completeness and to enable incorporating this function
+## easily into scripts or functions that may offer the ability to choose
+## to apply no filtering.
+##
+## This function is designed to be used by the function @command{iradon},
+## but has been exposed to facilitate custom inverse radon transforms
+## and to more clearly break down the process for educational purposes.
+## The operations
+## @example
+##     filtered = rho_filter (proj);
+##     reconstruction = iradon (filtered, 1, 'linear', 'none');
+## @end example
+## are exactly equivalent to
+## @example
+##     reconstruction = iradon (proj, 1, 'linear', 'Ram-Lak');
+## @end example
+##
+## Usage example:
+## @example
+##   P = phantom ();
+##   projections = radon (P);
+##   filtered_projections = rho_filter (projections, 'Hamming');
+##   reconstruction = iradon (filtered_projections, 1, 'linear', 'none');
+##   figure, imshow (reconstruction, [])
+## @end example
+
+function [filtered_proj, filt] = rho_filter (proj, type, scaling)
+
+  filtered_proj = proj;
+
+  if (strcmpi (type, 'none'))
+    return;
+  endif
+  
+  if (nargin < 2) || (size (type) == 0)
+    type = 'ram-lak';
+  endif
+  
+  if (nargin < 3)
+    scaling = 1;
+  endif
+  
+  if (scaling > 1)
+    error ('Scaling factor must be in [0,1]');
+  endif
+  
+  ## Extend the projections to a power of 2
+  new_len = 2 * 2^nextpow2 (size (filtered_proj, 1));
+  filtered_proj (new_len, 1) = 0;
+  
+  ## Scale the frequency response
+  int_len = (new_len * scaling);
+  if (mod (floor (int_len), 2))
+    int_len = ceil (int_len);
+  else
+    int_len = floor (int_len);
+  endif
+  
+  ## Create the basic filter response
+  rho = scaling * (0:1 / (int_len / 2):1);
+  rho = [rho'; rho(end - 1:-1:2)'];
+  
+  ## Create the window to apply to the filter response
+  if (strcmpi (type, 'ram-lak'))
+    filt = 1;
+  elseif (strcmpi (type, 'hamming'))
+    filt = fftshift (hamming (length (rho)));
+  elseif (strcmpi (type, 'hann'))
+    filt = fftshift (hann (length (rho)));
+  elseif (strcmpi (type, 'cosine'))
+    f = 0.5 * (0:length (rho) - 1)' / length (rho);
+    filt = fftshift (sin (2 * pi * f));
+  elseif (strcmpi (type, 'shepp-logan'))
+    f = (0:length (rho) / 2)' / length (rho);
+    filt = sin (pi * f) ./ (pi * f);
+    filt (1) = 1;
+    filt = [filt; filt(end - 1:-1:2)];
+  else
+    error ('rho_filter: Unknown window type');
+  endif
+  
+  ## Apply the window
+  filt = filt .* rho;
+  
+  ## Pad the response to the correct length
+  len_diff = new_len - int_len;
+  if (len_diff != 0)
+    pad = len_diff / 2;
+    filt = padarray (fftshift (filt), pad);
+    filt = fftshift (filt);
+  endif
+  
+  filtered_proj = fft (filtered_proj);
+  
+  ## Perform the filtering
+  for i = 1:size (filtered_proj, 2)
+    filtered_proj (:, i) = filtered_proj (:, i) .* filt;
+  endfor
+  
+  ## Finally bring the projections back to the spatial domain
+  filtered_proj = real (ifft (filtered_proj));
+  
+  ## Chop the projections back to their original size
+  filtered_proj (size (proj, 1) + 1:end, :) = [];
+  
+endfunction
+
+%!demo
+%! P = phantom ();
+%! projections = radon (P);
+%! filtered_projections = rho_filter (projections, 'Hamming');
+%! reconstruction = iradon (filtered_projections, 1, 'linear', 'none');
+%! figure, imshow (reconstruction, [])
+
+% $Log$
+% 2010-03-19 lxop
+% Function completed to Matlab compatible level.