Jak funguje JPEG - princip pod lupou Zajímalo vás někdy, jak napsat kompresní program na obrázky založený na stejném principu jako JPEG? Rád vám vysvětlím, jak to vlastně funguje a jak jednoduše se dá základ takového kompresoru napsat. Na začátku 80. let se na scéně objevil standard JPEG (Joint Photographic Experts Group) pojmenovný podle skupiny, kterou byl vytvořen. Byl to zároveň první krok pro ztrátové kompresních metody v grafice, který však přetrvává dodnes. Formát .JPG už zná prakticky každý, kdo vlastní digitální fotoaparát nebo serfuje po webu. O JPEGu (vyslovuj jay-pegu) bylo napsáno tolik, že už snad ani nemá cenu omílat to všechno dokola. Je ale nutno dodat, že JPEG zůstal populární hlavně kvůli své velké flexibilitě. Diskrétní kosinová transformace Formát JPEG je založen na tzv. DCT - diskrétní kosinové transformaci, která je realizována ve dvou rozměrech. Aplikuje se na čtvercové matice o velikosti 8x8 pixelů, obrázek se tedy rozřeže na čtverečky, které se potom zpracovávají jeden po druhém. Je ale nutno podotknout, že DCT sama o sobě není kompresní metoda. Jde pouze o transformaci, která je popsána jednou rovnicí a má i svou inverzní variantu. Pokud na matici aplikujete dopřednou DCT a potom inverzní DCT, nic se na ní nezmění. Původní pixelová matice sestává pouze z hodnot pixelů (pro černobílé obrázky to je stupeň šedi, pro barevné hodnota jedné složky), kdežto DCT matice obsahuje tzv. DC koeficienty tvořící jakousi frekvenční mapu. Při kompresi se využívá toho, že lidské oko není citlivé na funkce vysokých kmitočtů. Amplitudy těchto funkcí jsou soustředěné v jednom rohu DCT matice a naopak v protějším rohu jsou nízké frekvence (pomalejší barevné přechody). Pro kompresi se potom používá tzv. kvantizační matice, která má největší hodnoty právě v místech, kde jsou vysoké frekvence. Pomocí jednoduchého dělení matic se odřežou vysoké frekvence a v DCT matici tak vznikne velké množství nul. To je redundantní obsah, který se pak zredukuje nějakou neztrátovou metodou (obvykle Huffmanovo kódování). Při dekompresi se zredukovaná DCT matice opět vynásobí kvantizační maticí a pak se aplikuje I-DCT. Tím dostaneme aproximaci původního čtverečku (bloku). Míra komprese a tím taky kvalita výstupu je dána faktorem kvality (nazývá se Q-Faktor), který v podstatě stupňuje hodnoty v kvantizační matici. Funkce vyšších frekvencí obvykle nejsou tak viditelné, ovšem ostré hrany jsou často výjimkou. Tady potom dochází k nepříjemnému zvlnění (Gibbsův jev) a vůbec lehkému rozostření celého obrázku. Proto není JPEG ani žádná jiná metoda založená na DCT vhodná pro kompresi technických výkresů a prezentačních materiálů, naopak je výborná pro zpracování digitálních fotografií a videa (DCT využívají všechny MPEG kodeky). Ve specifikaci JPEG se samozřejmě mluví o mnoha jiných technologiích vedle kosinové transformace samotné, o kterých se zmíním, nicméně teď přejdu k tomu, jak lze DCT realizovat. Abychom si moc nekomplikovali život, budeme pracovat s černobílým obrázek v 256 stupních šedi (já jsem použil ovladač SuperVGA, takže se ten obrázek bude tvářit, že má 256 odstínů, ve skutečnosti bude mít 64, ale i to bohatě stačí). Pro začátek trocha teorie. Nebudu se pouštět do matematických rozborů, od toho jsou příslušné publikace. Vězme, že naším prvotním cílem je získat z pixelové matice o velikosti 8x8 onu frekvenční mapu, čili stejně velkou matici DC koeficientů. Pozice v pixelové matici se bude udávat souřadnicemi x,y; v DCT matici potom i,j. Vztah pro výpočet DC koeficientu DCT matice na pozici [i,j] je dán vztahem:
, kde
pro u=0, jinak C(u)=1 pro u>0
Malé n je rozměr čtvercového blokuv pixelech. V programu tenhle vztah interpretujete až překvapivě jednoduše. Celý vztah bude fungovat pro n=8. Sumy se samozřejmě realizují pomocí vnořeného for cyklu. Já jsem to řešil tak, že jsem nejprve posečítal ty "šílené" výrazy napravo a teprve potom celek vynásobil. Ostatně se podívejte sami: { DCT } for j:=0 to 7 do for i:=0 to 7 do begin DCT[i,j]:=0; for y:=0 to 7 do for x:=0 to 7 do DCT[i,j]:=DCT[i,j]+ ( A[x,y]*cos( ((2*x+1)*i*pi)/(2*8) )*cos( ((2*y+1)*j*pi)/(2*8) ) ); DCT[i,j]:=DCT[i,j]*C(j)*C(i)*( 1/sqrt(2*8) ); end;
Takže jak vidíte, není třeba hledat žádné složitosti. Inverzní vztah, který v podstatě provádí operace v opačném pořadí, je I-DCT. Připomínám, že se bavíme jen o dvourozměrné variantě:
pro C(u) platí totéž, co je uvedeno výše. Tímto vztahem dostanete z DCT matice zpátky hodnoty pixelů. První program Nyní už k realizaci komprese. První program provádí část JPEG komprese, přičemž vychází z výše uvedených vztahů. Komprese se totiž dá napsat tak, aby běžela mnohem rychleji. Jednak se to provádí vytvořením pomocné kosinové tabulky a potom také úpravou algoritmu tak, že procesor pracuje většinou v celých číslech. Komerční programy tak dosahují radikálně rychlejší (de)komprese za cenu triviální ztráty kvality. Následující kód používá černobílý ukázkový obrázek (níže) a ovladač svga256.bgi:
program DCT256; { image compression using DCT uses Crt,Graph; const kvalita = 15; { <1,+-25> } type Matice = array[0..7,0..7] of real; var d,m, xx,yy,x,y,i,j: integer; bmp: file of byte; b: byte; A,DCT,QM: Matice; procedure Pal(n,R,G,B: byte); begin Port[$3C8]:=n; Port[$3C9]:=R; Port[$3C9]:=G; Port[$3C9]:=B; end; function C(u: real): real; begin if u=0 then C:=1/sqrt(2) else C:=1;
by Libor Tinka }
end; begin { inicializace grafiky } d:=InstallUserDriver('svga256',nil);m:=0;InitGraph(d,m,'svga256.bgi'); { inicializace palety } for x:=0 to 255 do Pal(x,Round(x*63/255),Round(x*63/255),Round(x*63/255)); { nacteni obrazku } assign(bmp,'exmpl256.bmp');reset(bmp);Seek(bmp,1078); for y:=0 to 199 do for x:=0 to 319 do begin read(bmp,b); PutPixel(x,199-y,b); end; Close(bmp); { *** KOMPRESE *** } { sestaveni kvantizacni matice } for y:=0 to 7 do for x:=0 to 7 do QM[x,y]:=1+((1+x+y)*kvalita); { prochazeni bloku } for yy:=0 to 24 do for xx:=0 to 39 do begin { nacteni bodu do pixelove matice } for y:=0 to 7 do for x:=0 to 7 do A[x,y]:=GetPixel(xx*8+x,yy*8+y); { snizeni hodnot } for y:=0 to 7 do for x:=0 to 7 do A[x,y]:=A[x,y]-128; { DCT } for j:=0 to 7 do for i:=0 to 7 do begin DCT[i,j]:=0; for y:=0 to 7 do for x:=0 to 7 do DCT[i,j]:=DCT[i,j]+ ( A[x,y]*cos( ((2*x+1)*i*pi)/(2*8) )*cos( ((2*y+1)*j*pi)/(2*8) ) ); DCT[i,j]:=DCT[i,j]*C(j)*C(i)*( 1/sqrt(2*8) ); end; { kvantizace } for j:=0 to 7 do for i:=0 to 7 do DCT[i,j]:=Round(DCT[i,j]/QM[i,j]); { dekvantizace } for j:=0 to 7 do for i:=0 to 7 do DCT[i,j]:=DCT[i,j]*QM[i,j]; { iDCT } for y:=0 to 7 do for x:=0 to 7 do begin A[x,y]:=0; for i:=0 to 7 do for j:=0 to 7 do A[x,y]:=A[x,y]+ ( C(i)*C(j)*DCT[i,j]*cos( (2*x+1)*i*pi/(2*8) )*cos( (2*y+1)*j*pi/(2*8) ) ); A[x,y]:=A[x,y]*( 1/sqrt(2*8) ); end; { zvyseni hodnot } for y:=0 to 7 do for x:=0 to 7 do begin A[x,y]:=A[x,y]+128; if A[x,y]>255 then A[x,y]:=255; if A[x,y]<0 then A[x,y]:=0; end; { vykresleni } for y:=0 to 7 do for x:=0 to 7 do PutPixel(xx*8+x,yy*8+y,Round(A[x,y])); if KeyPressed then begin CloseGraph;Halt; end; end; readkey; CloseGraph; end.
Snažil jsem se to psát co možná nejpřehledněji. Jendotlivé kroky jsou odděleny komentáři, takže pěkně popořádku - po inicializaci grafiky, palety a načtení testovacího obrázku se připraví kvantizační matice. Ta bude později ovlivňovat poměr kvalita/komprese. Výslednou kvalitu si můžete sami regulovat tím, že nastavíte konstantu kvalita na nějakou hodnotu od 1 do 25 (můžete dát i víc, ale rozhodně ne míň - dělilo by se nulou). Aby jste viděli, jak taková kvantizační matice vypadá, takováto vznikne použitým vztahem při nastaveném faktoru kvality 2:
Kvantizace Při kompresi se budou dělit DC koeficienty odpovídajícími hodnotami v kvantizační matici. Největší ztrátu utrpí DC koeficienty v pravém spodním rohu matice, protože budou děleny největším číslem. Tam se také nacházejí vyšší frekvence. Uvedená kvantizační matice je ideální a nevychází tak úplně z charakteristiky lidského zraku. Skupina JPEG definovala jakousi optimální kvantizační matici, která je postavena na anatomii našeho oka a měla by znamenat rovnováhu mezi kvalitou a mírou komprese:
Hodnota 50 tu vystupuje jako střední. Určitě jste si sami všimli, že v nastavení formátu JPEG volíte obvykle kvalitu na stupnici od 0-100. Jestliže chcete využívat JPEG kvantizační matici, tak stačí vycházet z této, přičemž všechny ostatní z ní snadno odvodíte. Jestliže kvantizační faktor q nastavíte od 50 do 100, potom pro získání kvantizační matice upravte všechny hodnoty té vzorové takto:
Q-faktor nad 50 značí vyšší kvalitu obrazu, ale nižší kompresní poměr. Naopak, posunete-li jezdec q-faktoru pod padesát, upraví se kvantizační matice tímto způsobem:
Komprese a dekomprese Ale vraťme se k našemu programu. Kvantizační tabulka je sestavena a nyní jdeme vybírat jednotlivé bloky. Je nutno poznamenat, že tyto bloky, někdy označované jako makrobloky, mohou mít velikost i 16x16, 32x32... ale pixelové matice mohou být také obdélníkové. Údajně se při blocích 16x16 dosáhne trochu vyšší komprese, ale kosinová transformace je potom také zdlouhavější na výpočet. Asi proto
komise JPEG ustavila bloky 8x8 jako standard. Blok je načten pomocí funkce GetPixel do pole A, které představuje pixelovou matici. Jsou v ní hodnoty v rozmezí <0,255>, což se nehodí pro DCT. Hodnoty je nutno přesunout do intervalu <-128,127>, což se provede pouhým odečtením 128 od hodnoty pixelu. Tady samozřejmě k žádné ztrátě nedochází. Následuje vlastní DCT, prvními dvěma cykly pro hodnoty i a j se pohybujeme v DCT matici a počítáme její koeficienty na pozici [i,j] pomocí dalších dvou vnořených cyklů. To už je doslovná realizace vztahu pro výpočet DC koeficientů. Kdybychom nyní provedli inverzní DCT a získané hodnoty posunuli zpátky o 128 nahoru, dostali bychom přesně originál. Příkladně samotné odečetní je také transformací (lineární) a má svou inverzi (přičítání). Mezi DCT a I-DCT je tedy třeba vložit onen ztrátový krok a tím je kvantizace. Nyní už jde o dělení, ovšem se zaokrouhlením na celá čísla. Musíme uvažovat, že koeficienty bude kompresní utilita ukládat do souboru, kde by desetinná čísla zabrala více místa, než měl samotný originál - proto zaokrouhlení. { kvantizace } for j:=0 to 7 do for i:=0 to 7 do DCT[i,j]:=Round(DCT[i,j]/QM[i,j]); { dekvantizace } for j:=0 to 7 do for i:=0 to 7 do DCT[i,j]:=DCT[i,j]*QM[i,j];
Po kvantizaci už máme kompresi hotovou a hned se přehoupneme k dekompresi. Všechny kroky se provedou v opačném pořadí. Dekvantizace se provede opět roznásobením kompenzovaných DC hodnot hodnotami kvantizační matice. Obvykle se používá jedna kvantizační matice pro celý obrázek, ovšem v JPEG to je pravidlem. Následuje pochopitelně I-DCT a potom posunutí hodnot zpátky na původní hladinu. Po ztrátovém kroku nám nyní občas některá hodnota "odskočí" z chtěného intervalu <0,255>, ale tím si nemusíme lámat hlavu. Dekomprimovaný blok se znovu vykreslí na obrazovku. Takto se pokračuje, až je celý obrázek zpracovaný. A co dál? JPEG specifikace však nehovoří jenom o těchto několika krocích. Truecolorové obrázky jsou obvykle vyjadřovány pomocí RGB, tedy složek červené, modré a zelené. U JPEG se jedná vždy o 24-bitovou grafiku, takže na každou barevnou šložku připadá 8 bitů, tedy 256 odstínů. Model RGB však není zrovna nejvhodnější pro ztrátový formát. Hodí se spíš tam, kde jsou zapotřebí exaktní barvy. Člověk vnímá barvy díky zvláštním barvocitným receptorům na sítnici. Jsou dvou druhů, přičemž jeden je citlivý na modré a druhý na červené světlo. Valná většina ostatních receptorů příjimá pouze luminanci, tedy intenzitu bílého světla. Na základě tohoto poznatku vznikl barevný model YUV (často také označovaný jako YCbCr), používaný mimo jiné i v televizních normách. U komprese údajů jde o to, že oko je podstatně citlivější na změny světlosti (Y), než na barevné odchylky (U,V). Barevný obrázek je rozložený na tři barevné roviny a každá se komprimuje zvlášť. Jediný rozdíl je v tom, že roviny U a V se zredukují co do velikosti podvzorkují se v poměru 1:2. My si toho příliš nepovšimneme, ale výrazně se tak zvýší kompresní poměr. Další kroky jsou podobné těm, co jsem ukázal v programu, jenom se provádějí pro každou barevnou rovinu zvlášť. Je však zapotřebí uložit kvantifikované DC koeficienty bez redundance (nadbytečnosti). Je tu spousta nul po vydělení kvantizační maticí. DCT matice se převede na jednorozměrnou sekvenci způsobem, kterému se výstižně říká Zig-zag. V matici se pohybujeme po diagonálách a jednotlivé hodnoty vypisujeme do posloupnosti. Vysldkem je to, že dostaneme vedle sebe velké množství nul, takže redundace této posloupnosti je maximální. Obvykle se potom aplikuje jednoduché proudové kódování (RLE - Run Length Encoding) a potom ještě Huffmanovo kódování. Pro dosažení ještě o něco lepších výsledků je možné použít aritmetické kódování, ale zde jsou potřeba operace s čísli s dlouhým dekadickým rozvojem. To je nepraktické co do výkonu, a tudíž není v JFIF (souborový formát JPEG) aplikováno (ale asi také proto, že přišlo nějaký čas po standardizaci :). Jednoduše rychlá DCT Vypracoval jsem ještě mnohem svižnější verzi kompresního programu, která šetří procesor od otrockého počítání kosinů se stále stejným argumentem. Tentokrát se vytvoří substituční kosinová matice, takže při DCT se budou používat pouze základní aritmetické operace. Přibylo několik pomocných proměnných.
CT a CTt jsou kosinové tabulky, navzájem pouze převrácené. To proto, aby se zjednodušily na první pohled zvláštní přemety, které se uplatňují při násobení matic. T je pomocná matice (temp). program DCT256f; { fast image compression using DCT
by Libor Tinka }
uses Crt,Graph; const kvalita = 15; { <1,+-25> } type Matice = array[0..7,0..7] of real; var d,m, xx,yy,x,y,i,j,k: integer; bmp: file of byte; b: byte; A,DCT,QM,CT,CTt,T: Matice; T1: real; procedure Pal(n,R,G,B: byte); begin Port[$3C8]:=n; Port[$3C9]:=R; Port[$3C9]:=G; Port[$3C9]:=B; end; function C(u: real): real; begin if u=0 then C:=1/sqrt(2) else C:=1; end; begin { inicializace grafiky } d:=InstallUserDriver('svga256',nil);m:=0;InitGraph(d,m,'svga256.bgi'); { inicializace palety } for x:=0 to 255 do Pal(x,Round(x*63/255),Round(x*63/255),Round(x*63/255)); { nacteni obrazku } assign(bmp,'exmpl256.bmp');reset(bmp);Seek(bmp,1078); for y:=0 to 199 do for x:=0 to 319 do begin read(bmp,b); PutPixel(x,199-y,b); end; Close(bmp); { *** KOMPRESE *** } { sestaveni kvantizacni matice } for y:=0 to 7 do for x:=0 to 7 do QM[x,y]:=1+((1+x+y)*kvalita); { sestaveni kosinove matice } for j:=0 to 7 do begin CT[0,j]:=1/sqrt(8); CTt[j,0]:=CT[0,j]; end; for i:=1 to 7 do for j:=0 to 7 do begin CT[i,j]:=sqrt(2/8)*cos((2*j+1)*i*pi/(2*8)); CTt[j,i]:=CT[i,j]; end; { prochazeni bloku } for yy:=0 to 24 do for xx:=0 to 39 do begin { nacteni bodu do pixelove matice } for y:=0 to 7 do for x:=0 to 7 do A[x,y]:=GetPixel(xx*8+x,yy*8+y); { snizeni hodnot } for y:=0 to 7 do for x:=0 to 7 do A[x,y]:=A[x,y]-128; { ** DCT ** } { T, A, CTt } for i:=0 to 7 do begin for j:=0 to 7 do begin T[i,j]:=0; for k:=0 to 7 do T[i,j]:=T[i,j]+A[i,k]*CTt[k,j]; end; end; { DCT, CT, T } for i:=0 to 7 do for j:=0 to 7 do begin T1:=0; for k:=0 to 7 do T1:=T1+CT[i,k]*T[k,j];
DCT[i,j]:=T1; end; { kvantizace } for j:=0 to 7 do for i:=0 to 7 do DCT[i,j]:=Round(DCT[i,j]/QM[i,j]); { dekvantizace } for j:=0 to 7 do for i:=0 to 7 do DCT[i,j]:=DCT[i,j]*QM[i,j]; { *** I-DCT *** } { T, DCT, A } for i:=0 to 7 do for j:=0 to 7 do begin T[i,j]:=0; for k:=0 to 7 do T[i,j]:=T[i,j]+DCT[i,k]*CT[k,j]; end; { A, CTt, T } for i:=0 to 7 do for j:=0 to 7 do begin T1:=0; for k:=0 to 7 do T1:=T1+CTt[i,k]*T[k,j]; T1:=T1+128; if T1<0 then A[i,j]:=0 else if T1>255 then A[i,j]:=255 else A[i,j]:=Round(T1); end; { vykresleni } for y:=0 to 7 do for x:=0 to 7 do PutPixel(xx*8+x,yy*8+y,Round(A[x,y])); if KeyPressed then { preruseni } begin CloseGraph;Halt; end; end; readkey; CloseGraph; end.
Kosinovou matici vytvoříte podle tohoto pravidla (kódem je to vyjádřeno pod komenářem: sestavení kosinové matice):
Rutiny DCT a I-DCT pak mohou být ralizovány zcela jinak, a to násobením matic. Pozor! Zde platí trochu jiná pravidla, než v jakém smyslu se o nich hovořilo u kvantizačních matic. Tady se nenásobí odpovídající hodnoty, ale řádky a sloupce. Symbolicky potom můžeme DCT zapsat jako:
V kódu je T a T´ vyjádřeno jako CT a CTt (CTt je převrácená k CT - transposed). Ve stejném duchu vyjádříme i inverzní DCT:
Výsledky q=1
q=5
q = 10
q = 15
q = 25
q = 50 (q >> 25)
Závěrem Nevím, nakolik jsem vám osvětlil to sladké tajemství, které skrývá specifikace JPEG, ale snažil jsem se vyjádřit co nejvíce informací na co nejmenší prostor. Pokud jste všemu z mého výkladu neporozuměli (já sám nejsem na DCT odborník a moje znalosti jsou víceméně povrchní), rozhodně nahlédněte do některé z mnoha dokumentací na webu. Já jsem pochytil některé "fígle" ve skvělé knize "Data Compression Book" v elektronické podobě a úvodu do DCT: "Image Compression and the Discrete Cosine Transform" od Kena Cabeena a Petera Genta. Kód můžete volně šířit a modifikovat, od toho jsem ho napsal, jen prosím uveďte mé jméno jako původního autora.
