Home > src > main > matlab > s2let_ridgelet_compute_wav.m

s2let_ridgelet_compute_wav

PURPOSE ^

s2let_ridgelet_compute_wav - Compute ridgelets

SYNOPSIS ^

function [ridgelet_wav, ridgelet_scal] =s2let_ridgelet_compute_wav(L, varargin)

DESCRIPTION ^

 s2let_ridgelet_compute_wav - Compute ridgelets

 Compute ridgelet wavelet and scaling functions.

 Default usage :

   [ridgelet_wav, ridgelet_scal] = s2let_compute_wav(L, <options>)

 where L is harmonic band-limit for the reconstruction on the sphere,
 [ridgelet_wav, ridgelet_scal] are the reconstruct ridgelet wavelets and
 scaling functions on the sphere for each j.

 Options consist of parameter type and value pairs.
 Valid options include:

  'B'               = { Dilation factor; B > 1 (default = 2) }
  'Spin'            = { Spin number; Spin >= 0 (default = 0) }
  'J_min'           = { Minimum wavelet scale to consider;
                        0 <= J_min < log_B(L) (default=0) }
  'Sampling'        = { 'MW'           [McEwen & Wiaux sampling (default)],
                        'MWSS'         [McEwen & Wiaux symmetric sampling] }
  'Reality'         = { false        [do not assume f real (default)],
                        true         [assume f real (improves performance)] }

 S2LET package to perform Wavelet transform on the Sphere.
 Copyright (C) 2012-2015-2014  Boris Leistedt, Martin Büttner & Jason McEwen
 See LICENSE.txt for license details

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [ridgelet_wav, ridgelet_scal] = ...
0002    s2let_ridgelet_compute_wav(L, varargin)
0003 % s2let_ridgelet_compute_wav - Compute ridgelets
0004 %
0005 % Compute ridgelet wavelet and scaling functions.
0006 %
0007 % Default usage :
0008 %
0009 %   [ridgelet_wav, ridgelet_scal] = s2let_compute_wav(L, <options>)
0010 %
0011 % where L is harmonic band-limit for the reconstruction on the sphere,
0012 % [ridgelet_wav, ridgelet_scal] are the reconstruct ridgelet wavelets and
0013 % scaling functions on the sphere for each j.
0014 %
0015 % Options consist of parameter type and value pairs.
0016 % Valid options include:
0017 %
0018 %  'B'               = { Dilation factor; B > 1 (default = 2) }
0019 %  'Spin'            = { Spin number; Spin >= 0 (default = 0) }
0020 %  'J_min'           = { Minimum wavelet scale to consider;
0021 %                        0 <= J_min < log_B(L) (default=0) }
0022 %  'Sampling'        = { 'MW'           [McEwen & Wiaux sampling (default)],
0023 %                        'MWSS'         [McEwen & Wiaux symmetric sampling] }
0024 %  'Reality'         = { false        [do not assume f real (default)],
0025 %                        true         [assume f real (improves performance)] }
0026 %
0027 % S2LET package to perform Wavelet transform on the Sphere.
0028 % Copyright (C) 2012-2015-2014  Boris Leistedt, Martin Büttner & Jason McEwen
0029 % See LICENSE.txt for license details
0030 
0031 % Parse arguments.
0032 p = inputParser;
0033 p.addRequired('L', @isnumeric);
0034 p.addParamValue('B', 2, @isnumeric);
0035 p.addParamValue('Spin', 0, @isnumeric);
0036 p.addParamValue('J_min', 0, @isnumeric);
0037 p.addParamValue('Sampling', 'MW', @ischar);
0038 p.addParamValue('Reality', false, @islogical);
0039 p.parse(L, varargin{:});
0040 args = p.Results;
0041 
0042 N = 1;
0043 s = args.Spin; 
0044 
0045 % Compute ring in harmonic space.
0046 ring_lm = zeros(args.L^2,1);
0047 for el = max([0 abs(args.Spin)]):args.L-1
0048         
0049    % Spin zero setting.
0050 %    elon2 = el./2.0;
0051 %    logp = gammaln(2*elon2+1) - 2*elon2 * log(2) - 2 * gammaln(elon2+1);
0052 %    p0 = real((-1).^elon2) .* exp(logp);
0053    
0054    % Spin setting but computed via associated Legendre functions (slow).
0055 %    P = legendre(el, 0);
0056 %    p0 = P(abs(args.Spin)+1);
0057    
0058    % Spin setting but computed explicitly.
0059    
0060    logp2 = gammaln(el+s+1) - el * log(2) - gammaln((el+s)./2+1) - gammaln((el-s)./2+1);
0061    p0 = real((-1).^((el+s)./2)) .* exp(logp2);
0062    
0063    ind = ssht_elm2ind(el, 0);
0064    ring_lm(ind) = 2*pi * sqrt((2*el+1)/(4*pi)) * p0;      
0065    ring_lm(ind) = ring_lm(ind) .* ...   
0066      (-1).^args.Spin .* sqrt(exp(gammaln(el-s+1) - gammaln(el+s+1)));
0067    
0068 end
0069 
0070 
0071 %% Compute ridgets via wavelet transforms
0072 % Can only do this for scalar setting
0073 
0074 % % Compute ring on the sphere.
0075 % ring = ssht_inverse(ring_lm, L, ...
0076 %    'Reality', args.Reality, ...
0077 %    'Method', args.Sampling, ...
0078 %    'Spin', args.Spin);
0079 %
0080 % % Compute ridgelets.
0081 % [ridgelet_wav, ridgelet_scal] = ...
0082 %    s2let_transform_analysis_mw(ring, 'B', args.B, 'J_min', args.J_min, ...
0083 %    'N', N, 'Upsample', true, 'Spin', args.Spin, ...
0084 %    'Reality', args.Reality, 'Sampling', args.Sampling);
0085 
0086 % [ridgelet_wav, ridgelet_scal] = ...
0087 %    s2let_transform_axisym_analysis_mw(ring, 'L', L, 'B', args.B, 'J_min', args.J_min, ...
0088 %    'Upsample', true, ...
0089 %    'Reality', args.Reality);
0090 
0091 
0092 %% Compute ridgelets directly
0093 % Can then compute explicitly for spin functions.
0094 % But not that spin ridgelets are identical to scalar wavelets since can
0095 % always view as convolution of ring with scalar wavelets.
0096 
0097 J = s2let_jmax(args.L, args.B);
0098 [kappa, kappa0] = s2let_transform_axisym_tiling(args.B, args.L, args.J_min);
0099 
0100 ridgelet_scal_lm = zeros(args.L^2,1);
0101 for el = 0:args.L-1
0102    ind = ssht_elm2ind(el, 0);
0103    scal_lm = kappa0(el+1) .* sqrt((2*el+1)/(4*pi));
0104    ridgelet_scal_lm(ind) = scal_lm .* sqrt(4*pi / (2*el+1)) .* ring_lm(ind);
0105 end
0106 ridgelet_scal = ssht_inverse(ridgelet_scal_lm, L, 'Spin', args.Spin, ...
0107    'Reality', args.Reality, 'Method', args.Sampling);
0108 
0109 ridgelet_wav = cell(J-args.J_min+1,1);
0110 for j = args.J_min:J
0111       
0112    ridgelet_wav_lm = zeros(args.L^2,1);
0113    for el = 0:args.L-1
0114       ind = ssht_elm2ind(el, 0);
0115       wav_lm = kappa(j+1, el+1) .* sqrt((2*el+1)/(8*pi.^2)); 
0116       ridgelet_wav_lm(ind) = wav_lm .* sqrt(4*pi / (2*el+1)) .* ring_lm(ind);
0117    end
0118    ridgelet_wav{j-args.J_min+1} = ssht_inverse(ridgelet_wav_lm, L, ...
0119       'Spin', args.Spin, 'Reality', args.Reality, 'Method', args.Sampling);
0120    
0121 end
0122 
0123 
0124 
0125 
0126 
0127

Generated on Fri 11-Nov-2016 11:50:36 by m2html © 2005