flaglet_axisym_analysis

 flaglet_axisym_analysis 
 Compute axisymmetric wavelet transform, output in pixel space.

 Default usage :

   [f_wav, f_scal] = flaglet_axisym_analysis(f, <options>)

 f is the input field -- MW sampling,
 f_wav contains the output wavelet contributions,
 f_scal contains the output scaling contributions,
 B_l is the wavelet parameter for angular space,
 B_p is the wavelet parameter for radial space,
 L is the angular band-limit,
 P is the radial band-limit,
 J_min_l the first angular wavelet scale to use,
 J_min_p the first radial wavelet scale to use,
 R is the radial boundary-limit.

 Options :
  'B_l'               = { Dilation factor; B_l > 1 (default=2) }
  'B_p'               = { Dilation factor; B_p > 1 (default=2) }
  'L'               = { Angular harmonic band-limit; L > 1 (default=guessed) }
  'P'               = { Radial harmonic band-limit; P > 1 (default=guessed) }
  'J_min_l'           = { Minimum needlet scale to consider;
                        0 <= J_min_l < log_B_l(L) (default=0) }
  'J_min_p'           = { Minimum needlet scale to consider;
                        0 <= J_min_p < log_B_p(P) (default=0) }
  'R'               = { Radial boundary; R > 0 (default=1.0) }
  'Downsample'      = { true        [multiresolution algorithm (default)],
                        false       [full resolution wavelets] }
  'Reality'         = { false        [do not assume f real (default)],
                        true         [assume f real (improves performance)] }

 flaglet package to perform Wavelets transform on the Solid Sphere.
 Copyright (C) 2012  Boris Leistedt & Jason McEwen
 See LICENSE.txt for license details

0001 function [f_wav, f_scal] = flaglet_axisym_analysis(f, varargin)
0002 
0003 % flaglet_axisym_analysis
0004 % Compute axisymmetric wavelet transform, output in pixel space.
0005 %
0006 % Default usage :
0007 %
0008 %   [f_wav, f_scal] = flaglet_axisym_analysis(f, <options>)
0009 %
0010 % f is the input field -- MW sampling,
0011 % f_wav contains the output wavelet contributions,
0012 % f_scal contains the output scaling contributions,
0013 % B_l is the wavelet parameter for angular space,
0014 % B_p is the wavelet parameter for radial space,
0015 % L is the angular band-limit,
0016 % P is the radial band-limit,
0017 % J_min_l the first angular wavelet scale to use,
0018 % J_min_p the first radial wavelet scale to use,
0019 % R is the radial boundary-limit.
0020 %
0021 % Options :
0022 %  'B_l'               = { Dilation factor; B_l > 1 (default=2) }
0023 %  'B_p'               = { Dilation factor; B_p > 1 (default=2) }
0024 %  'L'               = { Angular harmonic band-limit; L > 1 (default=guessed) }
0025 %  'P'               = { Radial harmonic band-limit; P > 1 (default=guessed) }
0026 %  'J_min_l'           = { Minimum needlet scale to consider;
0027 %                        0 <= J_min_l < log_B_l(L) (default=0) }
0028 %  'J_min_p'           = { Minimum needlet scale to consider;
0029 %                        0 <= J_min_p < log_B_p(P) (default=0) }
0030 %  'R'               = { Radial boundary; R > 0 (default=1.0) }
0031 %  'Downsample'      = { true        [multiresolution algorithm (default)],
0032 %                        false       [full resolution wavelets] }
0033 %  'Reality'         = { false        [do not assume f real (default)],
0034 %                        true         [assume f real (improves performance)] }
0035 %
0036 % flaglet package to perform Wavelets transform on the Solid Sphere.
0037 % Copyright (C) 2012  Boris Leistedt & Jason McEwen
0038 % See LICENSE.txt for license details
0039 
0040 sz = size(f);
0041 Pguessed = sz(1);
0042 Lguessed = sz(2);
0043 
0044 p = inputParser;
0045 p.addRequired('f', @isnumeric); 
0046 p.addParamValue('B_l', 2, @isnumeric);       
0047 p.addParamValue('B_p', 2, @isnumeric);          
0048 p.addParamValue('L', Lguessed, @isnumeric);   
0049 p.addParamValue('P', Pguessed, @isnumeric);    
0050 p.addParamValue('J_min_l', 0, @isnumeric);   
0051 p.addParamValue('J_min_p', 0, @isnumeric); 
0052 p.addParamValue('R', 1.0, @isnumeric); 
0053 p.addParamValue('Downsample', true, @islogical);
0054 p.addParamValue('Reality', false, @islogical);
0055 p.parse(f, varargin{:});
0056 args = p.Results;
0057 
0058 P = args.P;
0059 L = args.L;
0060 J_min_l = args.J_min_l;
0061 J_min_p = args.J_min_p;
0062 f_vec = zeros(P, L*(2*L-1));
0063 for p = 1:args.P
0064     temp(:,:) = f(p,:,:);
0065     f_vec(p,:) = flag_mw_arr2vec( temp );
0066 end
0067 
0068 [f_wav_vec, f_scal_vec] = flaglet_axisym_analysis_mex(f_vec, args.B_l, args.B_p, args.L, args.P, args.J_min_l, args.J_min_p, args.R, args.Reality, args.Downsample);
0069 
0070 J_l = ceil(log(L) ./ log(args.B_l));
0071 J_p = ceil(log(P) ./ log(args.B_p));
0072 f_wav = cell(J_l+1-J_min_l, J_p+1-J_min_p);
0073 
0074 offset = 0;
0075 for jp = J_min_p:J_p
0076     if args.Downsample == true
0077         band_limit_p = min([ ceil(args.B_p^(jp+1)) P ]);
0078      else
0079         band_limit_p = P;
0080     end
0081     for jl = J_min_l:J_l
0082         if args.Downsample == true
0083             band_limit_l = min([ ceil(args.B_l^(jl+1)) L ]);
0084         else
0085             band_limit_l = L;
0086         end
0087         temp = zeros(band_limit_p, band_limit_l, 2*band_limit_l-1);
0088         for r = 0:band_limit_p-1
0089             for t = 0:band_limit_l-1
0090                 for p = 0:2*band_limit_l-2
0091                     ind = offset + r * band_limit_l *( 2 * band_limit_l - 1 ) + t * ( 2 * band_limit_l - 1) + p + 1;
0092                     temp(r+1,t+1,p+1) = f_wav_vec(1,ind);
0093                 end
0094             end
0095         end
0096         f_wav{jl+1-J_min_l, jp+1-J_min_p} = temp;
0097         offset = offset + band_limit_l * (2*band_limit_l-1) * band_limit_p;
0098     end
0099 end
0100 
0101 f_scal = zeros(P, L, (2*L-1));
0102 for p = 1:P
0103     temp = f_scal_vec(p,:);
0104     size(flag_mw_vec2arr( temp ));
0105     f_scal(p,:,:) = flag_mw_vec2arr( temp );
0106 end
0107 
0108 end