PCA, Eigenfaces i rozpoznawanie twarzy

Czym jest PCA?

PCA czyli Principal Component Analysis czyli analiza głównych składowych jest jedną z metod statystycznych służących do analizy zbioru danych. Jeśli kogoś interesuje dokładniejsze omówienie tematu - polecam poszperać choćby po sieci (Wikipedia dla leniwych). Na potrzeby zaś tego wpisu i naszych rozważań uświadomimy sobie jedynie, że PCA umożliwia określenie w jakim stopniu dany klasyfikator jest zgodny z badanym obiektem. Innymi słowy -  jak daleko wartości elementu wzorowego do wartości badanego obiektu.

Czym jest rozpoznawanie twarzy?

Niby oczywistość, ale lepiej ją uściślić teraz, niż potem gubić się w zeznaniach. Z grubsza chodzi o to aby mając zbiór zdjęć z zidentyfikowanymi osobami móc określić czy nowe zdjęcie należy do jakiejś osoby z naszego zbioru (a jeśli tak, to której), czy też jest ono fotografią kogoś zupełnie nowego.

Etapy działania

1. Przygotować zbiór uczący składający się z obrazów o tych samych wymiarach/rozdzielczości (szerokość i wysokość), w których oczy i usta osób znajdują się na podobnej wysokości.
2.  Poddać obrazy operacji wyrównywania histogramów (histogram equalization) w celu redukcji wpływu różnicy oświetlenia.
3. Obliczyć wartości pikseli dla średniego obrazu i odjąć ten obraz od każdego obrazu ze zbioru uczącego.
4. Utworzyć macierz zawierającą wszystkie obrazy (każdy obraz to jeden wektor).
5. Obliczyć wektory własne powyższej macierzy.
6. Obliczyć PCA na podstawie wektorów własnych (jeśli decydujemy się na pominięcie niektórych wektorów to należy odrzucić te z najmniejszymi wartościami własnymi).
7. Wykonać projekcję z wykorzystaniem PCA dla każdego obrazu ze zbioru uczącego.
8. Projekcje nowych zdjęć porównywać z projekcjami obiektów już rozpoznanych i wyszukiwać wyników najbliższych.

Accord.NET

Czym jest Accord.NET? Z grubsza jest to framework oferujący bogactwo funkcji matematycznych i statystycznych, w tym m.in. rozkład na wektory własne, analizę PCA, LDA oraz wykrywanie twarzy (które działa trochę topornie, ale działa). Dlaczego wybrałem ten framework do realizacji tego zadania? Cóż popularną alternatywą byłoby użycie biblioteki OpenCV, która dość szeroko jest do tego celu (rozpoznawania twarzy) rekomendowana. Niestety cud ten technologii nie doczekał się sensownego wsparcia pod C# (wrappery działają różnie i generalnie w mojej opinii nie nadają się do użycia w produktach komercyjnych ze względu na moc różnych drobnych i nie tylko problemów z ich działaniem). Natomiast z racji tego, że Accord.NET jest bezproblemowy w użyciu (i ma wiele przykładów i opisów działania różnych zaimplementowanych funkcji) to praca z nim okazała się bardzo przyjemna.

Eigenfaces = Eigenvectors = Wektory własne

Szperając w poszukiwaniu sposobu użycia PCA w rozpoznawaniu twarzy dość szybko napotkamy na te określenia. Określenia, które z początku wydają się nam sprzeczne z rzeczywistością - bo jak to twarz może być wektorem? Jednak warto pamiętać, że obraz to nic innego jak macierz, której poszczególne wartości odpowiadają poszczególnym pikselom. Każdą zaś macierz możemy zapisać jako wektor jeśli będziemy kolejno "sklejać" jej wiersze.
W poniższych przykładach wykorzystywane będą obrazy zapisane w formacie .pgm (o Netpbm i jego wsparciu dla C# więcej tutaj), który to format umożliwia nam bezproblemowy (i szybki) dostęp do numerycznych wartości każdego piksela.

Tak więc zaczynajmy!
Na początku musimy załadować nasz zbiór uczący, tak aby PCA mogło w ogóle na czymś działać. To w jakiej formie te obrazy będziemy przechowywać na dobrą sprawę nie ma znaczenia - może to być tablica, lista, mogą one być na bieżąco pobierane z bazy. W tym przypadku założyłem, że po prostu wykorzystujemy zapisane na dysku pliki .pgm.

            PgmImage[] rawImage;
            OpenFileDialog ofd = new OpenFileDialog();
            ofd.Multiselect = true;
            ofd.DefaultExt = "pgm";

            if (ofd.ShowDialog() == DialogResult.OK)
            {
                rawImage = new PgmImage[ofd.FileNames.Length];
                for (int i = 0; i < ofd.FileNames.Length; i++)
                {
                    PgmImage aux = (PgmImage)PnmImage.Create(ofd.FileNames[i]);
                    rawImage[i] = HistogramOperation.Equalize(aux);
                }
            }

Tak więc cała nasza baza ląduje w tablicy obrazów typu PgmImage o nazwie rawImage. Warto zauważyć, że od razu poddajemy obrazy wyrównywaniu histogramów, tak więc realizujemy dwa punkty z naszej listy w jednym fragmencie kodu.

Baza zdjęć uczących nasz algorytm - po 2 zdjęcia dla 10 osób

Kolejny krok to utworzenie "obrazu uśrednionego", czyli wyliczenie średniej arytmetycznej dla wszystkich pikseli znajdujących się na danej pozycji. Bardziej opisowo - na nowym "uśrednionym" obrazie wartość piksela znajdującego się w punkcie 0,0 (piksel w górnym lewym rogu) będzie równa średniej arytmetycznej wartości pikseli z punktów 0,0 na wszystkich obrazach zbioru.

            PgmImage avgImage;            
            int mult = (int)(rawImage[0].Size.Width * rawImage[0].Size.Height);
            avgImage = new PgmImage(rawImage[0].Size.Width, rawImage[0].Size.Height, rawImage[0].MaxColor);
            UInt64[] arr = new UInt64[mult];
            UInt16[] avgArr = new UInt16[mult];
            for (int i = 0; i < rawImage.Length; i++)
            {
                for (int j = 0; j < mult; j++)
                {
                    arr[j] += rawImage[i].GetRawData()[j];
                }
            }
            for (int j = 0; j < mult; j++)
            {
                avgArr[j] = (UInt16)(arr[j] / (UInt64)rawImage.Length);
            }
            avgImage.SetRawData(avgArr);

Tak wygląda "przeciętna osoba" z bazy uczącej

Kolejny krok to odjęcie wartości średniej od każdego z obrazów oraz utworzenie macierzy grupującej wszystkie te obrazy (wektory).

            Double[,] T = new Double[rawImage.Length, avgImage.GetRawData().Length];

            for (int i = 0; i < rawImage.Length; i++)
            {
                for (int j = 0; j < avgImage.GetRawData().Length; j++)
                {
                    T[i, j] = rawImage[i].GetRawData()[j] - avgImage.GetRawData()[j];
                }
            }

Tutaj warto przystanąć na chwilkę i pochylić się nad matematyczną stroną tego zadania, albowiem dochodzi tutaj do pewnych przemian, które mają za zadanie zoptymalizowane całego procesu (a konkretniej procesu rozkładu na wartości własne) poprzez ograniczenie danych przesyłanych do algorytmu.
Generalnie potrzebujemy obliczyć wektory własne dla macierzy T. Wszystko pięknie wygląda póki myślimy o małej bazie małych obrazów, jednak jeżeli przyjmiemy, że operujemy na obrazach 92x112 czyli każdy z nich składa się z 10304 punktów i mamy 100 takich obrazów, to nagle dostajemy macierz 100x10304, a każdy nowy obraz dorzuca dodatkowe 10304 wartości, które musimy uwzględnić. Do tego z tej macierzy tworzymy macierz kowariancji C, tak więc wynikiem naszych działań jest potwór 10304x10304, który to trafia do dekompozycji.. a ta będzie trwać.
Tak więc: C*v_i = T*T^T*v_i = λ_i*v_i
Jednak my zdekomponujemy: T^T*T*u_i = λ_i*u_i
co dzięki: T*T^T*u_i = λ_i*T*u_i (co oznacza, że  jeśli u_i jest wektorem własnym T^T*T to v_i = T*u_i jest wektorem własnym C) daje nam macierz o wymiarach 100x100, która porównując z 10304x10304 jest ogromną oszczędnością.

            Double[,] TT = T.Transpose();
            Double[,] C = T.Multiply(TT);
            EigenvalueDecomposition dec = new Accord.Math.Decompositions.EigenvalueDecomposition(C, true);
            Double[,] V = dec.Eigenvectors;
            Double[,] E = TT.Multiply(V);

W tym miejscu możemy przeanalizować również wartości własne poszczególnych wektorów i odrzucić te, które są najmniejsze. Sprawi to, że PCA będzie działało wydajniej nie tracąc jednocześnie ze swojej efektywności. Ważne jest jednak by nie odrzucić zbyt wielu wektorów, gdyż każde odrzucenie wpływa na jakość rozpoznawania.

Z racji tego, że ograniczyliśmy macierz z 10304x10304 do 100x100 koniecznym jest przywrócenie jej "pierwotnych wymiarów" co realizujemy przez mnożenie TT * V. Wartości w macierzy E to w zasadzie wektory własne czyli nasze Eigenfaces, z tą różnicą, że ich wartości nie pokrywają się z przedziałem wartości dostępnych dla pikseli (0-255), tak więc wymagają normalizacji, którą to wykonujemy.

            for (int i = 0; i < E.GetLength(1); i++)
            {
                PgmImage img = (PgmImage)avgImage.Clone();
                Double[] data = new Double[img.GetRawData().Length];
                for (int j = 0; j < avgImage.GetRawData().Length; j++)
                {
                    data[j] = E[j, i];
                }

                data = data.Normalize(0, 255);

                UInt16[] dataRaw = new UInt16[img.GetRawData().Length];
                for (int j = 0; j < avgImage.GetRawData().Length; j++)
                {
                    dataRaw[j] = (UInt16)Math.Round(data[j]);
                }

                img.SetRawData(dataRaw);
                img.Save(outDirectory + "\\eigen-"+i+".pgm", true);
            }

Eigenfaces możemy przechować w naszym programie, ale możemy je również zapisać na dysku. Poniżej klika takich "twarzy własnych" (są rzeczy, których na język polski tłumaczyć nie warto...):

Eigenfaces wygenerowane na podstawie bazy testowej

Jak widać nie odsiano tutaj "gorszych" wektorów własnych, czego efektem jest telewizyjny śnieg na 8 wektorze własnym. Wspomniany wektor oczywiście nie będzie wiele wnosił do jakości rozpoznawania (a wręcz może szkodzić), ale dla celów tej prezentacji warto było go tutaj zostawić.

Mając już wektory własne czyli de facto przetworzoną bazę uczącą czas przejść do wykorzystania PCA.

            Double[,] eigenvectors;            
            PrincipalComponentAnalysis pca;
            pca = new PrincipalComponentAnalysis(eigenvectors);
            pca.Compute();

gdzie w eigenvectors[i, j] i odpowiada za liczbę wektorów,  a j za kolejne wartości punktów w obrazie. Następnie dla każdego obrazu ze zbioru uczącego wykonujemy:

Double[,] projection = pca.Transform(image);

Operacja ta utworzy nam projekcję obrazu opartą o bazę uczącą. Innymi słowy zwróci nam macierz (w naszym przypadku de facto wektor), w której kolejne wiersze zawierać będą informację w jakim stopniu obraz poddany projekcji jest zgodny z daną wartością własną (dla 10 wartości własnych otrzymamy 10 elementów w wektorze projekcji). Zapisując wynikową macierz dla każdego z obrazów jesteśmy następnie w stanie porównać nowe obrazy z tymi już istniejącymi. Porównanie oczywiście można przeprowadzić na pierwszej wartości własnej (jeśli były one uporządkowane to istnieje szansa, że wyniki będą "znośne") lub na całym ich zestawie, przy czym warto pamiętać. że więcej znaczy dokładniej, ale nie szybciej.

A tu taki przykład porównywania projekcji oparty o dość prostą metodę - metrykę euklidesową

            Int64 best = 0;
            Double distance = Double.MaxValue;
            for (int i = 0; i < recognisedPCA.GetLength(0); i++)
            {
                double sum = 0;

                for (int j = 0; j < recognisedPCA.GetLength(1); j++)
                {
                    sum += Math.Pow(projection[0, j] - recognisedPCA[i, j], 2);
                }

                if (sum < distance)
                {
                    distance = sum;
                    best = i;
                }
            }

recognisedPCA - zbiór projekcji rozpoznanych już twarzy
projection - projekcja obecnie wyszukiwanego obrazu
distance - odległość od obecnie najlepszego wzoru
best - id najlepszego wzoru

Bibliografia

• http://www.shervinemami.info/faceRecognition.html
• http://crsouza.blogspot.com/2009/09/principal-component-analysis-in-c.html
• http://en.wikipedia.org/wiki/Eigenface 
• http://www.bytefish.de/blog/eigenfaces
• http://www.bytefish.de/wiki/pca_lda_with_gnu_octave
• http://www.cognotics.com/opencv/servo_2007_series/part_4/index.html
• http://www.cognotics.com/opencv/servo_2007_series/part_5/index.html

Keywords: Eigenvectors, PCA, Principal Component Analysis, Eigenfaces, Eygenwectors, Ejgenwektors, Eygen, Ejgen, Face recognition, Rozpoznawanie twarzy