Codificacion RLE/VLC y análisis de compresión espacial

jpeg.m

%
% Tema: Codificación RLE/VLC y análisis de compresión espacial
% Catedra: Sistemas Avanzados de Televisión y Técnicas Audiovisuales
% Profesor: Miguel Seoane
% Alumno: Tomás Agustín González
% Referencias:
%   - Apuntes de cátedra
%   - http://www.digicamsoft.com/itu/itu-t81-1.html
%   - http://www.stanford.edu/class/ee398a/handouts/lectures/08-JPEG.pdf
%   - http://www.mathworks.com/matlabcentral/newsreader
%   - http://www.mathworks.com/help/index.html
%   - http://www.globalspec.com/reference/39556/203279/appendix-b-huffman-tables-for-the-dc-and-ac-coefficients-of-the-jpeg-baseline-encoder
%   - http://stackoverflow.com/questions/6409233/matlab-conditional-assignment
%   - http://stackoverflow.com/questions/3024939/matrix-zigzag-reordering/3025478#3025478
%   - http://en.wikipedia.org/wiki/JPEG
%   - http://en.wikipedia.org/wiki/Run-length_encoding
%
% 2014
 
%
% Punto de entrada
%
function jpeg()
    
    inicializarTablas();
    nombreImagen = 'lena512color.tif';
    
    % Levanto la imagen en memoria
    imagen = imread(nombreImagen);
 
    % Convierto la imagen a escala de grises
    imagen = rgb2gray(imagen);
    
    % Muestro la imagen original como referencia
    subplot(1, 2, 1);
    imshow(imagen);
    title('Imagen original en escala de grises');
 
    % Obtengo dimensiones de la imagen
    infoImagen = imfinfo(nombreImagen);
    ancho = infoImagen.Width;
    alto = infoImagen.Height;
 
    % Convierto las dimensiones de la imagen en multiplos de 8
    [imagen, ancho, alto] = dimensionarMultiplo8(imagen, ancho, alto);
 
    % Obtengo un array de bloques de 8x8 de la imagen
    bloques = obtenerBloques(imagen, ancho, alto);
 
    % Llamo a la funcion codificarBloque por cada bloque de 8x8 de la
    % imagen. Obtengo un cell array de arrays de bits. Uno por cada bloque 
    % codificado
    bloquesCodificados = cellfun(@codificarBloque, bloques, 'UniformOutput', 0);
 
    % Obtengo la cantidad de bits que conforman cada bloque codificado
    mapaCompresion = cellfun(@numel, bloquesCodificados);
    
    % Calculo la cantidad de bits de la imagen comprimida, como suma de los
    % bits de cada bloque
    bitsImagenComprimida = sum(mapaCompresion(:));
    
    % Calculo el factor de compresión de toda la imagen. Se consideran
    % solamente los bits de datos de la imagen y no aquellos asociados a
    % metadata JPEG, que escapan a esta implementación
    factorDeCompresion = bitsImagenComprimida / (ancho * alto * 8);
    
    % Muestro el factor de compresión
    disp(['Factor de compresión: ', num2str(factorDeCompresion * 100), '%']);
 
    % Calculo el menor y mayor tamaño de bloque codificado. Serán los
    % valores correspondientes al negro y blanco en el mapa de compresión
    rangoDeCompresion = [min(mapaCompresion(:)) max(mapaCompresion(:))];
    
    % Muestro el mapa de compresión de la imagen. Representa la 
    % distribución del factor de compresión en el espacio. Como la 
    % compresión es por bloque, habra un píxel por cada bloque de la imagen
    % original. Se magnifica 8 veces para favorecer la comparación con la 
    % imagen original. Aquellas zonas mas oscuras de la imagen serán las 
    % que tienen un mayor grado de compresion. Se ven en blanco los
    % bordes de la imagen, que son aquellas zonas donde hay mayor contenido
    % de altas fecuencias y por lo tanto el nivel de compresion es menor
    subplot(1,2,2);
    imshow(mapaCompresion, rangoDeCompresion, 'InitialMagnification', 800);
    title('Mapa de compresión');
end
 
 
%
% Genera una imagen cuyas dimensiones son multiplo de 8 en base a la imagen
% recibida
%
function [imagenCompatible, anchoMultiplo, altoMultiplo] = dimensionarMultiplo8(imagen, ancho, alto)
 
    % Calculo el multiplo de 8 inmediato superior a ancho
    anchoEnBloques = ceil(ancho / 8);
    anchoMultiplo = anchoEnBloques * 8;
    
    % Calculo el multiplo de 8 inmediato superior a alto
    altoEnBloques = ceil(alto / 8);
    altoMultiplo = altoEnBloques * 8;
 
    % Creo una matriz de ceros de tamano multiplo de 8
    imagenCompatible = zeros(altoEnBloques * 8, anchoEnBloques * 8, 'uint8');
 
    % Copio los datos de la imagen a la matriz de las dimensiones multiplo 
    % de bloque
    imagenCompatible(1 : alto, 1 : ancho) = imagen(1 : alto, 1 : ancho);
end
 
 
%
% Genera un array de bloques de 8x8 en base a imagen
%
function bloques = obtenerBloques(imagen, columnas, filas)
 
    % Convierto el vector de datos de la imagen en una matriz de dimension 
    % filas x columnas
    matrizImagen = reshape(imagen, filas, columnas);
    
    % Calculo las dimensiones expresadas en bloques
    filasDeBloques = filas / 8;
    columnasDeBloques = columnas / 8;
    
    % Obtengo el array de bloques. mat2cell devuelve un array de celdas.
    % Las dimensiones de cada celda estan definidas en el segundo y tercer
    % parametros. En este caso siempre son de 8x8.
    bloques = mat2cell(matrizImagen, ...
        8 * ones(1, filasDeBloques), ...
        8 * ones(1, columnasDeBloques));
end
 
%
% Codifica un bloque de 8x8 y devuelve un vector de bits
%
function bits = codificarBloque(bloque)
    global lohscheller;
    
    % Ajusta el rango de luminancia de -128 a 127
    bloque = double(bloque) - 128;
    
    % Transformada coseno discreto en dos dimensiones
    cosenoDiscreto = dct2(bloque);
    
    % Aplica la matriz de cuantizacion
    cuantizado = round(cosenoDiscreto ./ lohscheller);
    
    % Recorre la matriz del bloque en zigzag y retorna un vector con los 
    % valores ordenados
    zigzag = leerEnZigZag(cuantizado);
    
    % Codifica los valores en RLE/VLC
    bits = longitudVariable(zigzag);
end
 
%
% Recorre la matriz en zigzag y devuelve un vector con los valores en ese
% orden
%
function vector = leerEnZigZag(matriz)
 
    % Genero un vector con las posiciones de la matriz en el orden en que 
    % se desea recorrer. Uso persistent para que se calcule solo la primera
    % vez que se llama a la funcion, luego se re-utiliza.
    % Referencia: http://stackoverflow.com/questions/3024939/matrix-zigzag-reordering/3025478#3025478
    persistent indices;
    if isempty(indices)
        indices = reshape(1:numel(matriz), size(matriz));
        indices = fliplr(spdiags(fliplr(indices)));
        indices(:, 1 : 2 : end) = flipud(indices(:,1:2:end));
        indices(indices == 0) = [];
    end
    
    % Recorro la matriz en zigzag y guardo los valores en un vector
    vector = matriz(indices);
end
 
%
% Codifica los coeficientes en RLE/VLC
%
function bits = longitudVariable(vector)
    global huffmanAC;
    global huffmanDC;
    
    function s = codificarEnComplementoUno(e)
        s = abs(de2bi(abs(e),'left-msb') - (e < 0));
    end
 
    % Obtengo el valor de continua
    valorDC = vector(1);
    
    % Codifico el valor en complemento a uno
    valorDC = codificarEnComplementoUno(valorDC);
    
    % Busco la cantidad de bits necesaria para representar el valorDC en
    % complemento a uno en la tabla de Huffman para DC y lo convierto de
    % acuerdo a ella
    sizeDC = huffmanDC(strcmp(numel(valorDC), huffmanDC), 2);
    
    % Junto las palabras tipo A y B del coeficiente de continua
    bitsDC = [cell2mat(sizeDC)' valorDC];
 
    % Obtengo las posiciones de los valores no nulos dentro del vector
    indicesNoNulos = find(vector);
    
    % En caso de haber coeficientes de alterna no nulos, los codifico
    if numel(indicesNoNulos) > 1
        
        % Obtengo un vector con los valores de los coeficientes no nulos
        valoresAC = vector(indicesNoNulos);
        
        % Remuevo el coeficiente de DC
        valoresAC(1) = [];
        
        % Diferenciando las posiciones de los valores no nulos obtengo la
        % cantidad de ceros anteriores a cada una de ellas. Por ejemplo, si
        % las posiciones fuesen [1 5 ...], la diferenciacion seria 
        % [0 4 ...] y la cantidad de ceros que preceden a la posicion 5 
        % serian 3
        run = diff(indicesNoNulos) - 1;
        
        % Calculo para cada coeficiente de AC la cantidad de bits
        % necesarios para su representacion en complemento a uno
        size = floor(log2(abs(valoresAC)) + 1);
        
        % Formo las palabras tipo A como dos caracteres hexadecimales,
        % donde el primero es la cantidad de ceros que lo preceden y el
        % segundo es la cantidad de bits necesarios para representarlo en
        % complemento a uno
        palabraTipoA = cellstr(strcat(dec2hex(run), dec2hex(size)));
        
        % Busco las palabras tipo A en la tabla de huffman y obtengo su
        % codificacion    
        palabraTipoA = cellfun(@(a) cell2mat(huffmanAC(strcmp(a, huffmanAC), 2)), palabraTipoA, 'UniformOutput', 0);
        
        palabraTipoB = arrayfun(@codificarEnComplementoUno, valoresAC, 'UniformOutput', 0);
 
        bitsAC = reshape([palabraTipoA'; palabraTipoB], 1, numel(palabraTipoA) + numel(palabraTipoB));
    else
        bitsAC = {};
    end
    
    % Agrego a los bits de AC el terminador de bloque EOB
    bitsAC = [bitsAC [1 0 1 0]];
    
    % Junto los bits resultantes de la codificacion del coeficiente DC y
    % los coeficientes AC
    bits = cell2mat([bitsDC bitsAC]);
end
 
% function decodificada = decodificar(bloque)
%     global lohscheller;
%     bloque = bloque .* lohscheller;
%     decodificada = uint8(idct2(bloque) + 128);
% end
 
%
% Inicializa los valores de las matrices utilizadas
%
function inicializarTablas()
    global lohscheller;
    global huffmanAC;
    global huffmanDC;
    
    lohscheller = [...
        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
    ];    
    huffmanDC = {...
        0 [0 0]
        1 [0 1 0]
        2 [0 1 1]
        3 [1 0 0]
        4 [1 0 1]
        5 [1 1 0]
        6 [1 1 1 0]
        7 [1 1 1 1 0]
        8 [1 1 1 1 1 0]
        9 [1 1 1 1 1 1 0]
        10 [1 1 1 1 1 1 1 0]
        11 [1 1 1 1 1 1 1 1 0]
    };
    huffmanAC = {...
        '00' [1 0 1 0]
        '01' [0 0]
        '02' [0 1]
        '03' [1 0 0]
        '04' [1 0 1 1]
        '05' [1 1 0 1 0]
        '06' [1 1 1 1 0 0 0]
        '07' [1 1 1 1 1 0 0 0]
        '08' [1 1 1 1 1 1 0 1 1 0]
        '09' [1 1 1 1 1 1 1 1 1 0 0 0 0 0 1 0]
        '0A' [1 1 1 1 1 1 1 1 1 0 0 0 0 0 1 1]
        '11' [1 1 0 0]
        '12' [1 1 0 1 1]
        '13' [1 1 1 1 0 0 1]
        '14' [1 1 1 1 1 0 1 1 0]
        '15' [1 1 1 1 1 1 1 0 1 1 0]
        '16' [1 1 1 1 1 1 1 1 1 0 0 0 0 1 0 0]
        '17' [1 1 1 1 1 1 1 1 1 0 0 0 0 1 0 1]
        '18' [1 1 1 1 1 1 1 1 1 0 0 0 0 1 1 0]
        '19' [1 1 1 1 1 1 1 1 1 0 0 0 0 1 1 1]
        '1A' [1 1 1 1 1 1 1 1 1 0 0 0 1 0 0 0]
        '21' [1 1 1 0 0]
        '22' [1 1 1 1 1 0 0 1]
        '23' [1 1 1 1 1 1 0 1 1 1]
        '24' [1 1 1 1  1 1 1 1 0 1 0 0]
        '25' [1 1 1 1 1 1 1 1 1 0 0 0 1 0 0 1]
        '26' [1 1 1 1 1 1 1 1 1 0 0 0 1 0 1 0]
        '27' [1 1 1 1 1 1 1 1 1 0 0 0 1 0 1 1]
        '28' [1 1 1 1 1 1 1 1 1 0 0 0 1 1 0 0]
        '29' [1 1 1 1 1 1 1 1 1 0 0 0 1 1 0 1]
        '2A' [1 1 1 1 1 1 1 1 1 0 0 0 1 1 1 0]
        '31' [1 1 1 0 1 0]
        '32' [1 1 1 1 1 0 1 1 1]
        '33' [1 1 1 1 1 1 1 1 0 1 0 1]
        '34' [1 1 1 1 1 1 1 1 1 0 0 0 1 1 1 1]
        '35' [1 1 1 1 1 1 1 1 1 0 0 1 0 0 0 0]
        '36' [1 1 1 1 1 1 1 1 1 0 0 1 0 0 0 1]
        '37' [1 1 1 1 1 1 1 1 1 0 0 1 0 0 1 0]
        '38' [1 1 1 1 1 1 1 1 1 0 0 1 0 0 1 1]
        '39' [1 1 1 1 1 1 1 1 1 0 0 1 0 1 0 0]
        '3A' [1 1 1 1 1 1 1 1 1 0 0 1 0 1 0 1]
        '41' [1 1 1 0 1 1]
        '42' [1 1 1 1 1 1 1 0 0 0]
        '43' [1 1 1 1 1 1 1 1 1 0 0 1 0 1 1 0]
        '44' [1 1 1 1 1 1 1 1 1 0 0 1 0 1 1 1]
        '45' [1 1 1 1 1 1 1 1 1 0 0 1 1 0 0 0]
        '46' [1 1 1 1 1 1 1 1 1 0 0 1 1 0 0 1]
        '47' [1 1 1 1 1 1 1 1 1 0 0 1 1 0 1 0]
        '48' [1 1 1 1 1 1 1 1 1 0 0 1 1 0 1 1]
        '49' [1 1 1 1 1 1 1 1 1 0 0 1 1 1 0 0]
        '4A' [1 1 1 1 1 1 1 1 1 0 0 1 1 1 0 1]
        '51' [1 1 1 1 0 1 0]
        '52' [1 1 1 1 1 1 1 0 1 1 1]
        '53' [1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 0]
        '54' [1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1]
        '55' [1 1 1 1 1 1 1 1 1 0 1 0 0 0 0 0]
        '56' [1 1 1 1 1 1 1 1 1 0 1 0 0 0 0 1]
        '57' [1 1 1 1 1 1 1 1 1 0 1 0 0 0 1 0]
        '58' [1 1 1 1 1 1 1 1 1 0 1 0 0 0 1 1]
        '59' [1 1 1 1 1 1 1 1 1 0 1 0 0 1 0 0]
        '5A' [1 1 1 1 1 1 1 1 1 0 1 0 0 1 0 1]
        '61' [1 1 1 1 0 1 1]
        '62' [1 1 1 1 1 1 1 1 0 1 1 0]
        '63' [1 1 1 1 1 1 1 1 1 0 1 0 0 1 1 0]
        '64' [1 1 1 1 1 1 1 1 1 0 1 0 0 1 1 1]
        '65' [1 1 1 1 1 1 1 1 1 0 1 0 1 0 0 0]
        '66' [1 1 1 1 1 1 1 1 1 0 1 0 1 0 0 1]
        '67' [1 1 1 1 1 1 1 1 1 0 1 0 1 0 1 0]
        '68' [1 1 1 1 1 1 1 1 1 0 1 0 1 0 1 1]
        '69' [1 1 1 1 1 1 1 1 1 0 1 0 1 1 0 0]
        '6A' [1 1 1 1 1 1 1 1 1 0 1 0 1 1 0 1]
        '71' [1 1 1 1 1 0 1 0]
        '72' [1 1 1 1 1 1 1 1 0 1 1 1]
        '73' [1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 0]
        '74' [1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1] 
        '75' [1 1 1 1 1 1 1 1 1 0 1 1 0 0 0 0]
        '76' [1 1 1 1 1 1 1 1 1 0 1 1 0 0 0 1] 
        '77' [1 1 1 1 1 1 1 1 1 0 1 1 0 0 1 0] 
        '78' [1 1 1 1 1 1 1 1 1 0 1 1 0 0 1 1] 
        '79' [1 1 1 1 1 1 1 1 1 0 1 1 0 1 0 0] 
        '7A' [1 1 1 1 1 1 1 1 1 0 1 1 0 1 0 1]
        '81' [1 1 1 1 1 1 0 0 0] 
        '82' [1 1 1 1 1 1 1 1 1 0 0 0 0 0 0] 
        '83' [1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 0] 
        '84' [1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 1] 
        '85' [1 1 1 1 1 1 1 1 1 0 1 1 1 0 0 0] 
        '86' [1 1 1 1 1 1 1 1 1 0 1 1 1 0 0 1] 
        '87' [1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 0] 
        '88' [1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 1] 
        '89' [1 1 1 1 1 1 1 1 1 0 1 1 1 1 0 0] 
        '8A' [1 1 1 1 1 1 1 1 1 0 1 1 1 1 0 1]
        '91' [1 1 1 1 1 1 0 0 1]
        '92' [1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0]
        '93' [1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1]
        '94' [1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0]
        '95' [1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 1]
        '96' [1 1 1 1 1 1 1 1 1 1 0 0 0 0 1 0]
        '97' [1 1 1 1 1 1 1 1 1 1 0 0 0 0 1 1]
        '98' [1 1 1 1 1 1 1 1 1 1 0 0 0 1 0 0]
        '99' [1 1 1 1 1 1 1 1 1 1 0 0 0 1 0 1]
        '9A' [1 1 1 1 1 1 1 1 1 1 0 0 0 1 1 0]
        'A1' [1 1 1 1 1 1 0 1 0]
    };
end