s2let_transform_curvelet_synthesis_cur2lm

 s2let_transform_curvelet_synthesis_cur2lm
 Compute curvelet transform:  
 input in curvelet space (i.e. hamonic to Wigner space via SO3_forward_direct)
 output in harmonic space (i.e. Wigner to harmonic space via synthesis_lmn2lm) .

 Default usage :

 flm_rec = s2let_transform_curvelet_synthesis_lm2cur(f_cur, f_scal,  <options>)

 f_cur contains the input curvelet contributions -- MW sampling,
 f_scal contains the input scaling contributions -- MW sampling,
 flm_rec is the output field = flm_cur_syn+ flm_scal_syn

 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 flm_rec = s2let_transform_curvelet_synthesis_cur2lm(f_cur, f_scal,  varargin)
0002 
0003 % s2let_transform_curvelet_synthesis_cur2lm
0004 % Compute curvelet transform:
0005 % input in curvelet space (i.e. hamonic to Wigner space via SO3_forward_direct)
0006 % output in harmonic space (i.e. Wigner to harmonic space via synthesis_lmn2lm) .
0007 %
0008 % Default usage :
0009 %
0010 % flm_rec = s2let_transform_curvelet_synthesis_lm2cur(f_cur, f_scal,  <options>)
0011 %
0012 % f_cur contains the input curvelet contributions -- MW sampling,
0013 % f_scal contains the input scaling contributions -- MW sampling,
0014 % flm_rec is the output field = flm_cur_syn+ flm_scal_syn
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 % S2LET package to perform Wavelet Transform on the Sphere.
0036 % Copyright (C) 2012-2016  Boris Leistedt, Jennifer Chan & Jason McEwen
0037 % See LICENSE.txt for license details
0038 % -----------------------------------------------------------
0039 
0040 len = size(f_cur);
0041 temp = f_cur{len};
0042 sz = size(temp);
0043 Lguessed = sz(2);
0044 
0045 p = inputParser;
0046 p.addRequired('f_cur');
0047 p.addRequired('f_scal', @isnumeric);
0048 p.addParamValue('B', 2, @isnumeric);
0049 p.addParamValue('L', Lguessed, @isnumeric);                           
0050 p.addParamValue('J_min', 0, @isnumeric);
0051 p.addParamValue('Spin', 0, @isnumeric);
0052 p.addParamValue('Upsample', false, @islogical);
0053 p.addParamValue('Reality', false, @islogical);
0054 p.addParamValue('SpinLowered', false, @islogical);
0055 p.addParamValue('SpinLoweredFrom', 0, @isnumeric);
0056 p.addParamValue('Sampling', 'MW', @ischar);
0057 p.parse(f_cur, f_scal, varargin{:});
0058 args = p.Results;
0059 
0060 J = s2let_jmax(args.L, args.B);
0061 
0062 % ---------------
0063 % Tile curvelets:
0064 % ---------------
0065 % Call curvelet- and scaling-function- generating functions
0066 [cur_lm scal_l] = s2let_curvelet_tiling(args.B, args.L, args.J_min,...
0067                                         'Spin', args.Spin,...
0068                                         'SpinLowered', args.SpinLowered, ...
0069                                         'SpinLoweredFrom', args.SpinLoweredFrom);
0070 
0071 % -----------------
0072 % Signal synthesis: (Transform to lmn space, then reconstruct the signal in harmonic space)
0073 % -----------------
0074 % Scaling function contribution:
0075 if (args.Upsample == 0)  
0076      band_limit = min([ s2let_bandlimit(args.J_min-1,args.J_min,args.B,args.L) args.L ]);
0077 else
0078      band_limit = args.L ;
0079 end
0080 f_scal_lm_syn = ssht_forward(f_scal, band_limit, ...
0081                             'Method', args.Sampling,...
0082                             'Spin', 0, ...
0083                             'Reality', args.Reality);
0084   
0085 % Curvelet contribution:
0086 for j = args.J_min:J,
0087     band_limit = min([ s2let_bandlimit(j,args.J_min,args.B,args.L) args.L ]);
0088     % Nj = orientational band-limit at j-th scale
0089     Nj = band_limit;
0090     if (args.Upsample == 0)  
0091         f_cur_lmn_syn{j-args.J_min+1} = so3_forward_direct(f_cur{j-args.J_min+1} , band_limit, Nj, ...
0092                                                     'Reality', args.Reality, 'Sampling', args.Sampling);
0093     else  
0094         f_cur_lmn_syn{j-args.J_min+1} = so3_forward_direct(f_cur{j-args.J_min+1} , args.L, Nj, ...
0095                                                     'Reality', args.Reality, 'Sampling', args.Sampling);
0096     end
0097 end
0098 
0099 % Rotate the Wigner coefficients f_cur_lmn (such the curvelets centered at the North pole)
0100 % Exploit the property of curvelets that cur_ln = (cur_ll)*(delta_ln)
0101 % such that cur_lmk_rotated = cur_lml*conj(Dlkl(evaluated at the desired rotation angle))
0102 % ---------------
0103 % Define Euler angles
0104 % (for rotating the curvelets to the north pole):
0105 % ---------------
0106 alpha = pi ;
0107 gamma = 0 ;
0108 
0109 for j = args.J_min:J,
0110     band_limit = min([ s2let_bandlimit(j,args.J_min,args.B,args.L) args.L ]);
0111     Nj = band_limit;
0112 % ---------------
0113 % Precompute Wigner small-d functions
0114 % denoted here as d (in the paper: d_lmn for all el, m, n evaluated at beta).
0115 % They are indexed d(el,m,n).
0116 % Alpha and gamma are the other two rotation angles.
0117 % ---------------
0118     if (args.Upsample ~= 0)
0119         beta = pi-acos(-args.Spin/args.B^j);
0120         d = zeros(args.L, 2*args.L-1, 2*args.L-1);
0121         d(1,:,:) = ssht_dl(squeeze(d(1,:,:)), args.L, 0, beta);
0122         for el = 1:args.L-1
0123             d(el+1,:,:) = ssht_dl(squeeze(d(el,:,:)), args.L, el, beta);
0124         end
0125     else
0126         beta = pi-acos(-args.Spin/args.B^j);
0127         d = zeros(band_limit, 2*band_limit-1, 2*band_limit-1);
0128         d(1,:,:) = ssht_dl(squeeze(d(1,:,:)), band_limit, 0, beta);
0129         for el = 1:band_limit-1
0130            d(el+1,:,:) = ssht_dl(squeeze(d(el,:,:)), band_limit, el, beta);
0131         end
0132     end
0133   % for the case SO3_STORAGE_PADDED:
0134     if (args.Reality == 0)  
0135         if (args.Upsample ~= 0) 
0136             f_cur_lmn_syn_rotated{j-args.J_min+1} = zeros((2*Nj-1)*args.L^2,1); 
0137         else 
0138             f_cur_lmn_syn_rotated{j-args.J_min+1} = zeros((2*Nj-1)*band_limit^2,1);
0139         end 
0140         for el = abs(args.Spin):(band_limit-1) 
0141             for m = -el:el
0142                 if (args.Upsample == 0)  
0143                          ind_lml = so3_elmn2ind(el,m,el,band_limit,Nj);
0144                          ind_lmnl = so3_elmn2ind(el,m,-el,band_limit,Nj);
0145                 else
0146                          ind_lml = so3_elmn2ind(el,m,el,args.L,Nj);
0147                          ind_lmnl = so3_elmn2ind(el,m,-el,args.L,Nj);
0148                 end 
0149                 en_max = min(el, Nj-1); 
0150                 for en = -en_max:en_max
0151                      %  Dlmn = exp(-1i*m*alpha) * d(el+1,m+L,n+L) * exp(-1i*n*gamma);
0152                     if (args.Upsample ~= 0)
0153                         Dlln = exp(-1i*el*alpha) * d(el+1,el+args.L,en+args.L) * exp(-1i*en*gamma);
0154                         Dl_nl_n = exp(-1i*-el*alpha) * d(el+1,-el+args.L,en+args.L) * exp(-1i*en*gamma);
0155                         ind_lmn = so3_elmn2ind(el,m,en,args.L,Nj);
0156                     else
0157                         Dlln = exp(-1i*el*alpha) * d(el+1,el+band_limit,en+band_limit) * exp(-1i*en*gamma);
0158                         Dl_nl_n = exp(-1i*-el*alpha) * d(el+1,-el+band_limit,en+band_limit) * exp(-1i*en*gamma);
0159                         ind_lmn = so3_elmn2ind(el,m,en,band_limit,Nj);
0160                      end  
0161                      f_cur_lmn_syn_rotated{j-args.J_min+1}(ind_lml)=  f_cur_lmn_syn_rotated{j-args.J_min+1}(ind_lml)+ ...
0162                                                                       conj(Dlln) * f_cur_lmn_syn{j-args.J_min+1}(ind_lmn); 
0163                      f_cur_lmn_syn_rotated{j-args.J_min+1}(ind_lmnl)= f_cur_lmn_syn_rotated{j-args.J_min+1}(ind_lmnl) + ...
0164                                                                       conj(Dl_nl_n) * f_cur_lmn_syn{j-args.J_min+1}(ind_lmn); 
0165                 end % end n-loop
0166             end  % end m-loop
0167         end % end el-loop
0168     else % real setting:
0169         if (args.Upsample ~= 0) 
0170             f_cur_lmn_syn_rotated{j-args.J_min+1} = zeros(Nj*args.L^2,1);  
0171         else
0172             f_cur_lmn_syn_rotated{j-args.J_min+1} = zeros(Nj*band_limit^2,1);  
0173         end
0174         for el = 0:(band_limit-1) 
0175             for m = -el:el
0176                 % (n=0) terms
0177                 if (args.Upsample == 0)  
0178                     ind_lml = so3_elmn2ind(el,m,el,band_limit,Nj, 'Reality', args.Reality);
0179                     ind_lmnzero = so3_elmn2ind(el,m,0,band_limit,Nj, 'Reality', args.Reality);
0180                     Dl_l_nzero = exp(-1i*el*alpha) * d(el+1,el+band_limit,0+band_limit) * exp(-1i*0*gamma);
0181                 else
0182                     ind_lml = so3_elmn2ind(el,m,el,args.L,Nj, 'Reality', args.Reality);
0183                     ind_lmnzero = so3_elmn2ind(el,m,0,args.L,Nj, 'Reality', args.Reality);
0184                     Dl_l_nzero = exp(-1i*el*alpha) * d(el+1,el+args.L,0+args.L) * exp(-1i*0*gamma);
0185                 end
0186                 f_cur_lmn_syn_rotated{j-args.J_min+1}(ind_lml)=  conj(Dl_l_nzero)*f_cur_lmn_syn{j-args.J_min+1}(ind_lmnzero);                           
0187                 % (n> 0) terms
0188                 en_max = min(el, Nj-1); 
0189                 for en = 1:en_max
0190                     if (args.Upsample ~= 0)
0191                         Dl_l_n = exp(-1i*el*alpha) * d(el+1,el+args.L,en+args.L) * exp(-1i*en*gamma);
0192                         Dl_l_nn = exp(-1i*el*alpha) * d(el+1,el+args.L,-en+args.L) * exp(-1i*-en*gamma);
0193                         ind_lmn = so3_elmn2ind(el,m,en,args.L,Nj, 'Reality', args.Reality);
0194                         ind_l_nm_n = so3_elmn2ind(el,-m,en,args.L,Nj, 'Reality', args.Reality);
0195                     else
0196                         Dl_l_n = exp(-1i*el*alpha) * d(el+1,el+band_limit,en+band_limit) * exp(-1i*en*gamma);
0197                         Dl_l_nn = exp(-1i*el*alpha) * d(el+1,el+band_limit,-en+band_limit) * exp(-1i*-en*gamma);
0198                         ind_lmn = so3_elmn2ind(el,m,en,band_limit,Nj, 'Reality', args.Reality);
0199                         ind_l_nm_n = so3_elmn2ind(el,-m,en,band_limit,Nj, 'Reality', args.Reality);
0200                     end
0201                     if (mod((m+en),2) == 1)
0202                            sign = -1; 
0203                     else 
0204                            sign = 1; 
0205                     end 
0206                     f_cur_lmn_syn_rotated{j-args.J_min+1}(ind_lml)=  f_cur_lmn_syn_rotated{j-args.J_min+1}(ind_lml)+ ...
0207                                                                      conj(Dl_l_n) * f_cur_lmn_syn{j-args.J_min+1}(ind_lmn)+ ... 
0208                                                                      sign*conj(Dl_l_nn)*conj(f_cur_lmn_syn{j-args.J_min+1}(ind_l_nm_n)); 
0209                 end % end en-loop
0210             end  % end m-loop
0211         end % end el-loop
0212     end % end if Reality loop
0213 end %end j-loop
0214 
0215 % Reconstruct the signals in harmonic space and perform Wigner transform
0216 flm_rec = s2let_transform_curvelet_synthesis_lmn2lm(f_cur_lmn_syn_rotated, f_scal_lm_syn, cur_lm, scal_l,...
0217                                                    'B', args.B, 'L', args.L, ...
0218                                                    'Spin', args.Spin, ...
0219                                                    'J_min', args.J_min, ...
0220                                                    'Upsample', args.Upsample,...
0221                                                    'Reality', args.Reality,...
0222                                                    'SpinLowered', args.SpinLowered, ...
0223                                                    'SpinLoweredFrom', args.SpinLoweredFrom,...
0224                                                    'Sampling', args.Sampling );
0225 % Clear array
0226 cur_lm = 0; 
0227 scal_l = 0; 
0228 f_cur_lmn_syn =0; 
0229 f_scal_lm_syn =0;
0230 f_cur_lmn_syn_rotated =0;                     
0231 end