曲波变换重构

function x = ifdct_wrapping（C is_real M N）

% ifdct_wrapping.m - Inverse Fast Discrete Curvelet Transform via wedge wrapping - Version 1.0
% This is in fact the adjoint also the pseudo-inverse
%
% Inputs
%   C           Cell array containing curvelet coefficients （see
%               description in fdct_wrapping.m）
%   is_real     As used in fdct_wrapping.m
%   M N        Size of the image to be recovered （not necessary if finest
%               = 2）
%
% Outputs
%   x           M-by-N matrix
%
% See also fdct_wrapping.m
%
% By Laurent Demanet 2004

% Initialization
nbscales = length（C）;
nbangles_coarse = length（C{2}）;
nbangles = [1 nbangles_coarse .* 2.^（ceil（（nbscales-（nbscales:-1:2））/2））];
if length（C{end}） == 1 finest = 2; else finest = 1; end;
if finest == 2 nbangles（nbscales） = 1; end;
if nargin < 2 is_real = 0; end;
if nargin < 4
    if finest == 1 error（‘Syntax: IFCT_wrapping（CMN） where the matrix to be recovered is M-by-N‘）; end;
    [N1N2] = size（C{end}{1}）;
else
    N1 = M;
    N2 = N;
end;

M1 = N1/3;
M2 = N2/3;

if finest == 1;
    
    bigN1 = 2*floor（2*M1）+1;
    bigN2 = 2*floor（2*M2）+1;
    X = zeros（bigN1bigN2）;

    % Initialization: preparing the lowpass filter at finest scale
    window_length_1 = floor（2*M1） - floor（M1） - 1 - （mod（N13）==0）;
    window_length_2 = floor（2*M2） - floor（M2） - 1 - （mod（N23）==0）;
    coord_1 = 0:（1/window_length_1）:1;
    coord_2 = 0:（1/window_length_2）:1;
    [wl_1wr_1] = fdct_wrapping_window（coord_1）;
    [wl_2wr_2] = fdct_wrapping_window（coord_2）;
    lowpass_1 = [wl_1 ones（12*floor（M1）+1） wr_1];
    if mod（N13）==0 lowpass_1 = [0 lowpass_1 0]; end;
    lowpass_2 = [wl_2 ones（12*floor（M2）+1） wr_2];
    if mod（N23）==0 lowpass_2 = [0 lowpass_2 0]; end;
    lowpass = lowpass_1‘*lowpass_2;

    scales = nbscales:-1:2;
   
else

    M1 = M1/2;
    M2 = M2/2;
    
    bigN1 = 2*floor（2*M1）+1;
    bigN2 = 2*floor（2*M2）+1;
    X = zeros（bigN1bigN2）;
    
    window_length_1 = floor（2*M1） - floor（M1） - 1;
    window_length_2 = floor（2*M2） - floor（M2） - 1;
    coord_1 = 0:（1/window_length_1）:1;
    coord_2 = 0:（1/window_length_2）:1;
    [wl_1wr_1] = fdct_wrapping_window（coord_1）;
    [wl_2wr_2] = fdct_wrapping_window（coord_2）;
    lowpass_1 = [wl_1 ones（12*floor（M1）+1） wr_1];
    lowpass_2 = [wl_2 ones（12*floor（M2）+1） wr_2];
    lowpass = lowpass_1‘*lowpass_2;
    hipass_finest = sqrt（1 - lowpass.^2）;
    
    scales = （nbscales-1）:-1:2;
    
end;

% Loop: pyramidal reconstruction

Xj_topleft_1 = 1;
Xj_topleft_2 = 1;
for j = scales

    M1 = M1/2;
    M2 = M2/2;
    window_length_1 = floor（2*M1） - floor（M1） - 1;
    window_length_2 = floor（2*M2） - floor（M2） - 1;
    coord_1 = 0:（1/window_length_1）:1;
    coord_2 = 0:（1/window_length_2）:1;
    [wl_1wr_1] = fdct_wrapping_window（coord_1）;
    [wl_2wr_2] = fdct_wrapping_window（coord_2）;
    lowpass_1 = [wl_1 ones（12*floor（M1）+1） wr_1];
    lowpass_2 = [wl_2 ones（12*floor（M2）+1） wr_2];
    lowpass_next = lowpass