X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Fimage-1.0.15%2Frho_filter.m;fp=octave_packages%2Fimage-1.0.15%2Frho_filter.m;h=0e10bedea689245a354d0f087ebe9d455ca92e12;hp=0000000000000000000000000000000000000000;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hpb=1705066eceaaea976f010f669ce8e972f3734b05

diff --git a/octave_packages/image-1.0.15/rho_filter.m b/octave_packages/image-1.0.15/rho_filter.m
new file mode 100644
index 0000000..0e10bed
--- /dev/null
+++ b/octave_packages/image-1.0.15/rho_filter.m
@@ -0,0 +1,183 @@
+## 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 (nargin < 3)
+    scaling = 1;
+  endif
+  if (nargin < 2) || (size (type) == 0)
+    type = 'ram-lak';
+  endif
+
+  if (strcmpi (type, 'none'))
+    filt = 1;
+    return;
+  endif
+  
+  if (scaling > 1) || (scaling < 0)
+    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-05-27 lxop
+% Function now returns a value in `filt' when `none' is specified
+% as the filter type.  Fixes the warning that was occuring with
+% the demo.
+%
+% 2010-04-08 lxop
+% Fixed default argument settings and checking.
+%
+% 2010-03-19 lxop
+% Function completed to Matlab compatible level.