function [img,pyr] = CurvWrapPyramid(Image,scale)

% CurvWrapPyramid - returns an image containing the curvelet coefficients pyramid and the locations of the 4 cardinal angular sectors.
%		    All angular wedge are packed into one angular subband. All atoms are normalized to a unit l2 norm.
%
% Inputs
%     Image         Input image
%     scale	    Coarsest decomposition scale 
%
% Outputs
%     img       Image containing all the curvelet coefficients pyramid. The coefficents are rescaled so that
%       	all atoms in each subband have unit l2 norm.
%     pyr	Cell array of angular sectors pyramid (EAST, WEST, NORTH, SOUTH).
%
  
  C = fdct_wrapping(Image,1,scale);
  
  % Compute L2 norms of curvelets (for noise std computations)
  E = computeL2norm(C);  
  
  m = C{1}{2}(1);
  n = C{1}{2}(2);
  nbscales = length(C);
  
  img = C{1}{1};  img = img/max(max(abs(img))); %normalize
  for sc=2:nbscales
    for w = 1:length(C{sc})
    	C{sc}{w} = C{sc}{w}/E{sc}{w};
    end
    nd = length(C{sc})/4;
    wcnt = 0;
    
    ONE = [];
    for w=1:nd
      ONE = [ONE, C{sc}{wcnt+w}];
    end
    wcnt = wcnt+nd;
    
    TWO = [];
    for w=1:nd
      TWO = [TWO; C{sc}{wcnt+w}];
    end
    wcnt = wcnt+nd;
    
    THREE = [];
    for w=1:nd
      THREE = [C{sc}{wcnt+w}, THREE];
    end
    wcnt = wcnt+nd;
    
    FOUR = [];
    for w=1:nd
      FOUR = [C{sc}{wcnt+w}; FOUR];
    end
    wcnt = wcnt+nd;
    
    [p,q] = size(img);
    [a,b] = size(ONE);
    [g,h] = size(TWO);
    m = 2*a+g;    n = 2*h+b; %size of new image
    scale = max(max( max(max(abs(ONE))),max(max(abs(TWO))) ), max(max(max(abs(THREE))), max(max(abs(FOUR))) )); %scaling factor
    
    pyr{sc-1}=struct('EAST',[],'WEST',[],'NORTH',[],'SOUTH',[]);
    pyr{sc-1}.NORTH = ONE;
    pyr{sc-1}.EAST = TWO;
    pyr{sc-1}.SOUTH = THREE;
    pyr{sc-1}.WEST = FOUR;
    %pind{sc-1}.NORTH{1} = [1,a];  % ONE.
    %pind{sc-1}.NORTH{2} = [h+1,h+b];
    
    %pind{sc-1}.EAST{1} = [a+1,a+g]; % TWO.
    %pind{sc-1}.EAST{2} = [h+b+1,2*h+b];
    
    %pind{sc-1}.SOUTH{1} = [a+g+1,2*a+g]; % THREE.
    %pind{sc-1}.SOUTH{2} = [h+1,h+b];
    
    %pind{sc-1}.WEST{1} = [a+1,a+g]; % FOUR.
    %pind{sc-1}.WEST{2} = [1,h];
    
    new = 0.5 * ones(m,n); %background value
    new(a+1:a+g,1:h) = FOUR /scale;
    new(a+g+1:2*a+g,h+1:h+b) = THREE /scale;
    new(a+1:a+g,h+b+1:2*h+b) = TWO /scale;
    new(1:a,h+1:h+b) = ONE /scale; %normalize
    
    dx = floor((g-p)/2);    dy = floor((b-q)/2);
    
    new(a+1+dx:a+p+dx,h+1+dy:h+q+dy) = img;
    
    img = new;
  end

pyr{end+1}=C{1}; % LP and size of original image.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function E=computeL2norm(coeffs)
% Compute norm of curvelets (exact)

F = ones(coeffs{1}{2});
X = fftshift(ifft2(F)) * sqrt(prod(size(F))); 	% Appropriately normalized Dirac
C = fdct_wrapping(X,1,length(coeffs)); 		% Get the curvelets
nb=prod(size(C{1}{1}));

E = cell(size(C));
for j=1:length(C)
  E{j} = cell(size(C{j}));
  for w=1:length(C{j})
    A = C{j}{w};
    E{j}{w} = sqrt(sum(sum(A.*conj(A))) / prod(size(A)));
    nb=nb+prod(size(C{j}{w}));
  end
end
%nb/prod(size(X))

%
% Part of DBlockToolbox Version:100
% Written by: Jalal Fadili, GREYC CNRS ENSICAEN
%             Christophe Chesneau, LMNO University of Caen
%             Jean-Luc Starck, CEA-CNRS
% Created June 2008
% This material is distributed under CeCiLL licence.
% E-mail: Jalal.Fadili@greyc.ensicaen.fr
%