function c = mefcv2(x,N1,N2,ns,nag)% mefcv2 - 2D forward mirror-extended curvelet transform% -----------------% INPUT% --% x is an N1-by-N2 matrix. % --% N1, N2 are positive integers.% --% ns is the number of levels, including the coarsest level. ns =% ceil(log2(min(N1,N2)) - 3) is commonly used.% --% nag. 2*pi/nag is the size of the spanning angle of each wedge.% nag is required to be a multiple of 8 and nag = 16 is often used.% -----------------% OUTPUT% --% c is a cell array which contains the curvelets coefficients. If% tp=='ortho', then c{j}{l}(n1,n2) is the coefficient at scale j,% direction l and spatial index (n1,n2). The directional index l% iterates through the wedges in the first quadrant. Notice that, for% the mirror-extended wave atoms, the spatial indices wrap around once.% -----------------% Written by Lexing Ying and Laurent Demanet, 2007    if(mod(nag,4)~=0)    error('wrong');  end    %1. dct  fd = dct2(x);    %2. scatter  E1 = ceil(N1/3);  E2 = ceil(N2/3);  %E1 = 0;  E2 = 0;  fd = mescatter(fd,E1);  fd = mescatter(fd',E2)';  A1 = size(fd,1);    A2 = size(fd,2);    G1 = 4/3*N1;  G2 = 4/3*N2;    c = cell(ns,1);  ttl = 0;    for s=ns:-1:2    %get ring    R1 = 2^(s-ns)*G1;    R2 = 2^(s-ns)*G2;        idx1 = [ceil(-R1):floor(R1)];    [wl,wr] = cvwindow((idx1+R1/1)/(R1/2));    tmpa = wl;    [wl,wr] = cvwindow((idx1-R1/2)/(R1/2));    tmpb = wr;    coef1 = tmpa.*tmpb;    idx2 = [ceil(-R2):floor(R2)];    [wl,wr] = cvwindow((idx2+R2/1)/(R2/2));    tmpa = wl;    [wl,wr] = cvwindow((idx2-R2/2)/(R2/2));    tmpb = wr;    coef2 = tmpa.*tmpb;    lowpass = coef1'*coef2;        idx1 = [ceil(-R1):floor(R1)];    [wl,wr] = cvwindow((idx1+R1/2)/(R1/4));    tmpa = wl;    [wl,wr] = cvwindow((idx1-R1/4)/(R1/4));    tmpb = wr;    coef1 = tmpa.*tmpb;    idx2 = [ceil(-R2):floor(R2)];    [wl,wr] = cvwindow((idx2+R2/2)/(R2/4));    tmpa = wl;    [wl,wr] = cvwindow((idx2-R2/4)/(R2/4));    tmpb = wr;    coef2 = tmpa.*tmpb;    tmppass = coef1'*coef2;    hghpass = sqrt(1-tmppass.^2);        pass = lowpass.*hghpass;    [M1,M2] = size(pass);    fh = zeros(M1,M2);    fh(mod(idx1,M1)+1,mod(idx2,M2)+1) = pass .* fd(mod(idx1,A1)+1,mod(idx2,A2)+1);    %ggg = ggg + norm(fh(:))^2;        nbangles = nag*2^(ceil((s-2)/2));        %c{s} = fcv2_sepangle(R1,R2,fh,nbangles);    %---------    [M1,M2] = size(fh);    nd = nbangles/4;    cs = cell(nd,1);    W1 = 2*R1/nd;  W2 = 2*R2/nd;        %take only first quadrant    cnt = 1;    for g=nd/2:nd-1      xs = R1/4-(W1/2)/4;    xe = R1;      ys = -R2 + (2*g-1)*W2/2;		ye = -R2 + (2*g+3)*W2/2;      xn = ceil(xe-xs);    yn = ceil(ye-ys);      if(g==0)        thts = atan2(-1.0, 1.0-1.0/nd);        thtm = atan2(-1.0+1.0/nd, 1.0);        thte = atan2(-1.0+3.0/nd, 1.0);      elseif(g==nd-1)        thts = atan2(-1.0+(2.0*g-1.0)/nd, 1.0);        thtm = atan2(-1.0+(2.0*g+1.0)/nd, 1.0);        thte = atan2(1.0, 1.0-1.0/nd);      else        thts = atan2(-1.0+(2.0*g-1.0)/nd, 1.0);        thtm = atan2(-1.0+(2.0*g+1.0)/nd, 1.0);        thte = atan2(-1.0+(2.0*g+3.0)/nd, 1.0);      end      %fprintf(1,'%d %d %d\n',thts,thtm,thte);      R21 = R2/R1;      wpdata = zeros(xn,yn);      for xcur=ceil(xs):xe        yfm = ceil( max([-R2, R21*xcur*tan(thts)]) );        yto = floor( min([R2, R21*xcur*tan(thte)]) );        ycur = yfm:yto;        thtcur = atan2(ycur/R2,xcur/R1);        [al,ar] = cvwindow((thtcur-thts)/(thtm-thts));        [bl,br] = cvwindow((thtcur-thtm)/(thte-thtm));        pou = al.*br;        wpdata(mod(xcur,xn)+1,mod(ycur,yn)+1) = fh(mod(xcur,M1)+1,mod(ycur,M2)+1) .* pou;      end      cs{cnt} = ifft2(wpdata) * sqrt(numel(wpdata));      ttl = ttl + norm(cs{cnt}(:))^2;      cnt = cnt+1;    end        for f=nd-1:-1:nd/2      ys = R2/4-(W2/2)/4;		  ye = R2;      xs = -R1 + (2*f-1)*W1/2;		  xe = -R1 + (2*f+3)*W1/2;      xn = ceil(xe-xs);		  yn = ceil(ye-ys);      if(f==0)        phis = atan2(-1.0, 1.0-1.0/nd);        phim = atan2(-1.0+1.0/nd, 1.0);        phie = atan2(-1.0+3.0/nd, 1.0);      elseif(f==nd-1)        phis = atan2(-1.0+(2.0*f-1.0)/nd, 1.0);        phim = atan2(-1.0+(2.0*f+1.0)/nd, 1.0);        phie = atan2(1.0, 1.0-1.0/nd);      else        phis = atan2(-1.0+(2.0*f-1.0)/nd, 1.0);        phim = atan2(-1.0+(2.0*f+1.0)/nd, 1.0);        phie = atan2(-1.0+(2.0*f+3.0)/nd, 1.0);      end      %fprintf(1,'%d %d %d\n',phis,phim,phie);      R12 = R1/R2;      wpdata = zeros(xn,yn);      for ycur=ceil(ys):ye        xfm = ceil( max([-R1, R12*ycur*tan(phis)]) );        xto = floor( min([R1, R12*ycur*tan(phie)]) );        xcur = xfm:xto;        phicur = atan2(xcur/R1, ycur/R2);        [al,ar] = cvwindow((phicur-phis)/(phim-phis));        [bl,br] = cvwindow((phicur-phim)/(phie-phim));        pou = al.*br;        wpdata(mod(xcur,xn)+1,mod(ycur,yn)+1) = fh(mod(xcur,M1)+1,mod(ycur,M2)+1) .* pou';      end      cs{cnt} = ifft2(wpdata) * sqrt(numel(wpdata));      ttl = ttl + norm(cs{cnt}(:))^2;      cnt = cnt+1;    end    c{s} = cs;      %fprintf(1,'%d %d\n', ttl, ggg/4);  end      %do the first level  if(1)    s = 1;    R1 = 2^(s-ns)*G1;    R2 = 2^(s-ns)*G2;    idx1 = [ceil(-R1):floor(R1)];    [wl,wr] = cvwindow((idx1+R1/1)/(R1/2));    tmpa = wl;    [wl,wr] = cvwindow((idx1-R1/2)/(R1/2));    tmpb = wr;    coef1 = tmpa.*tmpb;    idx2 = [ceil(-R2):floor(R2)];    [wl,wr] = cvwindow((idx2+R2/1)/(R2/2));    tmpa = wl;    [wl,wr] = cvwindow((idx2-R2/2)/(R2/2));    tmpb = wr;    coef2 = tmpa.*tmpb;    pass = coef1'*coef2;    [M1,M2] = size(pass);        fh = zeros(M1,M2);    fh(mod(idx1,M1)+1,mod(idx2,M2)+1) = pass .* fd(mod(idx1,A1)+1,mod(idx2,A2)+1);    %ggg = ggg + norm(fh(:))^2;        cs = cell(1,1);    K1 = M1+1;    K2 = M2+1;    tmp = zeros(K1,K2);    tmp(mod(idx1,K1)+1,mod(idx2,K2)+1) = fh(mod(idx1,M1)+1,mod(idx2,M2)+1);    tmp = mecombine(tmp,0);    tmp = mecombine(tmp',0)';    tmp = tmp/4;    cs{1} = ifft2(tmp) * sqrt(numel(tmp));    ttl = ttl + norm(cs{1}(:))^2;    c{s} = cs;        %fprintf(1,'%d\n', ttl);  end              %   f0 = dct2(x);%   o = zeros(N,1);%   f1 = [f0 o -f0(:,end:-1:2)];%   o = zeros(1,2*N);%   f2 = [f1; o; -f1(end:-1:2,:)];%   p = ones(2*N,1);%   p(1) = sqrt(2);%   f = (p * p') .* f2; %frequency content  %   tx = ifft2(f) * sqrt(prod(size(f)));  %   %call curvelet transform%   addpath /Users/lexing/projects/curve_lab/CurveLab-2.0.2/fdct_wrapping_cpp/mex;%   m = 2*N;  n = 2*N;  nbscales = floor(log2(min(m,n)))-3;  nbangles_coarse = 16;  allcurvelets = 1;%   tc = fdct_wrapping_mex(m,n,nbscales,nbangles_coarse,allcurvelets,double(tx));%   rmpath /Users/lexing/projects/curve_lab/CurveLab-2.0.2/fdct_wrapping_cpp/mex;  %   c = cell(1,length(tc));%   for s=2:length(tc)%     nw = length(tc{s});%     c{s} = tc{s}(1+3/8*nw:5/8*nw);%   end%   tmp0 = fft2(tc{1}{1}) / sqrt(prod(size(tc{1}{1})));%   [m,n] = size(tmp0);%   mp = (m+1)/2;  np = (n+1)/2;%   tmp = tmp0(1:mp,1:np);%   q = ones(mp,1);%   q(1) = 1/sqrt(2);%   tmp = (q*q') .* tmp;%   c{1}{1} = idct2(tmp);