mefcv2.m

This package is a free collection of Matlab routines for computing wave atom transforms in one, two

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);