DCT to skrót od Discrete Cosine Transform – czyli dyskretna transformacja kosinusowa. Jest podstawą wielu koderów obrazu i wideo,  przede wszystkim ciągle najszerzej stosowanego kodera do zdjęć JPEG, a ponadto koderów wideo MPEG-1, MPEG-2, VC-1,  H.263. Popularny H.264 używa wariantu DCT pod nazwą Integer Transform, czasem zwaną ICT (Integer Cosine Transform). Transformacja ta w skrócie, (podobnie jak inne transformacje, na przykład oparte o falki) pozwala na pewnego rodzaju rozkład obrazu na składowe o „różnym stopniu szczegółowości”. Masło maślane, więc w dalszej części postaram się wyjaśnić o co chodzi z DCT i przy okazji jak działają współczesne kodery obrazu działające wg schematu transformacja ↣ kwantyzacja ↣ kodowanie entropijne.

Reprezentacja częstotliwościowa

A gdyby tak dało się rozłożyć gotowy obraz – taki malowany, na płótnie, przez artystę – na poszczególne etapy jego tworzenia. Najpierw pewnie byłoby tło, długimi pociągnięciami szerokiego pędzla, niezbyt dokładnie, tak, żeby złapać kolor – nieba, łąki, morza, lasu. Potem byłyby duże obiekty – góry, lasy, jezioro, brzeg morza, chmury. Potem przyszłaby kolej na mniejsze przedmity, a także na doszczegółowienie tych obiektów, które już się znajdują na obrazie. Aż dochodzimy do najdrobniejszych szczegółów, ostatnich pociągnięć pędzla – tak jak tutaj – u Kevina Hilla na jego kanale youtube.

A potem takie nagranie malarza malującego obraz gdyby tak zagrać od końca. W pewnym sensie byłby to rozkład obrazu na elementy coraz mniej szczegółowe.

W pewnym sensie DCT realizuje właśnie taką funkcję. Przy pomocy tej transformacji możemy rozłożyć dowolny obraz na sumę wzorów graficznych utworzonych ze złożenia funkcji kosinus o różnych częstotliwościach. Ponieważ kodery obrazu i wideo operują na  blokach o boku 4-8-16 pikseli takim obrazem rozkładanym na składowe kosinusowe jest pojedynczy blok.

Poniżej zamieszczam znany z Wikipedii obrazek, zawierający wszystkie funkcje bazowe transformacji DCT dla bloku 8×8 pikseli. Obraz – blok 8×8 – zrekonstruowany z jego reprezentacji częstotliwościowej DCT będzie ważoną sumą wszystkich 64 obrazków.

W lewym górnym rogu funkcja bazowa to po prostu wypełniony kwadrat, odpowiada on składowej stałej bloku. Potem, wzdłuż osi, poziomo i pionowo, mamy wzory powstałe z funkcji kosinus, o rosnącej częstotliwości w miarę oddalania się od skłądowej stałej. Po przekątnej mamy iloczyny kosinusów z obu kierunków, które tworzą charakterystyczne kratki. W prawym dolnym rogu znajduje się po prostu kratka o oczku wielkości jednego piksela – odpowiada ona maksymalnej częstotliwości funkcji kosinus wzdłuż obu osi. Jak widać, najwięcej „szczegółów”, czyli największą zmienność prezentuje funkcja bazowa z prawego dolnego rogu, zaś najmniej „detali” znajdziemy w górnym lewym rogu, czyli w sąsiedztwie składowej stałej (częstotliwość równa zero). Kolor biały oznacza tu wartość maksymalną (np 255 przy 8 bitach precyzji), a czerń minimalną (np 0). Ale to są jedynie funkcje bazowe. Transformacja kosinusowa ma za zadanie ustalić „udział” każdej (każdej z osobna) z tych funkcji bazowych w transformowanym bloku tak, aby nakładając na siebie wszystkie te pofalowane kwadraciki pomnożone przez odpowiadające im udziały otrzymać pierwotny obraz.

Wzór

Czy musiało do tego dojść? Widocznie musiało, skoro doszło. Czyli wklejam wzór na wyliczenie DCT 8×8, wprost z jednego z dokumentów JPEG:

Gdzie Cu, Cv = 1/sqrt(2) dla u,v = 0
oraz Cu, Cv = 1 dla pozostałych

syx to wartość (jasność) piksela w punkcie (x,y), Svu to wartość udziału odpowiedniej funkcji bazowej, o odpowiednich częstotliwościach w obu kierunkach w obrazie pierwotnym (małe s). Tak więc duże S jest czymś w rodzaju dwuwymiarowego widma małego s, czyli pierwotnego obrazu. Widzimy tu też ten iloczyn dwóch kosinusów, który najlepiej widać na przekątnej obrazka z funkcjami bazowymi. Widzimy też, dlaczego wzdłuż krawędzi pionowej i poziomej, biegnących od lewego górnego rogu mamy czyste kosinusy, ponieważ albo u, albo v są dla nich równe 0, a cos(0) = 1, więc jeden z kosinusów znika.

Co z tego wyjdzie?

Pani Lena Sjööblom, Miss Listopada 1972 magazynu Playboy, prawdopodobnie nie przypuszczała, że stanie się osobą tak sławną i to niemal wyłącznie w środowiskach naukowych, a fragment zdjęcia z poświęconej jej rozkładówki magazynu zagości na tysiącach publikacji naukowych i stron internetowych traktujących o cyfrowym przetwarzaniu obrazów. Tradycja zobowiązuje:

Osoby zafascynowane Panią Leną zapraszam na poświęconą jej stronę. Można się tam nawet doklikać do oryginalnego zdjęcia, które nie jest tak ubodze skadrowane i tym samym tak brutalnie pozbawione szerszego kontekstu jak zdjęcie powyżej.

Pozwolę sobie na nieco większy blok, 16×16 pikseli. Wytnę go z kawałka kawałka Pani Leny, a dokładnie z nosa.

Oczywiście DCT policzę dla składowej luminancji, lub jak kto woli – w skali szarości.

I dla tego kawałka wyliczę współczynniki DCT. Przy okazji wzór na DCT  NxN :

Gdzie Cu, Cv = 1/sqrt(2) dla u,v = 0
oraz Cu, Cv = 1 dla pozostałych

Dla powyższego bloku 16×16 pikseli otrzymałem macierz 256 współczynników DCT:

 1984.94    43.52   -29.61     8.26    77.53    10.88    18.19    12.02     3.06     3.61    -3.94    -2.52    -1.47    -3.05    -0.79     0.41 
   50.46  -106.30   -12.71    29.95   -19.49     8.93    -0.20     3.17    -0.75     7.53     2.22     2.47     2.60    -0.63    -1.56     1.42 
   55.00   -22.37    53.24    37.75   -37.65     1.41    -8.37   -16.23    -9.11    -1.46    -1.49     1.77     4.71    -1.02    -1.95    -3.64 
 -116.22   138.91   -34.81   -25.08    36.91   -18.72     2.74     1.43     3.51    -5.74    -2.03     1.82     1.65    -1.56     0.32     2.88 
   42.46   -23.04     9.34     9.02   -22.54    -1.57     5.40     2.27     1.38     7.05     0.46    -1.24    -0.48    -0.19     1.46     4.53 
   35.91   -47.99    16.84    16.11    -9.18     6.34     1.40    -5.89    -1.63     4.44    -0.58    -3.13     0.44    -5.60    -1.77    -3.11 
  -36.92    54.94   -14.27   -20.48    13.39    -4.90     0.09     4.85     1.92    -9.56     1.12    -1.78     2.97    -6.16    -1.24     1.72 
    6.11   -10.87    -2.60     7.73    -2.28    -7.03     5.81     5.12    -3.93    -1.76     1.85    -2.08    -5.26    -0.08     0.29     1.51 
   14.31   -22.79     6.57    -6.35    -7.82     3.29     8.79    -3.17    -5.31     5.30    -3.16    -0.10     3.17    -0.28     3.74    -0.34 
   -6.34     5.74    -0.28    -1.62    -0.21     2.16    -8.08     1.34     2.74     0.04    -1.35     3.40    -1.59     2.82    -0.07    -2.53 
   -6.87    14.18    -2.71    -4.60     2.59    -2.61    -0.33     2.60    -0.31    -2.73    -0.29    -2.46    -1.50     0.06     0.45    -1.22 
    4.91    -8.50     2.50     1.54    -2.06    -3.36     4.42    -0.51    -0.08     1.78    -0.48    -0.55     2.61     0.10     0.65     2.23 
    1.42    -4.91     1.85     1.12    -4.73     3.71    -0.38    -2.31    -0.10     2.71    -0.79     1.99    -2.83    -2.25    -0.25     1.12 
   -9.80    11.18     1.33    -5.37    -0.63    -4.15    -0.80    -0.10     3.08    -2.42     0.54    -3.10    -0.73     0.01    -1.39     1.37 
   -1.73     2.93    -1.91     0.86     4.46    -4.35    -0.23     1.77    -1.36     0.01     0.59     0.09     4.07    -3.56     1.71    -0.21 
    5.23    -9.57     4.38    -1.54    -2.54     5.39    -1.24    -0.37     3.31     3.25     1.57    -2.11     2.97     1.12     2.80    -1.58

 

Pewne charakterystyczne własności DCT nie są może widoczne na pierwszy rzut oka, może oprócz tego, że współczynnik w lewym górnym rogu jest liczbą znacznie większą od pozostałych. To składowa stała, czyli miara „jasności” bloku.
Aby zwizualizować jakoś wartości DCT dla omawianego bloku postanowiłem zamienić macierz DCT w obraz, w tym celu najpierw wyciągnąłem wartośc bezwzględną z każdego współczynnika. Ponieważ różnice między współczynnikami bywają znaczne, wyciągnąłem logarytm z każdego współczynnika, żeby przedstawić macierz w skali logarytmicznej, następnie „ręcznie” zmniejszyłem współczynnik odpowiadający składowej stałej, a na końcu przeskalowałem powstałą macierz do przedziału 0-255, tworząc z macierzy bitmapę-obraz. A tak to wygląda:

Im jaśniejszy piksel, tym większa wartość współczynnika (pamiętamy o skali logarytmicznej, różnice są w rzeczywistości dużo większe). Tutaj od razu widać, że cała „energia”, czyli wielkości modułów współczynników koncentrują się w lewym górnym rogu. Innymi słowy, większość informacji obrazowej zawarta jest w składowych o niższych częstotliwościach, w stronę prawego dolnego rogu mamy coraz większe częstotliwości, czyli coraz drobniejsze szczegóły, których zawartość jest jednak bardzo niewielka. Jest to jednocześnie właściwość samej transformacji DCT (zwana energy compaction), a także właściwość typowych rzeczywistych obrazów, że informacja obrazowa w zasadzie skupiona jest w czestotliwościach niskich i średnich. I tutaj otwiera się pole dla właściwej kompresji. A jeśli dodamy do tego fakt, że oko ludzkie jest dość słabo wyczulone na drobne szczegóły, praktycznie nie zauważa niewielkich zmian kontrastu i przy okazji ma tendencje do „wygładzania” tego co widzi, to można przy okazji do wspomnianej kompresji wprowadzić straty. Nasze oko prawdopodobnie w ogóle nie zarejestruje różnic.

Poniżej znajdziecie wynik transformacji całego obrazka blok po bloku – DCT 16x16px:

Transformacja odwrotna

A co najważniejsze, z macierzy współczynników DCT możemy odtworzyć obraz. Służy do tego odwrotna transformacja kosinusowa, oznaczana IDCT (I jak inverse). Obrazek, który otrzymujemy po zastosowaniu odwrotnego DCT na macierzy współczynników niczym nie będzie się różnił od wyjściowego obrazka.

Jeśli chodzi o wzór, to tutaj niespodzianka:

Gdzie Cu, Cv = 1/sqrt(2) dla u,v = 0
oraz Cu, Cv = 1 dla pozostałych

Niewiele różni się od wzoru na DCT, poza tym, że do sumy trafiły indeksy z domeny współczynników DCT, przez co stałych Cu i Cv nie dało się wywlec przed znaki sum.

Filtrowanie w dziedzinie DCT

Poniżej postaram się zademonstrować to, co napisałem powyżej, czyli właściwość transformaty DCT polegającą na koncentracji informacji w składnikach o niskich częstotliwościach. Z każdego bloku DCT będę kolejno wycinał składniki od coraz wyższych częstotliwości do coraz niższych. Tutaj link do obrazka. Widać, że pomimo całkowitego wycięcia ponad połowy współczynników DCT i zastąpienia ich zerami, obraz po transformacji odwrotnej na pierwszy rzut oka niewiele się różni od oryginalnego obrazu. Po wyrzuceniu około 3/4 współczynników pogorszenie obrazu staje się nieco wyraźniejsze, ale nadal obraz jest bardzo czytelny.

Wygląda na to, że wystarczy zapisać jedynie 1/4 całej informacji obrazowej w postaci współczynników DCT, a widz nie zorientuje się, że obrabowaliśmy go z dużej części informacji zawartej w obrazie.
Na tym mechanizmie opiera się właśnie kompresja obrazu, z tą różnicą, że zamiast wyrzucania współczynników stosuje się odpowiednią kwantyzację, a następnie jakąś metodę kodowania bezstratnego (entropijnego). Dzięki niedoskonałościom naszej percepcji udaje się osiągnąć znaczny stopień kompresji bez dostrzegalnego pogorszenia jakości, dodatkowo w koderach wideo stosuje się kompensację ruchu, aby dodatkowo zwiększyć stopień kompresji wykorzystując naturalną nadmiarowość, jaką niesie ze sobą podobieństwo między dwiema sąsiednimi klatkami wideo.

Kody źródłowe

Liczenie DCT 16×16 wprost ze wzoru podanego wcześnie. Można łatwo wykorzystać poniższy kod do liczenia DCT NxN o dowolnym N, wystarczy zmienić wartość zmiennej BLOCK_SIZE

#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <cmath>
#include <cstdio>

static const char* LENA = "lena_std.tif";
static const char* DCT_BLK = "dct_blk.png";

// Size of the DCT block
static const unsigned BLOCK_SIZE = 16;
// Coordinates of the block to be DCT'ed
static const unsigned BLOCK_Y = 313;
static const unsigned BLOCK_X = 284;

static void printDct(cv::Mat& dctMat);
static void dctToImg(cv::Mat& dctMat, const char* filename);

int main()
{
    cv::Mat img0 = cv::imread(LENA);
    cv::Mat lenaNose = img0(cv::Range(BLOCK_Y, BLOCK_Y + BLOCK_SIZE), cv::Range(BLOCK_X, BLOCK_X + BLOCK_SIZE));

    cv::Mat lenaNoseGrey(cv::Size(BLOCK_SIZE, BLOCK_SIZE), CV_8UC1);
    cv::cvtColor(lenaNose, lenaNoseGrey, CV_BGR2GRAY);

    // 32-bit float for DCT
    cv::Mat lenaNoseDct(cv::Size(BLOCK_SIZE, BLOCK_SIZE), CV_64FC1);
    // Calculate DCT
    for (unsigned v = 0; v < BLOCK_SIZE; ++v) {
        for (unsigned u = 0; u < BLOCK_SIZE; ++u) {
            const double cu = (u == 0) ? 1.0 / sqrt(2) : 1.0;
            const double cv = (v == 0) ? 1.0 / sqrt(2) : 1.0;
            double dctCoeff = 0;

            for (unsigned y = 0; y < BLOCK_SIZE; ++y) {
                for (unsigned x = 0; x < BLOCK_SIZE; ++x) {
                    double uCosFactor = cos((double)(2 * x + 1) * M_PI * (double)u / (2 * (double) BLOCK_SIZE));
                    double vCosFactor = cos((double)(2 * y + 1) * M_PI * (double)v / (2 * (double) BLOCK_SIZE));
                    double pixel = (double)(lenaNoseGrey.at<unsigned char>(cv::Point(x,y)));
                    dctCoeff +=  pixel * uCosFactor * vCosFactor;
                }
            }
            dctCoeff *= (2 / (double) BLOCK_SIZE) * cu * cv;
            lenaNoseDct.at<double>(cv::Point(u,v)) = dctCoeff;
        }
    }

    // Prints all the coefficients to STDOUT
    printDct(lenaNoseDct);

    // Adjust coefficient values and save as image
    dctToImg(lenaNoseDct, DCT_BLK);

    return 0;
}

static void printDct(cv::Mat& dctMat)
{
    for (int yy = 0; yy < BLOCK_SIZE; ++yy) {
        printf("\n");
        for (int xx = 0; xx < BLOCK_SIZE; ++xx) {
            printf("%8.2f ", dctMat.at<double>(cv::Point(xx,yy)));
        }
    }
}

static void dctToImg(cv::Mat& dctMat, const char* filename)
{
    cv::Mat theDct = dctMat.clone();
    // get ABS coefficients
    theDct = cv::abs(theDct);

    cv::log(theDct, theDct);

    // preserve DC component
    double dc = theDct.at<double>(cv::Point(0,0));
    theDct.at<double>(cv::Point(0,0)) = 0;

    // extract max value
    double min, max;
    cv::minMaxLoc(theDct, &min, &max, NULL, NULL);
    // A little fiddling with DC if too high
    if (dc > 1.5 * max) {
        dc = 1.5 * max;
        max = dc;
    } else {
        max = dc > max ? dc : max;
    }
    // Restore DC component
    theDct.at<double>(cv::Point(0,0)) = dc;

    // rescale coefficients to make 8 bit greyscale image
    theDct.convertTo(theDct, -1, 255.0 / max);
    // Write the image
    cv::imwrite(filename, theDct);
}

Przy okazji kawałek liczący odwrotną transformatę:

void idct(cv::Mat& dct, cv::Mat& mat, unsigned blkSize)
{
    // Calculate DCT
    for (unsigned x = 0; x < blkSize; ++x) {
        for (unsigned y = 0; y < blkSize; ++y) {

            double pixel = 0;

            for (unsigned u = 0; u < blkSize; ++u) {
                for (unsigned v = 0; v < blkSize; ++v) {
                    const double cu = (u == 0) ? 1.0 / sqrt(2) : 1.0;
                    const double cv = (v == 0) ? 1.0 / sqrt(2) : 1.0;
                    double uCosFactor = cos((double)(2 * x + 1) * M_PI * (double)u / (2 * (double) blkSize));
                    double vCosFactor = cos((double)(2 * y + 1) * M_PI * (double)v / (2 * (double) blkSize));
                    double coeff = dct.at<double>(cv::Point(u,v));
                    pixel +=  coeff * uCosFactor * vCosFactor * cu * cv;
                }
            }
            pixel *= (2 / (double) blkSize) ;
            mat.at<double>(cv::Point(x,y)) = pixel;
        }
    }
}