function JPG = jpeg_encode(im_uint8, subsample_ratio, quant_depth)
%JPEG_ENCODE  encode image into jpeg-like format
%   JPG = JPEG_ENCODE(im_uint8, subsample_ratio, quant_depth)
%  Parameters:
%    im_uint8           H x W x 3  size uint8 array
%    subsample_ratio    ratio of chroma subsampling (integer)
%    quant_depth        quantization depth (0 to inf)

%% Global coding constants

% uint pixel data is represented in the interval 0 : 255
full_range = 256;

% Macro block size
block_size = 8;

% DCT is a 8x8 matrix. The direct cosine transform (DCT) of a column vector
% is computed as   DCT * a
DCT = dct_matrix(block_size);

% Q0 : Quantization matrix for default (level 1) quantizing
Q0 = [
    16 11 10 16  24  40  51  61
    12 12 14 19  26  58  60  55
    14 13 16 24  40  57  69  56
    14 17 22 29  51  87  80  62
    18 22 37 56  68 109 103  77
    24 35 55 64  81 104 113  92
    49 64 78 87 103 121 120 101
    72 92 95 98 112 100 103  99
    ];

%% Preprocess image for coding
% Convert to floating point for accurate signal processing
im = double(im_uint8);

%% Colorspace conversion from RGB to YCbCr
R = im(:,:,1);
G = im(:,:,2);
B = im(:,:,3);
[Y, Cb, Cr] = rgb_to_ycbcr(R, G, B);

%% Downsample chroma components
method = 'rect';
Cb_sub = downsample(Cb, subsample_ratio, method);
Cr_sub = downsample(Cr, subsample_ratio, method);

%% Divide into macroblocks and form contiguous sequence of blocks
Y_mb = matrix_to_blocks(Y, block_size);
Cb_mb = matrix_to_blocks(Cb_sub, block_size);
Cr_mb = matrix_to_blocks(Cr_sub, block_size);

Y_coded = encode_blocks(Y_mb);
Cb_coded = encode_blocks(Cb_mb);
Cr_coded = encode_blocks(Cr_mb);

%% Collect coded data into a structure
JPG.image_size = size(im_uint8);        % original image dimensions
JPG.subsample_ratio = subsample_ratio;  % chroma subsampling ratio
JPG.Y_coded = Y_coded;                  % the coded block data
JPG.Cb_coded = Cb_coded;                % the coded block data
JPG.Cr_coded = Cr_coded;                % the coded block data
JPG.Q = Q0 * quant_depth;               % the quantizing matrix

    function coded = encode_blocks(blocks)
        
        %% Perform DCT for each block and quantize
        % the actual quantizing matrix is a mutiple of the default
        Q = Q0 * quant_depth;
        
        % preallocate space for the coded blocks
        coded = cell(size(blocks));
        
        for k = 1 : numel(blocks)
            b = blocks{k};      % get block data
            b = b - full_range / 2;     % convert to zero offset
            b = DCT * b * DCT'; % perform DCT in both directions
            b = round(b ./ Q);  % quantize
            coded{k} = b;
        end
    end
end
