mecv.rar_curvelet_curvelet matlab_curvelet transform_curvelet 变换

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