s2let_transform_curvelet_analysis_lm2cur

 s2let_transform_curvelet_analysis_lm2cur
 Compute curvelet transform:
 input in harmonic space  (i.e. harmonic to Wigner space via analysis_lm2lmn),
 output in curvelet space (i.e. Wigner space to curvelet space via SO3_inverse_direct).

 Default usage :

   [f_cur, f_scal] = s2let_transform_curvelet_analysis_lm2cur(flm_init, <options>)

 flm_init is the input field in harmonic space,
 f_cur contains the output curvelet contributions,
 f_scal contains the output scaling contributions,

 Option :
  'B'               = { Dilation factor; B > 1 (default=2) }
  'L'               = { Harmonic band-limit; L > 1 (default=guessed from input) }
  'J_min'           = { Minimum curvelet scale to consider;
  'Spin'               = { Spin; (default=0) }
                        0 <= J_min < log_B(L) (default=0) }
  'Reality'         = { false        [do not assume corresponding signal f real (default)],
                        true         [assume f real (improves performance)] }
  'Upsample'        = { false        [multiresolution algorithm (default)],
                        true       [full resolution curvelets] }
  'SpinLowered'     = { true  [Apply normalisation factors for spin-lowered
                               curvelets and scaling function.],
                        false [Apply the usual normalisation factors such
                               that the curvelets fulfil the admissibility
                               condition (default)]}
  'SpinLoweredFrom' = [integer; if the SpinLowered option is used, this
                       option indicates which spin number the curvelets
                       should be lowered from (default = 0)]

 -----------------------------------------------------------
 S2LET package to perform Wavelet Transform on the Sphere.
 Copyright (C) 2012-2016  Boris Leistedt, Jennifer Chan & Jason McEwen
 See LICENSE.txt for license details
 -----------------------------------------------------------

0001 function [f_cur, f_scal] = s2let_transform_curvelet_analysis_lm2cur(flm_init, varargin)
0002 
0003 % s2let_transform_curvelet_analysis_lm2cur
0004 % Compute curvelet transform:
0005 % input in harmonic space  (i.e. harmonic to Wigner space via analysis_lm2lmn),
0006 % output in curvelet space (i.e. Wigner space to curvelet space via SO3_inverse_direct).
0007 %
0008 % Default usage :
0009 %
0010 %   [f_cur, f_scal] = s2let_transform_curvelet_analysis_lm2cur(flm_init, <options>)
0011 %
0012 % flm_init is the input field in harmonic space,
0013 % f_cur contains the output curvelet contributions,
0014 % f_scal contains the output scaling contributions,
0015 %
0016 % Option :
0017 %  'B'               = { Dilation factor; B > 1 (default=2) }
0018 %  'L'               = { Harmonic band-limit; L > 1 (default=guessed from input) }
0019 %  'J_min'           = { Minimum curvelet scale to consider;
0020 %  'Spin'               = { Spin; (default=0) }
0021 %                        0 <= J_min < log_B(L) (default=0) }
0022 %  'Reality'         = { false        [do not assume corresponding signal f real (default)],
0023 %                        true         [assume f real (improves performance)] }
0024 %  'Upsample'        = { false        [multiresolution algorithm (default)],
0025 %                        true       [full resolution curvelets] }
0026 %  'SpinLowered'     = { true  [Apply normalisation factors for spin-lowered
0027 %                               curvelets and scaling function.],
0028 %                        false [Apply the usual normalisation factors such
0029 %                               that the curvelets fulfil the admissibility
0030 %                               condition (default)]}
0031 %  'SpinLoweredFrom' = [integer; if the SpinLowered option is used, this
0032 %                       option indicates which spin number the curvelets
0033 %                       should be lowered from (default = 0)]
0034 %
0035 % -----------------------------------------------------------
0036 % S2LET package to perform Wavelet Transform on the Sphere.
0037 % Copyright (C) 2012-2016  Boris Leistedt, Jennifer Chan & Jason McEwen
0038 % See LICENSE.txt for license details
0039 % -----------------------------------------------------------
0040 
0041 sz = length(flm_init(:));
0042 Lguessed = sqrt(sz);
0043 
0044 p = inputParser;
0045 p.addRequired('flm_init', @isnumeric);
0046 p.addParamValue('B', 2, @isnumeric);
0047 p.addParamValue('L', Lguessed, @isnumeric);
0048 p.addParamValue('J_min', 0, @isnumeric);
0049 p.addParamValue('Spin', 0, @isnumeric);
0050 p.addParamValue('Reality', false, @islogical);
0051 p.addParamValue('Upsample', false, @islogical);
0052 p.addParamValue('SpinLowered', false, @islogical);
0053 p.addParamValue('SpinLoweredFrom', 0, @isnumeric);
0054 p.addParamValue('Sampling', 'MW', @ischar);
0055 p.parse(flm_init, varargin{:});
0056 args = p.Results;
0057 
0058 J = s2let_jmax(args.L, args.B);
0059 
0060 % ---------------
0061 % Tile curvelets:
0062 % ---------------
0063 % Call curvelet- and scaling-function- generating functions
0064 [cur_lm scal_l] = s2let_curvelet_tiling(args.B, args.L, args.J_min,...
0065                                         'Spin', args.Spin,...
0066                                         'SpinLowered', args.SpinLowered, ...
0067                                         'SpinLoweredFrom', args.SpinLoweredFrom);
0068 
0069 % -----------------
0070 % Signal analysis:
0071 % -----------------
0072 % Decompose the signals using curvelets and the scaling functions
0073 % Then perform Wigner transform (calling matlab function s2let_transform_analysis_lm2lmn)
0074 [f_cur_lmn, f_scal_lm] = s2let_transform_curvelet_analysis_lm2lmn(flm_init, cur_lm, scal_l,...
0075                                                                   'B',args.B, 'L', args.L, 'J_min', args.J_min, ...
0076                                                                   'Spin', args.Spin,'Reality', args.Reality,...
0077                                                                   'Upsample', args.Upsample, ...
0078                                                                   'SpinLowered', args.SpinLowered, ...
0079                                                                   'SpinLoweredFrom',  args.SpinLoweredFrom, ...
0080                                                                   'Sampling', args.Sampling);
0081 
0082 % Curvelet contribution:
0083 % Rotate the Wigner coefficients f_cur_lmn (such the curvelets centered at the North pole)
0084 % Exploit the property of curvelets that cur_ln = (cur_ll)*(delta_ln)
0085 % such that cur_lmk_rotated = cur_lml*conj(Dlkl(evaluated at the desired rotation angle))
0086 % ---------------
0087 % Define Euler angles for rotation:
0088 % ---------------
0089 alpha = 0;
0090 gamma = 0 ;
0091 
0092 for j = args.J_min:J,
0093     band_limit = min([ s2let_bandlimit(j,args.J_min,args.B,args.L) args.L ]);
0094     % Nj = orientational band-limit at j-th scale:
0095     Nj = band_limit;
0096 % ---------------
0097 % Specify beta and orecompute Wigner small-d functions, denoted here as d (in the paper: d_lmn for all el, m, n evaluated at beta).
0098 % They are indexed d(el,m,n). Alpha and gamma are the other two rotation angles.
0099 % ---------------
0100     if (args.Upsample ~= 0)
0101         beta = acos(-args.Spin/args.B^j);
0102         d = zeros(args.L, 2*args.L-1, 2*args.L-1);
0103         d(1,:,:) = ssht_dl(squeeze(d(1,:,:)), args.L, 0, beta);
0104         for el = 1:args.L-1
0105             d(el+1,:,:) = ssht_dl(squeeze(d(el,:,:)), args.L, el, beta);
0106         end
0107     else
0108         beta = acos(-args.Spin/args.B^j);
0109         d = zeros(band_limit, 2*band_limit-1, 2*band_limit-1);
0110         d(1,:,:) = ssht_dl(squeeze(d(1,:,:)), band_limit, 0, beta);
0111         for el = 1:band_limit-1
0112             d(el+1,:,:) = ssht_dl(squeeze(d(el,:,:)), band_limit, el, beta);
0113         end
0114     end
0115     % for the case SO3_STORAGE_PADDED:
0116     if (args.Reality == 0) 
0117         if (args.Upsample ~= 0) 
0118           f_cur_lmn_rotated{j-args.J_min+1} = zeros((2*Nj-1)*args.L^2,1); 
0119         else 
0120           f_cur_lmn_rotated{j-args.J_min+1} = zeros((2*Nj-1)*band_limit^2,1);
0121         end 
0122         for el = abs(args.Spin):(band_limit-1) 
0123             for m = -el:el
0124                 if (args.Upsample == 0)   
0125                    ind_lml = so3_elmn2ind(el,m,el,band_limit,Nj);
0126                    ind_l_m_nl = so3_elmn2ind(el,m,-el,band_limit,Nj);
0127                 else  
0128                    ind_lml = so3_elmn2ind(el,m,el,args.L,Nj);
0129                    ind_l_m_nl = so3_elmn2ind(el,m,-el,args.L,Nj);
0130                 end 
0131                 en_max = min(el, Nj-1); 
0132                 for k = -en_max:en_max 
0133                         % Dlmn = exp(-1i*m*alpha) * d(el+1,m+L,n+L) * exp(-1i*n*gamma);
0134                     if (args.Upsample ~= 0)
0135                         Dlkl = exp(-1i*k*alpha) * d(el+1,k+args.L,el+args.L) * exp(-1i*el*gamma);  
0136                         Dlknl = exp(-1i*k*alpha) * d(el+1,k+args.L,-el+args.L) * exp(-1i*(-el)*gamma);
0137                         ind_lmk = so3_elmn2ind(el,m,k,args.L,Nj);
0138                     else
0139                         Dlkl = exp(-1i*k*alpha) * d(el+1,k+band_limit,el+band_limit) * exp(-1i*el*gamma);
0140                         Dlknl = exp(-1i*k*alpha) * d(el+1,k+band_limit,-el+band_limit) * exp(-1i*(-el)*gamma);
0141                         ind_lmk = so3_elmn2ind(el,m,k,band_limit,Nj);
0142                     end
0143                     f_cur_lmn_rotated{j-args.J_min+1}(ind_lmk)= conj(Dlkl) * f_cur_lmn{j-args.J_min+1}(ind_lml)+ ...
0144                                                                     conj(Dlknl)* f_cur_lmn{j-args.J_min+1}(ind_l_m_nl);                                        
0145                 end % end k-loop
0146             end % end m-loop
0147         end % end el-loop
0148     else %i.e. real signals
0149         if (args.Upsample ~= 0) 
0150           f_cur_lmn_rotated{j-args.J_min+1} = zeros(Nj*args.L^2,1);  
0151         else
0152           f_cur_lmn_rotated{j-args.J_min+1} = zeros(Nj*band_limit^2,1);  
0153         end 
0154         for el = 0:(band_limit-1) 
0155             for m = -el:el
0156                 if (args.Upsample == 0)  
0157                    ind_lml = so3_elmn2ind(el,m,el,band_limit,Nj, 'Reality', args.Reality); 
0158                    ind_l_nm_l = so3_elmn2ind(el,-m,el,band_limit,Nj, 'Reality', args.Reality);
0159                 else
0160                    ind_lml = so3_elmn2ind(el,m,el,args.L,Nj, 'Reality', args.Reality) ; 
0161                    ind_l_nm_l = so3_elmn2ind(el,-m,el,args.L,Nj, 'Reality', args.Reality) ;
0162                 end 
0163                 if (mod((m+el),2) == 1) 
0164                     sign = -1; 
0165                 else     
0166                     sign = 1; 
0167                 end 
0168                 en_max = min(el, Nj-1); 
0169                 for k = 0:en_max
0170                      if (args.Upsample ~= 0)
0171                          Dl_k_l = exp(-1i*k*alpha) * d(el+1,k+args.L,el+args.L) * exp(-1i*el*gamma);
0172                          Dl_k_nl = exp(-1i*k*alpha) * d(el+1,k+args.L,-el+args.L) * exp(-1i*-el*gamma);
0173                          ind_lmk = so3_elmn2ind(el,m,k,args.L,Nj,'Reality', args.Reality);
0174                      else
0175                          Dl_k_l = exp(-1i*k*alpha) * d(el+1,k+band_limit,el+band_limit) * exp(-1i*el*gamma);
0176                          Dl_k_nl = exp(-1i*k*alpha) * d(el+1,k+band_limit,-el+band_limit) * exp(-1i*-el*gamma);
0177                          ind_lmk = so3_elmn2ind(el,m,k,band_limit,Nj, 'Reality', args.Reality);
0178                      end 
0179                      f_cur_lmn_rotated{j-args.J_min+1}(ind_lmk)= conj(Dl_k_l) * f_cur_lmn{j-args.J_min+1}(ind_lml)+ ...
0180                                                                  sign*conj(Dl_k_nl)* conj(f_cur_lmn{j-args.J_min+1}(ind_l_nm_l));
0181                 end  % end k-loop
0182             end % end m-loop
0183         end % end el-loop
0184     end % end if (reality)-loop
0185 end %end j-loop
0186 
0187 
0188 % -----------------
0189 % Transform to pixel space:
0190 % -----------------
0191 % Scaling functions in real space:
0192 if (args.Upsample == 0)  
0193      band_limit = min([s2let_bandlimit(args.J_min-1,args.J_min,args.B,args.L) args.L ]);
0194 else
0195      band_limit = args.L ;
0196 end
0197 f_scal = ssht_inverse(f_scal_lm, band_limit,  ...
0198                       'Method', args.Sampling, ...
0199                       'Spin', 0, ...
0200                       'Reality', args.Reality);          
0201                   
0202 % Rotated-curvelets contributions:
0203 % Compute the curvelet coefficients in real space
0204 for j = args.J_min:J,
0205     band_limit = min([s2let_bandlimit(j,args.J_min,args.B,args.L) args.L ]);
0206     Nj = band_limit; 
0207     if (args.Upsample == 0)  
0208         f_cur{j-args.J_min+1} = so3_inverse_direct(f_cur_lmn_rotated{j-args.J_min+1}, band_limit, Nj, ...
0209                                             'Sampling', args.Sampling, 'Reality', args.Reality) ;
0210     else
0211         f_cur{j-args.J_min+1} = so3_inverse_direct(f_cur_lmn_rotated{j-args.J_min+1}, args.L, Nj, ...
0212                                             'Sampling', args.Sampling, 'Reality', args.Reality) ;
0213     end
0214 end
0215 % size(f_cur_lmn_rotated{J-args.J_min+1})
0216 % size(f_cur{J-args.J_min+1})
0217 
0218 
0219 % Clear array
0220 cur_lm = 0; 
0221 scal_l = 0; 
0222 f_cur_lmn =0; 
0223 f_scal_lm =0;
0224 
0225 end