function [outIm] = uncompress(codedIm, lambda, imDims)

% [outIm] = UNCOMPRESS(codedIm, lambda, imDims)
%
% Reconstruct an image from its coded representation (as generated by
% the accompanying function COMPRESS) [Balle,Laparra,Simoncelli,
% ICLR2017].  CODEDIM is a tensor (3D matrix) containing 128 channels
% that form the coded representation. LAMBDA determines the amount of
% compression rate (tradeoff between bits and quality), and can only
% take values 2^k for k in the range [5:12]. It must be the same value
% used to generate CODEDIM.  IMDIMS specifies the dimensions (a [y,x]
% vector) of the reconstructed image, OUTIM, and should be the same as
% the size of the originally coded image.
%
% Relies on convolution code from matlabPyrTools, available at:
%    http://www.cns.nyu.edu/~lcv/software.php 
%
% Published in:
%   End-to-end optimized image compression
%   J Ballé, V Laparra and E P Simoncelli
%   Int'l. Conference on Learning Representations (ICLR2017), April 2017. 
%   https://arxiv.org/abs/1611.01704
%
% See http://www.cns.nyu.edu/~lcv/iclr2017/ for further information.
% 
% Released July, 2017.

% lambda must be 2^n, for n an integer in [5,12].
rl =  2.^min(max(round(log2(lambda)),5),12);  
if (abs(rl-lambda)>0.1)
    warning('Rounding lambda value %.2f to %d', lambda, rl);
end
lambda = rl;

% File directory containing parameter files, assumed relative to
% m-file source code of the current function
myPath = fileparts(mfilename('fullpath'));
paramDir = sprintf('%s/parameters/gray/lambda-%06d/', myPath, lambda);

inRep = codedIm;                        % Initialize input representation to codedIm argument
inDims = [size(inRep,1), size(inRep,2)];
inDepth = size(inRep,3);                % "depth" = number of channels

stage = 0;                              % stage 0: additive constant only
load(sprintf('%ssynthesis-%02d.mat', paramDir, stage), 'c');
const = double(c); clear c;

for n=1:inDepth
    inRep(:,:,n) = inRep(:,:,n) + const(n);
end

for stage = 1:3
    load(sprintf('%ssynthesis-%02d.mat', paramDir, stage), 'h', 'strides', 'c', 'beta', 'gamma');
    filters = double(h); clear h;  fDims = [size(filters,3),size(filters,4)];
    upSampling = double(strides); clear strides;
    const = double(c); clear c;
    beta = double(beta); gamma = double(gamma);
   
    if (inDepth ~= size(filters,1))
        error(sprintf('Stage %d: depth of input %d does not match that of filter set',stage,inDepth));
    end
    outDims = (inDims-1).*upSampling + fDims;  %**** minimal size for valid convolution
    outStart = 1+(fDims-1)/2;
    outStop = outDims-(fDims-1)/2;
    outDepth = size(filters,2);

    % multiplicative re-normalization
    norm2 = reshape(inRep.^2, [prod(inDims), inDepth]);
    norm2 = norm2 * gamma + ones(size(norm2,1),1) * beta;
    inRep = inRep .* reshape(sqrt(norm2), size(inRep));

    outRep = zeros([outDims, outDepth]);
    for n=1:outDepth
        outIm = zeros(outDims);
        for m=1:inDepth
            filt = squeeze(filters(m,n,:,:));  edges = 'reflect1';
            % Upsample and convolve, with params chosen to act like matlab 'valid' convolution
            % NOTE: convolution result is added (in place) into outIm.
            upConv(inRep(:,:,m), filt, edges, upSampling, outStart, outStop, outIm);
        end
        outRep(:,:,n) = outIm + const(n);
    end

    inRep = outRep;                     % set input for next stage to output of this stage
    inDims = outDims;
    inDepth = outDepth;
end

% Trim padding inserted by compress function:
bdry = [floor((size(outRep)-imDims)/2), ceil((size(outRep)-imDims)/2)];
outIm = outRep((bdry(1)+1):(size(outRep,1)-bdry(3)), ((bdry(2)+1):(size(outRep,2)-bdry(4))));

return
