s2let_radon_transform

 s2let_radon_transform
 Compute radon transform in harmonic space.

 Default usage:

   f_radon_lm = s2let_radon_transform(f_lm, <options>)

 where f_lm is the vector of L^2 harmonic coefficients and f_random_lm
 is the vector of harmonic coefficients of the spherical radon transform
 of f.

 Option :
  'Reality'         = { false        [do not assume f real (default)],
                        true         [assume f real (improves performance)] }
  'Spin'            = { Spin; (default=0) }
  'L'               = { Harmonic band-limit; L > 1 (default=guessed from input) }

0001 function [f_radon_lm] = s2let_radon_transform(f_lm, varargin)
0002 
0003 % s2let_radon_transform
0004 % Compute radon transform in harmonic space.
0005 %
0006 % Default usage:
0007 %
0008 %   f_radon_lm = s2let_radon_transform(f_lm, <options>)
0009 %
0010 % where f_lm is the vector of L^2 harmonic coefficients and f_random_lm
0011 % is the vector of harmonic coefficients of the spherical radon transform
0012 % of f.
0013 %
0014 % Option :
0015 %  'Reality'         = { false        [do not assume f real (default)],
0016 %                        true         [assume f real (improves performance)] }
0017 %  'Spin'            = { Spin; (default=0) }
0018 %  'L'               = { Harmonic band-limit; L > 1 (default=guessed from input) }
0019 
0020 
0021 % S2LET package to perform Wavelets transform on the Sphere.
0022 % Copyright (C) 2015  Boris Leistedt & Jason McEwen
0023 % See LICENSE.txt for license details
0024 
0025 L_guess = sqrt(length(f_lm));
0026 
0027 p = inputParser;
0028 p.addRequired('f_lm', @isnumeric);
0029 p.addParamValue('L', L_guess, @isnumeric);
0030 p.addParamValue('Reality', false, @islogical);
0031 p.addParamValue('Spin', 0, @isnumeric);
0032 p.parse(f_lm, varargin{:});
0033 args = p.Results;
0034 
0035 s = args.Spin; 
0036 
0037 % Compute radon transform
0038 ring_lm = zeros(args.L^2,1);
0039 f_radon_lm = zeros(args.L^2,1);
0040 for el = max([0 abs(args.Spin)]):args.L-1
0041                 
0042    logp2 = gammaln(el+s+1) - el * log(2) - gammaln((el+s)./2+1) - gammaln((el-s)./2+1);
0043    p0 = real((-1).^((el+s)./2)) .* exp(logp2);    
0044    ind = ssht_elm2ind(el, 0);
0045    ring_lm(ind) = 2*pi * sqrt((2*el+1)/(4*pi)) * p0;
0046    ring_lm(ind) = ring_lm(ind) .* ...
0047       (-1).^s .* sqrt(exp(gammaln(el-s+1) - gammaln(el+s+1)));
0048    
0049    if args.Reality
0050       m_min = 0;
0051    else
0052       m_min = -el;
0053    end
0054 
0055    for m = m_min:el
0056       ind_lm = ssht_elm2ind(el, m);
0057       f_radon_lm(ind_lm) = sqrt(4 * pi ./ (2*el+1)) ...
0058          * f_lm(ind_lm) * ring_lm(ind);
0059    end
0060        
0061 end
0062