octave_packages/image-1.0.15/rho_filter.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{filtered} =} rho_filter (@var{proj}, @var{type}, @var{scaling})
  20 ##
  21 ## Filters the parallel ray projections in the columns of @var{proj},
  22 ## according to the filter type chosen by @var{type}.  @var{type}
  23 ## can be chosen from
  24 ## @itemize
  25 ## @item 'none'
  26 ## @item 'Ram-Lak' (default)
  27 ## @item 'Shepp-Logan'
  28 ## @item 'Cosine'
  29 ## @item 'Hann'
  30 ## @item 'Hamming'
  31 ## @end itemize
  32 ##
  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.
  37 ##
  38 ## @defunx {[@var{filtered}, @var{filter}] =} rho_filter (...)
  39 ##
  40 ## This form also returns the frequency response of the filter in
  41 ## the vector @var{filter}.
  42 ##
  43 ## @end defun
  44 ##
  45 ## Performs rho filtering on the parallel ray projections provided.
  46 ##
  47 ## Rho filtering is performed as part of the filtered back-projection
  48 ## method of CT image reconstruction.  It is the filtered part of
  49 ## the name.
  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.
  62 ##
  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.
  66 ## The operations
  67 ## @example
  68 ##     filtered = rho_filter (proj);
  69 ##     reconstruction = iradon (filtered, 1, 'linear', 'none');
  70 ## @end example
  71 ## are exactly equivalent to
  72 ## @example
  73 ##     reconstruction = iradon (proj, 1, 'linear', 'Ram-Lak');
  74 ## @end example
  75 ##
  76 ## Usage example:
  77 ## @example
  78 ##   P = phantom ();
  79 ##   projections = radon (P);
  80 ##   filtered_projections = rho_filter (projections, 'Hamming');
  81 ##   reconstruction = iradon (filtered_projections, 1, 'linear', 'none');
  82 ##   figure, imshow (reconstruction, [])
  83 ## @end example
  84
  85 function [filtered_proj, filt] = rho_filter (proj, type, scaling)
  86
  87   filtered_proj = proj;
  88
  89   if (nargin < 3)
  90     scaling = 1;
  91   endif
  92   if (nargin < 2) || (size (type) == 0)
  93     type = 'ram-lak';
  94   endif
  95
  96   if (strcmpi (type, 'none'))
  97     filt = 1;
  98     return;
  99   endif
 100
 101   if (scaling > 1) || (scaling < 0)
 102     error ('Scaling factor must be in [0,1]');
 103   endif
 104
 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;
 108
 109   ## Scale the frequency response
 110   int_len = (new_len * scaling);
 111   if (mod (floor (int_len), 2))
 112     int_len = ceil (int_len);
 113   else
 114     int_len = floor (int_len);
 115   endif
 116
 117   ## Create the basic filter response
 118   rho = scaling * (0:1 / (int_len / 2):1);
 119   rho = [rho'; rho(end - 1:-1:2)'];
 120
 121   ## Create the window to apply to the filter response
 122   if (strcmpi (type, 'ram-lak'))
 123     filt = 1;
 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);
 134     filt (1) = 1;
 135     filt = [filt; filt(end - 1:-1:2)];
 136   else
 137     error ('rho_filter: Unknown window type');
 138   endif
 139
 140   ## Apply the window
 141   filt = filt .* rho;
 142
 143   ## Pad the response to the correct length
 144   len_diff = new_len - int_len;
 145   if (len_diff != 0)
 146     pad = len_diff / 2;
 147     filt = padarray (fftshift (filt), pad);
 148     filt = fftshift (filt);
 149   endif
 150
 151   filtered_proj = fft (filtered_proj);
 152
 153   ## Perform the filtering
 154   for i = 1:size (filtered_proj, 2)
 155     filtered_proj (:, i) = filtered_proj (:, i) .* filt;
 156   endfor
 157
 158   ## Finally bring the projections back to the spatial domain
 159   filtered_proj = real (ifft (filtered_proj));
 160
 161   ## Chop the projections back to their original size
 162   filtered_proj (size (proj, 1) + 1:end, :) = [];
 163
 164 endfunction
 165
 166 %!demo
 167 %! P = phantom ();
 168 %! projections = radon (P);
 169 %! filtered_projections = rho_filter (projections, 'Hamming');
 170 %! reconstruction = iradon (filtered_projections, 1, 'linear', 'none');
 171 %! figure, imshow (reconstruction, [])
 172
 173 % $Log$
 174 % 2010-05-27 lxop
 175 % Function now returns a value in `filt' when `none' is specified
 176 % as the filter type.  Fixes the warning that was occuring with
 177 % the demo.
 178 %
 179 % 2010-04-08 lxop
 180 % Fixed default argument settings and checking.
 181 %
 182 % 2010-03-19 lxop
 183 % Function completed to Matlab compatible level.