Special article
Computer-based analysis of human skin rashes

Wstęp
W dermatologii ocenia się u pacjenta charakter wykwitów, tzn. patologicznych zmian towarzyszących różnym chorobom, przed, w trakcie i po zakończonym leczeniu. W ocenie tej bierze się pod uwagę m.in. następujące parametry: wielkość (powierzchnię), liczbę i kolor. Do pomiaru tych parametrów można zastosować techniki tradycyjne poprzez liczenie i mierzenie takich wykwitów oddzielnie, wykorzystując na przykład planimetrię manualną. Procedura taka jest pracochłonna, a przy dużej liczbie chorych lub konieczności zmierzenia bardzo wielu wykwitów, na przykład plam bielaczych, błąd takiej metody może być zbyt duży. Potrzebne są zatem metody cyfrowej analizy obrazów skóry. Obecnie wykorzystuje się różne zjawiska fizyczne do generacji obrazów medycznych, a postacie obrazów i sposoby ich przetwarzania są bardzo różnorodne. Wspólnym celem technik ich wstępnego przetwarzania jest poprawa jakości obrazu poprzez metody filtrowania, segmentacji oraz redukcji zakłóceń szumowych w stopniu gwarantującym ich przydatność diagnostyczną. Praca ta jest efektem współpracy między Zakładem Algorytmiki i Programowania na Wydziale Matematyki i Informatyki UAM i Katedrą i Kliniką Dermatologii Akademii Medycznej, zapoczątkowanej na początku 2004 r. Dotychczasowe wspólne badania doprowadziły do opracowania pierwszej, podstawowej wersji oprogramowania komputerowego, które pozwala dokonywać analizy wykwitów skóry o względnie regularnej strukturze (mamy tutaj na uwadze zarówno kształt wykwitu, jak i jego zabarwienie). Zastosowanie specjalistycznego oprogramowania ma na celu z jednej strony istotne skrócenie procesu pomiaru i analizy stopnia zaawansowania choroby i skuteczności stosowanej terapii, a z drugiej strony, co należy wyraźnie podkreślić, zapewnienie obiektywności uzyskanych wyników. Przyjęto trzy podstawowe założenia praktycznej realizacji projektu. 1. Pomiary są wykonywane na podstawie zdjęć cyfrowych fragmentów skóry dotkniętych schorzeniami. 2. Zdjęcia wykonywane są ogólnie dostępnym aparatem cyfrowym o średniej lub wysokiej rozdzielczości, w warunkach odpowiadających tym w gabinecie lekarskim. 3. Projektowane oprogramowanie musi realizować swoje funkcje w akceptowalnym czasie, na przeciętnym komputerze klasy PC.
Trzecie założenie było nadrzędnym wymaganiem. Oznaczało to konieczność zastosowania metod i algorytmów, które działają w akceptowalnym czasie na komputerach o stosunkowo niewielkiej mocy obliczeniowej i przy znacznych ograniczeniach pamięciowych. Wymaganie to skłoniło nas do rozważenia innych metod niż tradycyjna analiza i klasyfikacja obrazów związana z segmentacją i porównywaniem poszczególnych próbek z przygotowaną wcześniej bazą danych próbek wzorcowych. Nasze podejście, oparte na metodach przyrostowych, istotnie różni się od stosowanego [1] w sieci oscylatorów, choć cel pracy jest podobny.
Sformułowanie problemu
Dane jest zdjęcie fragmentu skóry ludzkiej w postaci elektronicznej, czyli zbioru n punktów (zwanych dalej pikselami), przy czym n=w x h, co interpretowane jest jako h wierszy, z których każdy zawiera w pikseli. Każdemu pikselowi przyporządkowany jest trzyelementowy wektor [R, G, B] reprezentujący jego kolor (odpowiednio składową czerwoną, zieloną i niebieską). Jest to typowa reprezentacja obrazów w pamięci komputera [2]. Zakładamy, że analizowane zdjęcia są wykonane w rozdzielczości możliwej do uzyskania w dostępnych obecnie cyfrowych aparatach fotograficznych (co najmniej 3 mln pikseli). Zakładamy ponadto, że zdjęcia są wykonywane w warunkach odpowiadających gabinetowi lekarskiemu, przy różnym oświetleniu, różnym poziomie zbliżenia i ostrości. Zbiór pikseli tworzących zdjęcie, który oznaczać będziemy przez I, zawiera pewną liczbę rozłącznych podzbiorów pikseli S1, S2,..., Sk, które reprezentują obszary wykwitów (będziemy w skrócie nazywać je obszarami chorymi). Wynikiem analizy powinno być określenie liczby k chorych obszarów oraz powierzchni każdego z obszarów S1,... Sk, wyrażonej w dogodnych jednostkach (zwykle w mm2). Problem został podzielony na następujące etapy: 1) wymiarowanie – wyrażenie zdefiniowanych wyżej wielkości w oraz h w milimetrach, 2) identyfikacja obszarów – dla każdego obszaru S1,... Sk, znalezienie co najmniej jednego piksela należącego do danego obszaru, 3) ograniczenie obszaru – startując z punktu należącego do obszaru Si, wyznaczenie całego obszar Si i, korzystając z danych wyjściowych etapu 1, obliczenie jego powierzchni. Dodatkowo przyjęto następujące założenia: 1. Zastosowane metody muszą być wydajne w tym sensie, że ich implementacja będzie generować wyniki w akceptowalnym czasie na komputerze klasy PC. 2. Algorytmy muszą być stabilne i deterministyczne, co rozumiemy jako niewielkie zaburzenie wyników przy niewielkim zaburzeniu danych wejściowych oraz uzyskanie tych samych wyników przy wielokrotnej analizie tych samych danych wejściowych. 3. Program rozwiązujący postawiony problem musi dopuszczać ingerencję użytkownika (lekarza) w przypadku niezadowalających wyników uzyskanych metodą automatyczną.

Narzędzia i metody
Analiza zbioru zdjęć skłoniła autorów do sformułowania następujących wniosków, które miały istotny wpływ na dobór metod algorytmicznych rozwiązania postawionego wyżej problemu: • zarówno obszary zdrowe, jak i chore wykazują daleko posuniętą niejednorodność, co do koloru i struktury; • poziom jasności i nasycenia kolorów poszczególnych zdjęć istotnie się różni; • na zdjęciach, oprócz obszarów skóry ludzkiej, występują także inne obiekty, jak fragmenty ubrania, przyrządy pomiarowe, uchwycone w kadrze części tła; • zarówno w przypadku obszarów zdrowych, jak i chorych mamy często do czynienia z owłosieniem, dodatkowo zwiększającym stopień niejednorodności obszarów.
Wymiarowanie zdjęcia
Problem wymiarowania, czyli przejścia z abstrakcyjnych jednostek – pikseli – na jednostki metryczne, został rozwiązany poprzez umieszczenie na każdym zdjęciu charakterystycznego elementu o znanych z góry wymiarach. W prototypie oprogramowania przyjęto, że jest to kwadrat w kolorze zielonym o ustalonej długości boku. Automatyczne wymiarowanie polega na wyznaczeniu zbioru pikseli tworzących zielony kwadrat, a następnie znalezienie jego wierzchołków (za pomocą prostych równań algebraicznych). Umożliwia to przeliczenie długości boku kwadratu wyrażonej w pikselach na długość tego boku wyrażoną w milimetrach, co w rezultacie daje możliwość wyznaczenia rozmiarów pojedynczego piksela na zdjęciu.

Wyznaczanie granic obszaru
Załóżmy, że dany jest zbiór punktów {s₁,..., s_k}, o których wiemy, że należą do obszaru chorego, oraz zbiór punktów {h1,..., hm}, o których wiemy, że nie należą do obszaru chorego. Do rozwiązania problemu wyznaczenia chorego obszaru Si zawierającego punkty {s₁,..., s_k} użyte zostały algorytmy przyrostowe, tzn. takie metody, które startując z jednego bądź więcej punktów ze zbioru {s₁,..., s_k}, w kolejnych iteracjach powiększają szukany chory obszar, aż do momentu wyznaczenia całego obszaru. Algorytmy przyrostowe udowodniły swoją wysoką wydajność również w innych zagadnieniach związanych z grafiką komputerową (por. [3]). W celu zoptymalizowania wydajności algorytmu zastosowano metodę zmniejszającego się kroku wyznaczania obszaru. Jest to metoda iteracyjna, która w pierwszym etapie przyjmuje względnie duży krok przeszukiwania (rzędu 0,5–1,0 mm) i wyznacza zgrubnie (ryc. 1.) granice chorego obszaru. W kolejnych krokach zajmujemy się wyłącznie coraz dokładniejszym wyznaczeniem brzegu obszaru, zmniejszając w każdej iteracji krok o połowę (ryc. 2.), aż do momentu uzyskania kroku o długości nie przekraczającej rozmiaru pojedynczego piksela na zdjęciu. Zastosowanie takiego podejścia umożliwia pominięcie analizy dużego zbioru pikseli z wnętrza chorego obszaru i wykorzystanie mocy obliczeniowej komputera do dokładnej analizy jego obrzeża. Opisywana metoda iteracyjna ma dodatkowo tę zaletę, że pochłania drobne zakłócenia występujące wewnątrz chorego obszaru, dzięki czemu można w większości przypadków uniknąć konieczności stosowania operacji domykania obszaru (pochłaniania zakłóceń), np. za pomocą operacji morfologicznego domknięcia [4]. Po wyznaczeniu zewnętrznego konturu obszaru, przynależność do niego poszczególnych punktów wewnętrznych można wyznaczyć na podstawie jednego z algorytmów przynależności do wielokąta [5]. W niektórych przypadkach znaczne polepszenie skuteczności algorytmu można uzyskać przez wstępne przetworzenie obrazu, np. za pomocą filtrów medianowych, wygładzających i uśredniających [6].
Klasyfikacja punktów
Z poczynionych wcześniej założeń wynika, że dany jest co najmniej jeden punkt obrazu, o którym wiemy, że należy do obszaru chorego (nazwijmy go s), i co najmniej jeden, o którym wiemy, że nie należy do obszaru chorego (nazwijmy go h). Rozważając w danym kroku algorytmu pewien punkt p, musimy zakwalifikować go albo do klasy punktów chorych, albo do klasy punktów zdrowych. Musimy zatem wykonać porównanie punktu p (i ewentualnie jego otoczenia) z punktami s i h (i ewentualnie z ich otoczeniami) oraz stwierdzić, czy badany punkt p znajduje się bliżej (przy pewnej definicji odległości) punktu s, czy punktu h. Przypomnijmy, że punkty (piksele) reprezentowane są jako trzyelementowe wektory [R, G, B], definiujące kolor danego punktu. Jako miarę odległości w sensie koloru między punktami x i y przyjmuje się zwykle trójwymiarową odległość euklidesową (traktując składowe koloru R, G, B jako współrzędne wektora w przestrzeni trójwymiarowej) lub jej warianty ważone, w których nie wszystkie składowe wpływają w równym stopniu na odległość (por. [7, 8]). Prototyp oprogramowania użyty do wstępnej oceny metod wykorzystuje zarówno klasyczną odległość euklidesową, jak i jej warianty ważone (z których jeden uwzględnia dodatkowo luminancję – jasność badanego punktu). Już pierwsze testy wykazały, że badanie podobieństwa pojedynczych punktów może dać interesujące rezultaty tylko w przypadku bardzo wąskiej grupy obrazów. Metoda ogólna musi uwzględniać nie tylko samą odległość (w sensie koloru) między parami punktów (p, s) i (p, h), lecz również między otoczeniem tych punktów. Jako otoczenie punktu p o promieniu r przyjmujemy zbiór punktów N (p), których odległość (w klasycznym rozumieniu dystansu) od p nie przekracza r. Najprostszą metodą określania podobieństwa punktów p, s wraz z ich otoczeniem jest obliczenie średniego koloru po wszystkich punktach ze zbiorów N (p) i wyznaczenie odległości (w sensie koloru) tak uśrednionych wartości. W praktyce jednak dużo lepsze rezultaty udało się osiągnąć, stosując średnią ważoną koloru otoczenia. Ogólna zasada wyznaczania średniego koloru ważonego polega na pomnożeniu każdego koloru punktu ze zbioru N (p) przez pewną wagę, której wartość jest odwrotnie proporcjonalna do odległości (w zwykłym sensie) punktu z otoczenia od punktu p. W prototypie oprogramowania na otoczenia porównywanych punktów p i s (odp. p i h) nałożona została powierzchnia Gaussa (ryc. 3.), definiująca wagi kolorów poszczególnych punktów z otoczenia (powierzchnie tego typu wykorzystywane są szeroko m.in. w problemach skalowania obrazów [9]).

Wyznaczanie punktów startowych
We wcześniejszych rozważaniach zakładaliśmy, że dany jest co najmniej jeden punkt należący do szukanego chorego obszaru oraz co najmniej jeden punkt nie należący do tego obszaru. Takie dane wejściowe umożliwiają zastosowanie wydajnej metody przyrostowej. Do rozwiązania pozostaje jednak problem wyznaczenia punktów startowych. Podejście klasyczne (oparte na dopasowaniu wzorca) polega na porównywaniu obszarów przetwarzanego obrazu z bazą danych próbek i klasyfikowaniu poszczególnych obszarów na podstawie podobieństwa do konkretnych próbek znanego typu. Takie rozwiązanie może być zastosowane w połączeniu z algorytmem przyrostowym (algorytm ten nie czyni żadnych założeń, co do pochodzenia punktów startowych), jest ono jednak mało wydajne z punktu widzenia czasu obliczeń oraz wymaga wcześniejszego przygotowania bazy danych próbek. W naszych rozważaniach pozostawiamy problem wyznaczania punktów startowych jako otwarty, prowadząc jednocześnie intensywne badania nad zastosowaniem technik znanych z geometrii fraktalnej (dokładniej nad wykorzystaniem wymiaru fraktalnego obszaru – por. np. [2, 10]) w celu przyspieszenia procesu wyznaczania punktów startowych.
Wyniki
Z bazy ponad 7000 fotografii obrazów różnych stanów chorobowych wybrano 320 fotografii wykonanych aparatem cyfrowym. Wśród nich znalazły się przypadki chorych ze znamionami barwnikowymi, naczyniakami płaskimi i jamistymi, z bielactwem, łuszczycą, łysieniem plackowatym, keloidami, trądzikiem pospolitym, owrzodzeniami żylakowatymi podudzi i niektórymi postaciami nowotworów skóry (czerniak złośliwy, nabłoniaki podstawno-komórkowe lub kolczysto-komórkowe). Prototypowe oprogramowanie, wykorzystywane do oceny stosowanych metod i algorytmów, poprawnie wyznacza obszary w przypadku wyraźnej różnicy między strukturą obszaru zdrowego i chorego (ryc. 4.). Uzyskany czas obliczeń dla tego typu obszarów również jest zadowalający (zwykle poniżej 1 s na średniej klasy komputerze PC). Na przedstawionych ilustracjach (ryc. 5., 6.) kolorem czerwonym zostały oznaczone obszary wyznaczone komputerowo poprzez zastosowanie opisanej wyżej procedury iteracyjnej połączonej z porównywaniem otoczenia punktów opartym na średniej koloru, ważonej zgodnie z nałożoną powierzchnią Gaussa. W przypadkach słabo wyróżniających się obszarów (ryc. 5.) lub zdjęć wyjątkowo niskiej jakości, np. zawierających zakłócenia wywołane lampą błyskową aparatu o powierzchni zbliżonej do powierzchni całego szukanego obszaru (ryc. 6.) zastosowane algorytmy wyznaczają część chorego obszaru, jednak ich wydajność (w sensie czasu obliczeń) nadal pozostaje wysoka (zwykle poniżej 1 s), a wyznaczone podobszary są właściwymi podzbiorami faktycznego szukanego obszaru. Istnieje pewna klasa obrazów, dla których metoda przyrostowa w jej obecnej postaci nie daje zadowalających rezultatów w sensie jakości wyznaczonego obszaru. Problem ten dotyczy głównie obrazów, w których szukany obszar jest trudny do wyznaczenia nawet przez specjalistę (ryc. 7.).
W przypadku obrazów, dla których szukany obszar zajmuje znaczną część całego obrazu (ryc. 8.), czas działania algorytmu przyrostowego w naturalny sposób wydłuża się, jednak dzięki zastosowaniu opisanej wcześniej metody iteracyjnej nadal pozostaje akceptowalny (prototyp oprogramowania dla obrazu z ryc. 8. wyznacza obszar w czasie rzędu kilkunastu sekund z bardzo dobrą dokładnością).
Omówienie wyników
Opracowane oprogramowanie może znaleźć wielorakie zastosowanie w dermatologii i innych dziedzinach medycyny, gdzie ocenia się patologię skóry lub innych tkanek, np. w medycynie sądowej, w ocenie powierzchni oparzeń czy dla wspomożenia diagnostyki kardiologicznej [11–13]. W dermatologii, wśród wielu parametrów służących obiektywizacji stanów chorobowych, ich nasilenia, skuteczności podjętej terapii, są m.in. takie wskaźniki, jak PASI dla łuszczycy czy SCORAD dla atopowego zapalenia skóry. Wśród kilku ocenianych parametrów jest także powierzchnia zmienionej chorobowo skóry. Posługiwanie się regułą 9 czy innymi szacunkowymi metodami wiąże się z relatywnie wysokim marginesem błędu. W ocenie zdjęć fotograficznych stanów chorobowych skóry ważna jest nie tylko ostrość, ale wiele innych parametrów, w tym układ barw z wysyceniem. Należy zdawać sobie także sprawę z tego, że metody planimetryczne nie są w stanie dać w pełni obiektywnej oceny trójwymiarowych wykwitów skóry, np. pęcherzy, guzów czy owrzodzeń. Nawet fotografie plam barwnikowych wykonane w różnym oświetleniu (widzialnym, spolaryzowanym, UV i podczerwonym) były różnie oceniane tą samą metodą planimetryczną [14–15].
Podsumowując wyniki eksperymentów przeprowadzonych przy użyciu prototypowego oprogramowania, można stwierdzić, że iteracyjna metoda przyrostowa gwarantuje uzyskanie zadowalającej wydajności procesu wyznaczania obszarów nawet na średniej klasy komputerach PC, co spełnia postawione wcześniej założenie dotyczące efektywności w sensie czasu obliczeń. Jakość wyznaczanych obszarów zależy w głównej mierze od jakości przetwarzanego obrazu. Możemy wyróżnić trzy klasy obrazów z punktu widzenia zachowania się badanego algorytmu: • obrazy, w których metoda automatyczna daje dokładne (z błędem rzędu 1%) wyniki w porównaniu z pomiarem wykonanym ręcznie; do tej klasy należą głównie zdjęcia zawierające schorzenia typu znamiona barwnikowe, znamiona naczyniowe, czerniak złośliwy;
• obrazy, w których metoda automatyczna daje zadowalające wyniki po ustaleniu pewnych dodatkowych parametrów (np. czułości porównania otoczenia punktów, co w przypadku stosowania średniej ważonej powierzchnią Gaussa jest równoważne z dobraniem odpowiednich współczynników tej powierzchni); ta klasa zdjęć obejmuje różne schorzenia, m.in. bielactwo, łuszczycę, wyprysk; • obrazy, dla których wyniki uzyskane metodą przyrostową odbiegają znacznie od wyników uzyskanych za pomocą ręcznego pomiaru; dotyczy to głównie zdjęć o wyjątkowo niskiej jakości lub zawierających znaczną liczbę zakłóceń. Będą to przede wszystkim zdjęcia wykonane ze zmian umiejscowionych na owłosionej skórze głowy. Z przebadanej eksperymentalnie próbki 320 zdjęć większość została zakwalifikowana do pierwszej kategorii, co daje pewne przesłanki do twierdzenia, że iteracyjna metoda przyrostowa, po jej ewentualnie dalszym udoskonaleniu, może stać się przydatnym narzędziem w ocenie rozmiarów objawów schorzeń skóry ludzkiej. Na podkreślenie zasługuje fakt, że prototyp oprogramowania umożliwia ręczną korektę uzyskanych komputerowo wyników przez specjalistę (zarówno w procesie wymiarowania, jak i wyznaczania obszarów), co pozwala na jego zastosowanie, m.in. jako narzędzia elektronicznej planimetrii. Dalszym etapem prac będzie nie tylko pomiar samej powierzchni wykwitów, ale również ich liczby. Ma to np. znaczenie w ocenie skuteczności leczenia trądziku pospolitego czy łuszczycy. Bardziej wyrafinowanymi metodami są już stosowane aplikacje oparte na tzw. sieciach neuronowych (artificial neural networks), służące do cyfrowej analizy obrazów [16]. Nie wydaje się jednak na obecnym etapie wiedzy, aby metody te zastąpiły lekarzy dermatologów.
Piśmiennictwo
1. Zalewska-Janowska A, Strzelecki M, Kwiecień A. Wykorzystanie metod komputerowej analizy obrazu do oceny postępu procesu leczenia owrzodzeń podudzi – badania wstępne. Post Dermatol Alergol 2004; 21: 291-5. 2. Skarbek W. Metody reprezentacji obrazów cyfrowych. Akademicka Oficyna Wydawnicza PLJ. Warszawa 1993. 3. Foley JD, van Dam A, Feiner SK, et al. Wprowadzenie do grafiki komputerowej. Wydawnictwa Naukowo-Techniczne. Warszawa 1995. 4. Pratt W. Digital Imane Processing. 3rd ed. John Wiley & Sons. New York 2001. 5. Preparata FP, Shamos MI. Computational geometry. Springer-Verlag. New York 1985. 6. Gonzalez RC, Woods RE. Digital Image Processing. 2nd ed. Prentice Hall, Upper Sadle River. New Jersey 2002. 7. Alman DH. Industrial color difference evaluation. Color Research and Application 1993; 18: 137-9. 8. Ponton CA. Gamma and its disguises. Journal of the Society of Motion Picture and Television Engineers 1993; 12: 1099-108. 9. Jähne B. Digital Image Processing. Springer-Verlag. New York 2002. 10. Falconer K. Fractal Geometry – Mathematical Foundations and Applications. John Wiley & Sons. Baffis Lane, UK 1990. 11. Schulz-Menger J, Abdel-Aty H, Busjahn A, et al. Left ventricular outflow tract planimetry by cardiovascular magnetic resonance differentiates obstructive from non-obstructive hypertrophic cardiomyopathy. J Cardiovasc Magn Reson 2006; 8: 741-6. 12. Landthaler M, Plewig G, Hofmann C, et al. Autoradiographic and planimetric analyses of the inflammatory infiltrate in psoriatic lesions under PUVA-therapy. Arch Dermatol Res 1980; 269: 31-7. 13. Weissmann A, Bowden J, Frank B, et al. Morphometric studies of the hamster flank organ: an improved model to evaluate pharmacologic effects on sebaceous glands. J Invest Dermatol 1984; 82: 522-5. 14. Cutler NR, George R, Seifert RD, et al. Comparison of quantitative methodologies to define chronic pressure ulcer measurements. Decubitus 1993; 6: 22-30. 15. Korber A, Grabbe S, Dissemond J. 3-dimensional documentation of wound-healing. Zentralbl Chir 2006; 131 (Suppl 1): S174-6. 16. Antkowiak M. Artificial Neural Networks vs Support Vector Machines for Skin Diseases Recognition. 2006, Umea ° University of Computing Science SE-901 87 Umea ° Sweden Umnad 640/06, 2-36.