CTrecon-an improved parallel data image reconstruct function

Reconstruct an image from parallel projection data using the convolution-backprojection algorithm
function [img,H] = CTrecon(prj,varargin)

for a discussion of this function including examples see: http://aprendtech.com/wordpress/?p=74

Mathworks has created a workaround to fix the offset problem in iradon for R2010 and later users:
http://www.mathworks.com/support/bugreports/744762
CTrecon will work with earlier versions and also has the other features described here.

Has the following improvements over iradon:

1. Fixes a bug in iradon that results in offsets from
the true values in the reconstructed image. Does this by computing
the filter transfer function as transform of direct space function
instead of transforming a ramp in frequency.
See Ch. 3 of "Principles of Computerized Tomographic Imaging" by
A. C. Kak and Malcolm Slaney, IEEE Press, 1988.
available online at http://www.slaney.org/pct/pct-toc.html

2. Uses mex functions based on iradon_speedy to back project so it is faster

3. handles complex projections to return a complex reconstructed image - see zBackproject.c

4. regularizes the interface to ('parameter_name', parameter_value) syntax

inputs:
prj: (required) the projection data matrix of size (number of lines)X(number_of_angles)
angles: optional = angles are equally spaced from 0 to 180 degrees but do not wrap around
number of angles is size(prj,2)

interpolation: optional = 'linear', a string specifying the interpolation in the backprojection
'linear' - linear interpolation (default)
'nearest neighbor' - nearest neighbor interpolation
filter_type:(optional = 'Ram-Lak') a string that specifies the filter used. Possibilities are
'Ram-Lak' The cropped Ramachandran-Lakshimrayanan or ramp filter_type (default).
'Shepp-Logan' The Shepp-Logan filter_type: a ramp times a sinc
'Cosine' a ramp times cosine from 0 to the Nyquist frequency.
'Hamming' a ramp times the Hamming filter_type.
'Hann' a ramp times the Hann filter_type.

freq_cutoff: (optional = 1) sets the filter to 0 starting at freq_cutoff*freq_nyquist

outputs:
img: the reconstructed image
H: the transfer function of the filter

The image output size is determined from the length of the projections:

OUTPUT_SIZE = 2*floor(size(prj,1)/(2*sqrt(2)))

Example: see http://aprendtech.com/wordpress/?p=74 for other examples
-------
P = phantom(128);
R = radon(P,0:179);
img = CTrecon(R,'filter_type','cosine','freq_cutoff',0.9);
figure, imshow(img,[-.05 .3])