/usr/share/octave/packages/image-2.2.2/rho_filter.m is in octave-image 2.2.2-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 | ## 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; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {Function File} {@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.
##
## @end deftypefn
## @deftypefn {Function File} {[@var{filtered}, @var{filter}] =} rho_filter (@dots{})
##
## This form also returns the frequency response of the filter in
## the vector @var{filter}.
##
## 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
##
## @end deftypefn
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 (hanning (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, [])
|