Wykład nr 9
Transformata kosinusowa - zastosowanie transformaty kosinusowej w formacie jpeg
przygotowała Klaudia Wrzask

Off[General::"spell"]
Off[General::"spell1"]

Dyskretna transformata Kosinusowa - cztery reprezentacje

Dyskretna transfoemata kosinusowa została po raz pierwszy przedstawiona przez N. Ahmed, T. Natarajan i K. R. Rao w 1974 roku

DCTI
(Subscript[X, k])= 1/2 ((Subscript[x, 0])+(((-))^k)(Subscript[x, N-1]))+(Underoverscript[ ∑, (n=1), (N-2)])(Subscript[x, n]) cos((( π)/(N-1))n k)   k=0,...,N-1

DCTI[x_]:=
Table[
1/2 (x[[1]]+(-1)^(k-1)*x[[Length[x]]])+
Sum[x[[n]]*Cos[ Pi  (n-1) ( k-1)/(Length[x]-1)],{n,2,Length[x]-1}],
{k,1,Length[x]}]

data=Table[N[Sin[((12 π n)/240)]]+N[Sin[((20 π n)/240)]],{n,1,240}];
ListPlot[data,Joined->True]
DCTI[data];

W Mathematice transformata kosinusowa to FourierDCT[dane,typ] z dokładnością do czynnika normalizującego

{ListPlot[DCTI[data],Joined->True,PlotRange->All],ListPlot[FourierDCT[data,1],Joined->True,PlotRange->All]}

DCTII
(Subscript[X, k])= (Underoverscript[ ∑, (n=0), (N-1)])(Subscript[x, n]) cos((( π)/N)(n+(1/2)) k)   k=0,...,N-1

DCTII[x_]:=
Table[
Sum[x[[n]]*Cos[ Pi  (n-1/2) ( k-1)/(Length[x])],{n,1,Length[x]}],
{k,1,Length[x]}]

{ListPlot[DCTII[data],Joined->True,PlotRange->All],ListPlot[FourierDCT[data,2],Joined->True,PlotRange->All]}

DCTIII
(Subscript[X, k])= (1/2)(Subscript[x, 0])+(Underoverscript[ ∑, (n=0), (N-1)])(Subscript[x, n]) cos((( π)/N)(n+(1/2)) k)   k=0,...,N-1

DCTIII[x_]:=
1/2 x[[1]]+
Table[
Sum[x[[n]]*Cos[ Pi  (n-1) ( k-1/2)/(Length[x])],{n,2,Length[x]}],
{k,1,Length[x]}]

{ListPlot[DCTIII[data],Joined->True,PlotRange->All],ListPlot[FourierDCT[data,3],Joined->True,PlotRange->All]}

DCTIV
(Subscript[X, k])= (Underoverscript[ ∑, (n=0), (N-1)])(Subscript[x, n]) cos((( π)/N)(n+(1/2)) (k+(1/2)))   k=0,...,N-1

DCTIV[x_]:=
Table[
Sum[x[[n]]*Cos[ Pi  (n-1/2) ( k-1/2)/(Length[x])],{n,1,Length[x]}],
{k,1,Length[x]}]

{ListPlot[DCTIV[data],Joined->True,PlotRange->All],ListPlot[FourierDCT[data,4],Joined->True,PlotRange->All]}

I jeszcze test szybkości:

Clear[test];
n=2^10;
test[n_]:=Table[Random[],{i,1,n}];

Timing[DCTII[test[2^10]];][[1]]

Timing[FourierDCT[test[2^10]];][[1]]

Transformata kosinusowa ma złożoność obliczeniową O((N^2))

TableForm[Table[{2^n,Timing[DCTII[test[2^n]];][[1]],Timing[DCTII[test[2^n]];][[1]]/(2^(2 n) )},{n,1,10}]]

Szybka transformata kosinusowa zaimplementowana w Mathematice ma złożoność obliczeniową O(N Log[2,N])

TableForm[Table[{2^n,Timing[FourierDCT[test[2^n]];][[1]],Timing[FourierDCT[test[2^n]];][[1]]/(2^n Log[2,2^n])},{n,15,25}]]

Odwrotne transformaty kosinusowe

odwrtona DCTI jest równa DCTI z dokładnością do czynnika normalizującego 2/(N-1)

InvDCTI[x_]:=2/(Length[x]-1) *DCTI[x]

data2={1,2,3,4,5};

DCTI[data2]//MatrixForm

Simplify[InvDCTI[DCTI[data2]]]//MatrixForm

odwrotna DCTII jest równa DCTIII z dokładnością do czynnika normalizującego 2/N

InvDCTII[x_]:=2/(Length[x]) *DCTIII[x]

Simplify[InvDCTII[DCTII[data2]]]//MatrixForm

odwrotna DCTIII jest równa DCTII z dokładnością do czynnika normalizującego 2/N

InvDCTIII[x_]:=2/(Length[x]) *DCTII[x]

Simplify[InvDCTIII[DCTIII[data2]]]//MatrixForm

odwrotna DCTIV jest równa DCTIV z dokładnością do czynnika normalizującego 2/N

InvDCTIV[x_]:=2/(Length[x]) *DCTIV[x]

FullSimplify[InvDCTIV[DCTIV[data2]]]//MatrixForm

Własności transformaty kosinusowej

✔ Transformata kosinusowa jest rzeczywista
✔ Transformata kosinusowa ma własność skupiania energii w jej początkowych współczynnikach, bardzo przydatną z punktu widzenia kompresji. Poniżej dane do testów

data3=Table[8 n,{n,1,16}]

ListPlot[data3]

Na wykresach poniżej widać, że widmo transformaty kosinusowej jest bardziej skoncentrowane w początkowych współczynnikach w porównaniu z transformatą Fouriera.

ListPlot[Abs[Fourier[data3]],Filling->Axis,PlotLabel->"DFT"]
ListPlot[DCTII[data3],Filling->Axis,PlotLabel->"DCT",PlotRange->All]

Wykonujemy “naiwną” kompresję motodą progową, tzn. zerujemu najniższe współczynniki

ListPlot[Abs[Chop[Fourier[data3],50]],Filling->Axis,PlotLabel->"DFT"]
ListPlot[Chop[N[DCTII[data3]],200],Filling->Axis,PlotLabel->"DCT",PlotRange->All]

Poniżej wykonujemy odwrotną transformatę w obu przypadkach, widać że "kształt sygnału" po kompresji kosinusowej jest lepiej zachowany.

ListPlot[data3,Filling->Axis,PlotLabel->"dane wejściowe"]
ListPlot[Re[InverseFourier[Chop[Fourier[data3],50]]],Filling->Axis,PlotLabel->"kompresja  z DFT"]
ListPlot[InvDCTII[Chop[N[DCTII[data3]],200]],Filling->Axis,PlotLabel->"kompresja z DCT",PlotRange->All]

Dwuwymiarowa transformata kosinusowa

✔ Dwuwymiarowa transformata kosinusowa to transformata jednowymiarowa wykonana na wszystkich wierszach, po czym w kolejnym kroku wykonujemy transformatę jednowymiarową z tak otrzymanego sygnały na wszysktich kolumnach.
✔ Kolejność wiersze kolumny jest dowolna

Poniżej procedura na dwuwymiarową transformatę na podstawie DCTII

DCT2D[x_]:=
N[Table[
Sum[
x[[n,m]]*Cos[ Pi  (n-1/2) ( k-1)/(Length[x])]*Cos[ Pi  (m-1/2) ( l-1)/(Length[x])],
{n,1,Length[x]},{m,1,Length[x]}
],{k,1,Length[x]},{l,1,Length[x]}
]]

DCT2D[{{1,2,1},{2,3,4},{4,5,6}}]

I odwrotna

c[n_]:=If[n==1,1/2,1];
InvDCT2D[x_]:=
N[
4/(Length[x])^2* 
Table[
Sum[c[n]c[m]  x[[n,m]]*Cos[ Pi  (n-1) ( k-1/2)/(Length[x])] Cos[ Pi  (m-1) (l-1/2)/(Length[x])],{n,1,Length[x]},{m,1,Length[x]}],
{k,1,Length[x]},{l,1,Length[x]}]]

InvDCT2D[DCT2D[{{1,2,1},{2,3,4},{4,5,6}}]]

Kompresja grafiki z transfromatą kosinusową - metoda “naiwna”

letters=Import["c:/Users/wrzask/Desktop/PO/oldletters.jpg"]
dataletters=ImageData[letters];
ImageDimensions[letters]
Dimensions[dataletters]
{Min[dataletters],Max[dataletters]}

Poniżej wycinam sobie literkę A - w Mathematice pierwszy element tablicy tzn [[1,1]] do lewa górna krawędź grafiki:

datalettersA=dataletters[[Range[21,180],Range[24,151],1]];
Dimensions[datalettersA]
Image[datalettersA]

{Image[FourierDCT[datalettersA]],
Image[RotateRight[Abs[Fourier[datalettersA]],{80,64}]]}

ListPlot3D[FourierDCT[datalettersA],Mesh->False,PlotRange->All]
ListPlot3D[RotateRight[Abs[Fourier[datalettersA]],{80,64}],Mesh->False,PlotRange->All]

Histogram[Flatten[FourierDCT[datalettersA]]]
Histogram[Flatten[Abs[Fourier[datalettersA]]]]

Kompresja metodą naiwną
✔ Wykonujemy dwuwymiarowa transformatę
✔ Zerujemu wspołczynniki transformaty poniżej pewnego progu
✔ Wykomujemy odwrotną transformatę

Kompresje około 
13% dla f=0.01
37% dla f=0.03
55% dla f=0.05
78% dla f=0.1
99% dla f=0.5

f=0.5;
chopA=Chop[FourierDCT[datalettersA],f];
N[Dimensions[Select[Flatten[FourierDCT[datalettersA]],Abs[#]<f&]]/Dimensions[Flatten[FourierDCT[datalettersA]]]]

comprA=FourierDCT[chopA,3];
{Image[comprA],Image[datalettersA]}

Kompresja grafiki jpeg (Joint Photographic Experts Group)

Kodowanie jpeg przebiega następująco:

✔ Krok 1 obraz jest konwertowany z kanałów czerwony-zielony-niebieski (RGB) na jasność (luminancję) i 2 kanały barwy (chrominancje) YCbCr (Patrz wykład nr 1). Wynika to z tego, że dokładniej postrzegamy drobne różnice jasności od drobnych różnic barwy. Użyteczne staje się użycie różnych parametrów kompresji na różnych kanałach. Taka transformacja na kolorach nie jest w pełni odwracalna. Krok ten nie jest obowiązkowy.
✔ Krok 2 wstępnie odrzucana jest część pikseli kanałów barwy - zazwyczaj jest to co druga próbka. (tzw. downsampling)
✔ Krok 3 kanały są dzielone na bloki 8 × 8 pikseli. W przypadku kanałów kolorów, jest to 8 × 8 danych po redukcji rozdzielczości, a więc zwykle 16 × 16 pikseli obrazu początkowego.
✔ Krok 4 na blokach wykonywana jest dyskretna transformata kosinusowa (DCT). Zamiast wartości pikseli mamy teraz średnią wartość wewnątrz bloku oraz częstotliwości zmian wewnątrz bloku. Transformata DCT jest odwracalna.
✔ Krok 5 Kwantyzacja, czyli zastąpienie danych zmiennoprzecinkowych przez liczby całkowite. To właśnie tutaj występują straty danych. Zależnie od parametrów kompresora, odrzuca się mniej lub więcej danych. Zasadniczo większa dokładność jest stosowana do danych dotyczących niskich częstotliwości niż wysokich.
✔ Krok 6 współczynniki DCT są uporządkowywane zygzakowato, aby zera leżały obok siebie. (entropy coding)
✔ Krok 7 współczynniki niezerowe są kompresowane algorytmem Huffmana. Są specjalne kody dla ciągów zer.

Użyta transformata powoduje efekty blokowe w przypadku mocno skompresowanych obrazków.

Poniżej zajmiemy się krokami 3-5:

Krok 3 blokowanie

datalettersA=dataletters[[Range[21,180],Range[24,151],1]];
Dimensions[datalettersA]
Image[datalettersA]

W tym miejscu dzielę obrazek na bloki 8x8, w tym celu korzystam z Partition[]

blockA=Partition[datalettersA,{8,8}];
Dimensions[blockA]

Image[blockA[[10,12]]]

blockA[[10,12]]

W tym miejscu łączę bloki w całość, w tym celu korzystam z ArrayFlatten[]

Dimensions[ArrayFlatten[blockA]]
Image[ArrayFlatten[blockA]]

Krok 4 transformata kosinusowa
Krok 5 kwantyzacja

Poniżej: 
transformata kosinusowa na blokach
kwantyzacja - czyli kompresja polegajaca na zaokrągleniu wartości: f=0.01 do dwóch miejsc po przecinku, f=0.1 do jednego miejsca po przecinku,...
odwrotna transformata kosinusowa

f=1;
chopblockA=Table[Round[FourierDCT[blockA[[i,j]]],f],{i,1,20},{j,1,16}];
Dimensions[chopblockA]
comprblockA=Table[FourierDCT[chopblockA[[i,j]],3],{i,1,20},{j,1,16}];
Dimensions[comprblockA]

{Image[ArrayFlatten[comprblockA]],Image[datalettersA]}

Pytania na kolokwium

1. Cztery typy transformaty kosinusowej oraz odwrotnej transformaty kosinusowej - p.p.
2. Jakie są własności transformaty kosinusowej? - p.p.
3. W jakich krokach przebiega kodowanie jpeg, opisz bardziej szczegółowo jeden z kroków? - p.r.

Zadania na laboratorium

nie ma podstawowych zadań do tego wykładu

Dodatkowe zadania na laboratorium (projekty)

Zadanie D - ok 10p.
✔ Porównaj kompresję grafik metodą naiwną (mogą być literki, ale nie koniecznie) przy użyciu transformaty Fouriera i transformaty kosinusowej przy tych samych stopniach kompresji np. 50%, 70%, 90% 
✔ Wyeksportuj grafikę i sprawdź jak te stopnie kompresji mają się do rzeczywistej wagi plików.

Źródła/literatura

Andrew B. Watson, NASA Ames Research Center,
Image Compression Using the Discrete Cosine Transform
http://vision.arc.nasa.gov/publications/mathjournal94/mathjournal94.pdf