function x = meicv2(c,N1,N2,ns,nag)% meicv2 - 2D inverse mirror-extended curvelet transform% -----------------% INPUT% --% 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.% --% 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 is the number of angles used for the second coarsest level.% nag is required to be a multiple of 4 and nag = 16 is often used.% -----------------% OUTPUT% --% x is an N1-by-N2 matrix. % -----------------% Written by Lexing Ying and Laurent Demanet, 2007    E1 = ceil(N1/3);  E2 = ceil(N2/3);  %E1 = 0;  E2 = 0;  A1 = 2*(N1+E1);  A2 = 2*(N2+E2);  fd = zeros(A1,A2);    G1 = 4/3*N1;  G2 = 4/3*N2;    for s=ns:-1:2    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);    %get fh        nbangles = nag*2^(ceil((s-2)/2));    %---------    [M1,M2] = size(fh);    nd = nbangles/4;    cs = c{s};    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 = fft2(cs{cnt}) / sqrt(numel(cs{cnt}));      cnt = cnt+1;      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;        fh(mod(xcur,M1)+1,mod(ycur,M2)+1) = fh(mod(xcur,M1)+1,mod(ycur,M2)+1) + wpdata(mod(xcur,xn)+1,mod(ycur,yn)+1) .* pou;      end    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 = fft2(cs{cnt}) / sqrt(numel(cs{cnt}));      cnt = cnt+1;      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;        fh(mod(xcur,M1)+1,mod(ycur,M2)+1) = fh(mod(xcur,M1)+1,mod(ycur,M2)+1) + wpdata(mod(xcur,xn)+1,mod(ycur,yn)+1) .* pou';      end    end        %put back into fd    fd(mod(idx1,A1)+1,mod(idx2,A2)+1) = fd(mod(idx1,A1)+1,mod(idx2,A2)+1) + pass .* fh(mod(idx1,M1)+1,mod(idx2,M2)+1);  end    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);        cs = c{s};    tmp = fft2(cs{1}) / sqrt(numel(cs{1}));    tmp = tmp/4;    tmp = mescatter(tmp,0);    tmp = mescatter(tmp',0)';    [K1,K2] = size(tmp);    fh = zeros(M1,M2);    fh(mod(idx1,M1)+1,mod(idx2,M2)+1) = tmp(mod(idx1,K1)+1,mod(idx2,K2)+1);        fd(mod(idx1,A1)+1,mod(idx2,A2)+1) = fd(mod(idx1,A1)+1,mod(idx2,A2)+1) + pass .* fh(mod(idx1,M1)+1,mod(idx2,M2)+1);  end    fd = mecombine(fd,E1);  fd = mecombine(fd',E2)';  x = idct2(fd);        %   tc = cell(1,length(c));%   for s=2:length(c)%     nw = length(c{s})*4;%     [ni,nj] = size(c{s}{1});%     tc{s} = cell(1,nw);%     for w=[1:nw/4 nw/2+1:3*nw/4]%       tc{s}{w} = zeros(nj,ni);%     end%     for w=[1+nw/4:nw/2 1+3*nw/4:nw]%       tc{s}{w} = zeros(ni,nj);%     end%     tc{s}(1+3/8*nw:5/8*nw) = c{s};%   end%   tmp = dct2(c{1}{1});%   [mp,np] = size(tmp);%   q = ones(mp,1);%   q(1) = 1/sqrt(2);%   tmp = (q*q') .* tmp;%   tmp0 = zeros(mp*2-1,np*2-1);%   tmp0(1:mp,1:np) = tmp;%   tc{1}{1} = ifft2(tmp0) * sqrt(prod(size(tmp0)));  %   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;%   tx = ifdct_wrapping_mex(m,n,nbscales,nbangles_coarse,allcurvelets,tc);%   rmpath /Users/lexing/projects/curve_lab/CurveLab-2.0.2/fdct_wrapping_cpp/mex;  %   f = fft2(tx) / sqrt(prod(size(tx)));  %   p = ones(2*N,1); p(1) = sqrt(2);%   f2 = (p * p') .* f;%   f1 = f2(1:N,:);%   f1(end:-1:2,:) = f1(end:-1:2,:) - f2(end/2+2:end,:);%   f0 = f1(:,1:N);%   f0(:,end:-1:2) = f0(:,end:-1:2) - f1(:,end/2+2:end);%  %   x = idct2(f0);