Transformada Discreta de Cosseno:
uma aplicação da Álgebra Linear na
compressão de imagens do formato JPEG

Resumo

Apresentamos como é feita a compressão de imagens com perdas como implementada no formato JPEG, um dos formatos de imagens mais usados na atualidade para fotos. Faremos em detalhe a Transformada Discreta de Cosseno, uma interessante aplicação da Álgebra Linear na computação, o núcleo do processo de compressão do JPEG.

Notação: vamos denotar o produto escalar de dois vetores \(v\) e \(u\) por \(\left\langle v,u \right\rangle\), por esta notação ser mais cômoda para somatórios.

Caso as fórmulas abaixo não estejam sendo corretamente renderizadas, opte por ler o texto em pdf.

1. Passos da compressão JPEG

• A representação de cores é convertida de \(RGB\) (usada em monitores) para \(YC_bC_r\), consistindo em uma componente de brilho \(Y\) e duas componentes de cores \(C_b\) e \(C_r\) (\(b\) e \(r\) de blue e red). Neste passo aparece por primeira vez a conversão por transformações lineares, mas que não trataremos aqui. Este passo às vezes é omitido.
• A resolução dos dados de cores é reduzida, usualmente, pela metade. Isto reflete o fato de que o olho é menos sensível a detalhes finos de cor do que de brilho.

• A imagem é dividida em blocos de \(8\times8\) pixels, e a cada bloco é aplicada a Transformada Discreta de Cosseno com sigla DCT, do inglês. Este é o passo mais importante e será explicado em detalhe. Usa fortemente as bases ortogonais da Álgebra Linear.

  Na terminologia da transformada, o espaço dos dados iniciais são chamados de dados espaciais e os resultados após a transformada, em dados de frequência.

• A visão humana é muito mais sensível a pequenas variações de cor e brilho em áreas maiores do que à variações de brilho em altas frequências. Assim, as amplitudes dos componentes de frequência são “quantizadas”, isto é, reescaladas usando um sistema de pesos, onde pesos maiores são aplicados a frequências menores. Portanto a magnitude dos componentes de alta frequência são guardadas com precisão menor que componentes de baixas frequências. Aqui entra a qualidade da compactação. A configuração de qualidade do codificador (por exemplo 50% ou 95% numa escala de 0 a 100%) afeta a extensão com a qual cada componente de frequência é reduzida. Se uma configuração de baixa qualidade é usada, os componentes de alta frequência são todos descartados.

  A matriz de quantização é gravada no início do arquivo JPEG para posterior desquantização quando a imagem for lida por outra aplicação.

• Os dados resultantes para todos os blocos \(8\times8\), que tipicamente são matrizes onde os números significativos estão no canto superior esquerdo, e com aspecto de matriz esparsa, são novamente compactados com um algoritmo sem perdas, uma variante da codificação de Huffman e armazenados no arquivo JPEG.

2. Transformada Discreta de Cosseno (DTC)

A DCT foi desenvolvida em 1974 por N. Ahmed, T. Natarajan and K. R. Rao [1] e tornou-se, por sua eficiência em ser calculada e pela excelente compressão que proporciona, uma das ferramentas mais usadas em processamento de imagens; ela é semelhante à transformada de Fourier (para funções periódicas) no sentido de que ela produz um tipo de espectro de espaço-frequência. Existem várias variantes da DTC --- aqui estudaremos a DTC-II, mais usada.

Espera-se que a maioria das amostras de dados reais (por exemplo uma linha de cores de pixels em uma imagem) tenham valores com uma transição quase contínua de uma entrada para sua vizinha. O funcionamento da DCT está na escolha de uma base onde o primeiros elementos tem entradas vizinhas com pouca variação de valores (isto é, tem baixa frequência de variação) e os últimos tem entradas vizinhas com grandes variações (alta frequência). Assim, ao escrevermos os dados de amostra como combinação linear da base, é de se esperar que os coeficientes de maior magnitude sejam os primeiros, e que os últimos coeficientes sejam praticamente nulos, podendo ser desprezíveis.

Na linguagem da Álgebra linear, o que ocorre é uma mudança de bases da base canônica para outra base, ortonormal, em que cada elemento desta base representa a discretização do cosseno de uma determinada frequência.

A partir de agora, fixemos \(N\in\mathbb N\), que fixa a dimensão do espaço de amostragem.

2.1 DTC unidimensional

Começaremos com a DTC unidimensional. Dado \(n\in \{0,1,\dots,N-1\}\), à frequência \(n\pi/N\), associamos o vetor coluna que interpola a função \(\cos\bigl(\frac{n\pi}{N}t\bigr)\), dividido por sua norma. Mais precisamente, seja o vetor \begin{equation}\label{eq:v_n} v_n = \left( \cos\Bigl( \frac{n\pi}{N}\frac12 \Bigr) , \cos\Bigl( \frac{n\pi}{N}\Bigl( 1+\frac12 \Bigr) \Bigr) , \dots, \cos\Bigl( \frac{n\pi}{N}\Bigl( N-1+\frac12 \Bigr) \Bigr) \right)^T. \end{equation}

Na Figura 1 vemos as entradas de cada vetor \(v_n\) para o caso \(N=8\).

Figura 1: os vetores ortogonais \(v_n\) para o caso \(N=8\), obtidos pela interpolação dos cossenos nas diversas frequências.

Queremos mostrar que \(C=\{v_0,v_1,\dots,v_{N-1}\}\) é uma base ortogonal. Para isto, basta mostrar que \(\left\langle v_n,v_m \right\rangle=0\). Antes, precisamos de um resultado auxiliar.

Prova: O caso \(n=0\) é trivial. Suponha \(n\neq0\). Recordando as identidades \begin{gather*} \cos \theta = \operatorname{Re} e^{i\theta},\\ \operatorname{Re} \bar w = \operatorname{Re} w \implies \operatorname{Re} (\bar w - w) = 0, \quad \forall w\in \mathbb C,\\ z^N-1 = (z-1)(1+z+z^2+ \dots + z^{N-1}), \quad \forall z\in \mathbb C, \end{gather*} e definindo os números complexos \(w = \exp\bigl(\frac{in\pi}{2N}\bigr)\) e \(z = \exp\bigl(\frac{in\pi}{N}\bigr) = w^2\), de forma que \(z^N = e^{iN\pi} = (-1)^N\) e \(|w|=1\), obtemos \begin{align*} \sum_{k=0}^{N-1} \cos\biggl( \frac{n\pi}{N}\biggl( k+\frac12 \biggr) \biggr) & = \operatorname{Re} \sum_{k=0}^{N-1} \exp\biggl( \frac{in\pi}{N}\biggl( k+\frac12 \biggr) \biggr) \\ & = \operatorname{Re} \sum_{k=0}^{N-1} \underbrace{\exp\biggl( \frac{in\pi}{2N} \biggr) }_w \underbrace{\exp\biggl( \frac{in\pi }{N}k \biggr) }_{z^k} \\ & = \operatorname{Re} w(1+z+z^2+\dots+z^{N-1}) = \operatorname{Re} w \frac{z^n-1}{z-1} \\ & = \operatorname{Re} w \frac{(z^N-1)}{(z-1)}\frac{(\bar z-1)}{(\bar z-1)} = \frac{(-1)^N-1}{|z-1|^2} \operatorname{Re} w (\bar z-w) \\ & = \frac{(-1)^N-1}{|z-1|^2} \operatorname{Re} (\underbrace{w \bar w}_{|w|^2 =1} \bar w -w) = \frac{(-1)^N-1}{|z-1|^2} \operatorname{Re} (\bar w - w) = 0. \end{align*} c.q.d.

Prova: Sejam \(n,m \in{0,1,\dots,N-1}\). Podemos supor que \(n\geqslant m\). Recordemos a identidade \[ \cos a\cos b = \frac12\bigl(\cos(a+b)+\cos(a-b)\bigr). \]

Então \begin{align} \left\langle v_n,v_m \right\rangle & = \sum_{k=0}^{N-1} \cos\biggl( \frac{n\pi}{N}\biggl( k+\frac12 \biggr) \biggr) \cos\biggl( \frac{m\pi}{N}\Bigl( k+\frac12 \biggr) \Bigr) \nonumber\\ & = \frac12 \sum_{k=0}^{N-1} \cos\biggl( \frac{\pi}{N}\biggl( k+\frac12 \biggr) (n+m) \biggr) + \frac12 \sum_{k=0}^{N-1} \cos\biggl( \frac{\pi}{N}\Bigl( k+\frac12 \Bigr)(n-m) \biggr) \label{eq:v_n.v_m} \end{align}

Se \(n=m=0\), então \(n+m=n-m=0\) e portanto de \eqref{eq:v_n.v_m} e do Lema 1, obtemos \(\left\langle v_0,v_0 \right\rangle=N\). Se \(n=m\neq0\), então \(n+m=2n\neq0\) mas \(n-m=0\) e portanto de \eqref{eq:v_n.v_m} e do Lema 1, obtemos \(\left\langle v_0,v_0 \right\rangle=N/2\). Se \(n\neq m\), então \(n+m\neq 0\) e \(n-m\neq0\), de onde vem \(\left\langle v_n,v_m \right\rangle=0\). c.q.d.

Seja \(V=\mathbb R^N\) e \(x\in V\). Como \(B\) é base, existem escalares \(X_0, X_1,\dots,X_{N-1}\) tais que \[ x = \sum_{k=0}^{N-1} X_k u_k. \] Podemos calcular \(x_k\) pelo produto interno \begin{equation} \label{eq: x_n} \left\langle x,u_n \right\rangle = \left\langle \sum_{k=0}^{N-1} X_k u_k, u_n \right\rangle = \sum_{k=0}^{N-1} X_k \underbrace{\left\langle u_k,u_n \right\rangle}_{\begin{smallmatrix} 1,& k=n\\ 0,& k\neq n \end{smallmatrix}} = X_n. \end{equation}

Definimos então a DCT pela mudança de base da base canônica para a base \(C\), ou seja, a DCT do vetor \(x = [x_0,x_1,\dots,x_{N-1}]\in\mathbb R^N\) é o vetor \(DTC(x) = [X_0, X_1, \dots, X_{N-1}]\). Nas referências sobre DCT em artigos com enfoque computacional, a fórmula para os componentes de frequência \(X_k\) é dada com o produto interno expandido explicitamente como em \[ X_n = \frac1{\|v_n\|} \sum_{k=0}^{N-1} x_k \cos\Bigl( \frac{n\pi}{N}\Bigl( k+\frac12 \Bigr) \Bigr) , \quad n>0, \] onde \(\|v_0\|= \sqrt N\) e \(\|v_n\|=\sqrt{N/2}\), \(n>0\).

2.2 DCT bidimensional

Caso os dados sejam mais convenientemente representados em forma matricial, podemos definir a DTC bidimensional compondo a DCT unidimensional nas linhas e depois nas colunas. Por exemplo, na DCT usada no padrão de imagens JPEG, usa-se a DTC bidimensional dada pela base ortogonal \begin{equation}\label{eq:DCT 2d} \begin{gathered} C_2 = \{V_{n,m} \in M_{8\times 8}(\mathbb R): 0\leqslant n,m \lt 8\},\\ V_{n,m} = (a_{n,m}(i,j))\quad \text{ com} \\ a_{n,m}(i,j) = \cos\Bigl( \frac{n\pi}{8}\Bigl( i+\frac12 \Bigr) \Bigr) \cos\Bigl( \frac{m\pi}{8}\Bigl( k+\frac12 \Bigr) \Bigr) \end{gathered} \end{equation} onde a norma \(\|V\|\) e o produto interno \(\left\langle U,V \right\rangle\) para \(U,V\in M_{8\times8}(\mathbb R)\), \(U=(a_{i,j})\) e \(V=(b_{i,j})\), são os usuais, dados por \[ \left\langle U,V \right\rangle = \sum_{i=0}^7 \sum_{j=0}^7 a_{i,j}b_{i,j},\qquad \|U\| = \sqrt{\sum_{i=0}^7 \sum_{j=0}^7 (a_{i,j})^2} = \sqrt{\left\langle U,U \right\rangle}. \]

Na Figura 2 vemos cada uma das matrizes \(v_{n,m}\) de \(C_2\) representada como uma imagem de tamanho \(8\times8\).

Figura 2: A base ortogonal dos 64 vetores \(v_{n,m}\) da DTC bidimensional usadas no JPEG.

Também a base \(C_2\) é ortogonal.

Prova: \begin{align} \left\langle V_{n,m} , V_{k,\ell} \right\rangle & = \sum_{i=0}^7\sum_{j=0}^7 \cos\biggl( \frac{n\pi}{8}\biggl( i+\frac12 \biggr) \biggr) \cos\biggl( \frac{m\pi}{8}\Bigl( j+\frac12 \Bigr) \biggr) \cos\biggl( \frac{k\pi}{8}\biggl( i+\frac12 \biggr) \biggr) \cos\biggl( \frac{\ell\pi}{8}\Bigl( j+\frac12 \Bigr) \biggr) \nonumber\\ & = \left[ \sum_{i=0}^7 \cos\biggl( \frac{n\pi}{8}\biggl( i+\frac12 \biggr) \biggr) \cos\biggl( \frac{k\pi}{8}\biggl( i+\frac12 \biggr) \biggr) \right] \times \left[\sum_{j=0}^7 \cos\biggl( \frac{m\pi}{8}\Bigl( j+\frac12 \Bigr) \biggr) \cos\biggl( \frac{\ell\pi}{8}\Bigl( j+\frac12 \Bigr) \biggr) \right]. \label{eq:prod-2d} \end{align}

Cada fator em \eqref{eq:prod-2d} é um somatório idêntico ao de \eqref{eq:v_n.v_m}, e as discussões seguindo \eqref{eq:v_n.v_m} implicam o que está no enunciado. c.q.d.

3. Um exemplo

Na Figura 3 (abaixo) vemos selecionada uma região da imagem de \(8\times 8\) pixels. As escalas de cinza (de 0 a 255) estão mostradas na matriz \(A\) abaixo. \[ A= \begin{pmatrix} 89& 97& 101& 101& 113& 206& 231& 235\\ 85& 85& 89& 89& 109& 231& 235& 235\\ 85& 81& 85& 89& 105& 215& 239& 239\\ 105& 101& 101& 101& 113& 239& 243& 243\\ 178& 186& 190& 186& 194& 231& 239& 243\\ 235& 239& 239& 243& 243& 239& 243& 243\\ 235& 239& 235& 235& 231& 239& 231& 239\\ 231& 231& 231& 231& 227& 227& 227& 227 \end{pmatrix} \]

Figura 3: Imagem original (com um quadrado selecionado para aplicarmos o exemplo) e abaixo os resultados após a compressão com matriz de quantização iguais respectivamente a \(\frac12 Q\) e \(Q\), dada em \eqref{eq:Q}. Note que a segunda é praticamente igual à original. A terceira tem certos detalhes perdidos, como os parafusos nas vigas cinzas (sombras) embaixo da ponte, mas ainda assim com qualidade satisfatória.

Aplicando a DCT a \(A\), obtemos \(B=\operatorname{DCT}(A)=(b_{i,j})\) onde \[ b_{i,j} = \frac{\langle A, V_{i,j}\rangle}{\|V_{i,j}\|} \] ou em números \[ B= \left( \begin{array}{rrrrrrrr} 1500.2 & -263.7 & 91.0 & 26.8 & -61.0 & 16.1 & 28.3 & -32.3 \\ -311.9 & -211.7 & 71.8 & 25.3 & -47.2 & 15.4 & 20.8 & -20.3 \\ 16.9 & 38.3 & -13.0 & -8.6 & 9.5 & -2.8 & -5.7 & 4.4 \\ 103.8 & 62.4 & -27.1 & -17.6 & 15.7 & -8.7 & -11.4 & 10.2 \\ -25.7 & -8.3 & 5.9 & -3.2 & -5.0 & 1.2 & 0.7 & -0.2 \\ -24.8 & -16.4 & 13.4 & 1.6 & -4.0 & 2.4 & 2.6 & -6.3 \\ 10.8 & 11.6 & -4.5 & -8.5 & 9.3 & 3.7 & -10.2 & 8.0 \\ 2.5 & 10.6 & -2.7 & -11.0 & 6.0 & -2.4 & -5.2 & 5.9 \end{array} \right) \]

A esta matriz aplicaremos uma quantização por uma matriz típica \(Q=(q_{i,j})\) do padrão JPEG. Esta matriz de quantização para luminância é dada em [2, p. 143]. \begin{equation}\label{eq:Q} Q= \begin{pmatrix} 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 \end{pmatrix}. \end{equation}

A matriz \(C=(c_{i,j})\), gerada pela quantização de \(B=(b_{i,j})\), é dada por \[ c_{i,j} = \operatorname{round} \left(\frac{b_{i,j}}{q_{i,j}}\right) \] e aparece abaixo: \[ C = \begin{pmatrix} 94 & -24 & 9 & 2 & -3 & 0 & 1 & -1 \\ -26 & -18 & 5 & 1 & -2 & 0 & 0 & 0 \\ 1 & 3 & -1 & 0 & 0 & 0 & 0 & 0 \\ 7 & 4 & -1 & -1 & 0 & 0 & 0 & 0 \\ -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix} \]

Esta matriz esparsa será compactada com o algoritmo de Huffman (sem perdas) e armazenada no arquivo JPEG.

Para vermos a eficácia da compressão, aplicaremos as operações inversas para reconstruir uma aproximação de \(A\) (mantendo os valores em \([0,255]\)) e mostramos abaixo o retângulo \(8\times 8\) original e, ao lado, o produzido pela compressão descrita acima:

\(\qquad\)

A. Código fonte de tudo e créditos

Baixe o arquivo jpeg.zip que contém todos os arquivos necessários para gerar todos os “subprodutos” deste projeto, isto é, o pdf e o html do texto e os cálculos numéricos.

• Para gerar o PDF, foi usada a linguagem \(\rm{\LaTeX}\) (compilada com pdflatex). O arquivo é o jpeg.tex. Para gerar as figuras de alta qualidade, usamos o MetaPost, que tem como saída os arquivos vetoriais com extensão mps (eps que pode ser usado com pdflatex).
• O html foi convertido “no braço” a partir do arquivo \(\rm{\LaTeX}\), tendo em vista a possibilidade de usarmos MathJax para renderizar as fórmulas.
• Os cálculos numéricos foram feitos usando a linguagem GNU Octave, semelhante ao MatLab, mas aquela é software livre. O arquivo jpeg.m lê o arquivo ponte.pgm e gera diversos arquivos pgm (Portable Graymap), um formato de imagem fácil de programar por ser texto. Posteriormente, os arquivos foram convertidos para png com o editor de imagens Gimp.
• A imagem do riacho e ponte usada no exemplo (Figura 3) teve sua origem no antigo banco de dados de imagens teste da USC-SIPI. O arquivo original é a imagem 5.2.10, na seção Miscellaneous.

Referências

[1] Ahmed, N., Natarajan, T., Rao, K. R. Discrete Cosine Transform. IEEE Transactions on Computers C-23 (1): 90-93. 1974.

[2] ITU-T.81 – Recomendação T.81 do International Telecommunication Union (ITU) que define o formato JPEG, 1992.

[3] Wikipedia – verbete JPEG, em inglês.

[4] \(\rm{\LaTeX}\) – a document preparation system.

[5] The MetaPost system – a picture-drawing language.

[6] MathJax, an JavaScript display engine for mathematics.

[7] GNU Octave – a high-level interpreted language, primarily intended for numerical computations.

[8] Wikipedia – verbete Portable Graymap (pgm), em inglês.

[9] Gimp, the GNU Image Manipulation Program.