Identyfikacja kinematycznych modeli ruchu stawu kolanowego na podstawie

P OLITECHNIKA P ^{OZNA ´ NSKA}

R OZPRAWA D OKTORSKA

Identyfikacja kinematycznych modeli ruchu stawu kolanowego na podstawie

sekwencji obrazów RTG

Autor:

mgr in ˙z. Marta D

Promotor:

prof. dr hab. in ˙z. Krzysztof K

Instytut Automatyki i Robotyki Wydział Informatyki

17 czerwca 2019

iii

O´swiadczenie autora

Ja,

o´swiadczam, ˙ze niniejsza rozprawa doktorska, zatytuło- wana „Identyfikacja kinematycznych modeli ruchu stawu kolanowego na podsta- wie sekwencji obrazów RTG” oraz cała zaprezentowana w niej praca została wyko- nana osobi´scie i samodzielnie. Dodatkowo o´swiadczam, ˙ze:

• badania przedstawione w niniejszej rozprawie zostały przeprowadzone w ca- ło´sci w czasie moich studiów doktoranckich oraz działalno´sci naukowej na Politechnice Pozna ´nskiej

• ˙zadna z cz˛e´sci niniejszej rozprawy nie była wcze´sniej opublikowana w pracach kwalifikuj ˛ acych do nadania tytułu naukowego

• nie korzystałam ze ´zródeł innych ni ˙z wymienione w pracy

• cz˛e´sci niniejszej rozprawy, które powstały w wyniku kolaboracji z innymi oso- bami s ˛ a dokładnie oznaczone, a nazwiska wszystkich zaanga ˙zowanych osób zostały wymienione

Podpis:

Data:

v

POLITECHNIKA POZNA ´NSKA

Streszczenie

Instytut Automatyki i Robotyki Wydział Informatyki

Identyfikacja kinematycznych modeli ruchu stawu kolanowego na podstawie sekwencji obrazów RTG

Marta DR ˛A ˙ZKOWSKA

Niniejsza rozprawa ma na celu przedstawienie opracowanej metody identyfi- kacji modelu stawu kolanowego M , który umo˙zliwia opis kinematyki stawów na ró ˙znym etapie rozwoju, zarówno patologicznych jak i zdrowych. Zaproponowano podział modelu kinematyki na dwa podmodele: pozycji f

oraz orientacji f

. Po- prawny jest zatem zapis M = f

( _ ) , f

( _ )
. Dekompozycja jest mo ˙zliwa dzi˛eki przedstawionemu w rozprawie rozumowaniu, w którym orientacj˛e mo ˙zna jedno- znacznie wyznaczy´c dla danej pozycji stawu kolanowego.

Jako dane wej´sciowe wykorzystano zbiór zdj˛e´c rentgenowskich zarejestrowa- nych podczas wykonywania pasywnego ruchu zginania i/lub wyprostu. Zdj˛ecia dotycz ˛ a kolan patologicznych młodych pacjentów, u których proces kostnienia nie jest jeszcze uko ´nczony. Powoduje to znaczne trudno´sci zwi ˛ azane z analiz ˛ a obrazów i wykrywaniem ich cech charakterystycznych. W zwi ˛ azku z powy ˙zszym istotne wydaje si˛e utworzenie dedykowanego algorytmu identyfikacji modelu kinematyki stawu kolanowego.

Zaproponowany algorytm składa si˛e z nast˛epuj ˛ acych etapów: wyboru anato- micznych cech charakterystycznych ko´sci stawu kolanowego, wykrycia na obrazach punktów kluczowych jednoznacznie opisuj ˛ acych te cechy, wyznaczenia konfiguracji stawu kolanowego, wyboru struktury modelu oraz wyznaczenia parametrów mini- malizuj ˛ acych przyj˛ete kryterium, walidacji i analizy porównawczej w ´swietle ist- niej ˛ acych rozwi ˛ aza ´n. Zakłada si˛e, ˙ze algorytm b˛edzie wymagał wiedzy eksperckiej specjalistów z dziedziny ortopedii jedynie w fazie projektowania, a pozostałe etapy b˛ed ˛ a wykonywane automatycznie.

Praca prezentuje interdyscyplinarne podej´scie do rozwi ˛ azania problemu. Wyko- rzystano m.in. techniki z dziedziny przetwarzania obrazów, ł ˛ acznie z wykorzysta- niem splotowych sieci neuronowych, metody identyfikacji modelu, metody staty- styczne. Podkre´slenia wymaga fakt, ˙ze proces doboru struktury modelu jest opty- malizowany w oparciu o heurystyczny algorytm przeszukiwania grafu wa ˙zonego.

Zaproponowano równie ˙z przeprowadzenie optymalizacji estymatora pozycji punk- tów kluczowych na obrazach, dzi˛eki czemu mo ˙zliwa jest analiza setek struktur kan- dyduj ˛ acych. Kompleksowe podej´scie do zaprezentowanego rozwi ˛ azania powoduje,

˙ze algorytm wymaga strojenia bardzo małej liczby parametrów a wynikowe modele z du ˙z ˛ a dokładno´sci ˛ a odwzorowuj ˛ a geometryczne zale ˙zno´sci zachodz ˛ ace w stawie kolanowym.

Przedstawiona metoda wpisuje si˛e w trend wykorzystania narz˛edzi automatyki i robotyki do usprawnienia procesu diagnozy i leczenia opartych o dane graficzne.

A to z kolei wynika z bezpo´sredniej potrzeby personelu medycznego.

POZNAN UNIVERSITY OF TECHNOLOGY

Abstract

Institute of Automation and Robotics Faculty of Computing

Identification of kinematic models for the knee joint based on X-ray image sequences

Marta DR ˛A ˙ZKOWSKA

The aim of this work is to present the designed model M identification method of knee joint, which is able to describe the kinematics of joints at different maturity level, the normal as well as the pathological ones. It is proposed to divide the ki- nematic model into two sub-models: position f

and orientation f

. Therefore, the following relation can be given M = f

( _ ) , f

( _ )
. The decomposition is possible due to the presented idea, in which the orientation can be derived for a given posi- tion of the knee joint.

The sets of X-ray images gathered during passive flexion/extension movement are used as initial data. The images are of pathological knees of young children, for whom the ossification is not accomplished. It raises difficulties connected with image analysis and feature detection. Due to the aforementioned reasons the need to introduce a tailored knee joint kinematic model identification algorithm arises.

The proposed algorithm consists of the following steps: selection of anatomi- cal features of knee joint bones, detection of image points unambiguously denoting those features, estimation of knee joint configuration, selection of model structure and calculation of parameter values that minimise the chosen criterion, validation and comparison with existing methods. It is assumed, that the algorithm will requ- ire the orthopaedic specialists contribution only in design stage and the rest of the steps will be conducted automatically.

The work presents the interdisciplinary approach towards the solution of a stated problem. Among others, the following techniques had been applied: image proces- sing, including convolutional neural networks, model identification methods, sta- tistical methods. It is worth pointing out, that the process of the model structure selection is optimized by the heuristic graph search algorithm. The optimization of characteristic points estimator is proposed as well, allowing the analysis of do- zens of candidate structures. The complexity of the presented approach causes the limitation in the number of design parameters, and the resulting models are highly accurate.

The presented method fits within the trend of transforming the methods used in automation and robotics to improve image based medical diagnostics and treatment.

And the latter arise from the direct medical personnel needs.

ix

Podzi˛ekowanie

W ci ˛ agu ostatnich lat uzyskałam wsparcie wielu ˙zyczliwych mi osób i to dzi˛eki nim udało mi si˛e uko ´nczy´c prace, których zwie ´nczeniem jest niniejsza rozprawa doktorska.

Chciałabym przede wszystkim podzi˛ekowa´c moim Rodzicom, których udziału w kształtowaniu mnie jako osoby, któr ˛ a jestem dzisiaj, nie mo ˙zna podwa ˙zy´c. Moi kochani Rodzice, udzielili mi nieocenionej pomocy w najtrudniejszych momentach i bez w ˛ atpienia, bez ich wkładu, nie udało by mi si˛e zako ´nczy´c tego etapu mojego

˙zycia.

Szczególne podzi˛ekowanie kieruj˛e równie ˙z do dr n. med. Pawła Koczewskiego za po´swi˛ecenie wielu godzin na konsultacje oraz nieust˛epliwo´s´c w poszukiwaniu poprawnego rozwi ˛ azania, pomimo istnienia wielu fałszywych ´zródeł. Dzi˛ekuj˛e rów- nie ˙z dr n. med. Miludowi Shadi, za zgromadzenie wraz z dr Koczewskim niezb˛ed- nych danych do przeprowadzenia prac zebranych w tej rozprawie.

Chciałabym równie ˙z podzi˛ekowa´c Kole ˙zankom i Kolegom z Instytutu Automa- tyki i Robotyki, za codzienne wsparcie i motywacj˛e do pracy. Na szczególne wy- razy uznania zasługuje mgr in ˙z. Tomasz Gawron, który pomógł mi przezwyci˛e-

˙zy´c kryzys naukowy jakiego do´swiadczyłam podczas moich studiów doktoranc- kich. Chciałabym podzi˛ekowa´c równie ˙z dr hab. in ˙z. Dariuszowi Pazderskiemu oraz dr hab. in ˙z. Maciejowi M. Michałkowi za pomoc, której mogły mi udzieli´c je- dynie osoby z ich poziomem do´swiadczenia i wiedzy. Szczególne podzi˛ekowania kieruj˛e równie ˙z do Marcina Nowickiego, Mateusza Przybyły, Rafała Mado ´nskiego, Krzysztofa Łakomego i Joanny Gaw˛eckiej za du ˙z ˛ a dawk˛e optymizmu i motywacji.

Szczególne wyrazy podzi˛ekowania nale ˙z ˛ a si˛e równie ˙z mojemu Promotorowi,

Panu Profesorowi Krzysztofowi Kozłowskiemu, za po´swi˛econy czas oraz wsparcie

naukowe.

xi

Spis tre´sci

Streszczenie v

1 Wprowadzenie 1

1.1 Motywacja i przedstawienie problemu . . . . 1

1.2 Opis istniej ˛ acych rozwi ˛ aza ´n . . . . 8

1.3 Teza, cel i zakres pracy . . . . 9

2 Stan wiedzy 13

2.1 Model kinematyczny stawu kolanowego . . . 13

2.2 Metody identyfikacji modelu . . . 16

2.3 Analiza obrazów medycznych . . . 18

3 Wst˛epne przetwarzanie danych 23

3.1 Selekcja anatomicznych cech ko´sci stawu kolanowego . . . 23

3.2 Wykrywanie punktów kluczowych adaptacyjnym progowaniem . . . 31

3.2.1 Filtracja obrazu RTG . . . 31

3.2.2 Dedykowany algorytm wykrywania kraw˛edzi ko´sci . . . 36

3.2.3 Wyniki estymacji pozycji punktów kluczowych LA . . . 37

3.3 Wykrywanie punktów kluczowych w oparciu o CNN . . . 40

3.3.1 Struktura sieci i regularyzacja . . . 41

3.3.2 Przygotowanie zbiorów ucz ˛ acych . . . 44

3.3.3 Optymalizacja struktury estymatora . . . 47

3.3.4 Wyniki procesu uczenia i walidacji modelu . . . 48

4 Wyznaczenie konfiguracji stawu kolanowego 55

4.1 Wyznaczenie konfiguracji ko´sci . . . 55

4.2 Nakładanie obrazów z wykorzystaniem konfiguracji ko´sci udowej . . 60

5 Identyfikacja modelu kinematyki 67

5.1 Identyfikacja z przyj˛et ˛ a struktur ˛ a modelu . . . 68

5.2 Identyfikacja z szukaniem struktury modelu . . . 74

5.2.1 Rozpatrywane struktury modelu . . . 74

5.2.2 Algorytm wyszukiwania struktury modelu . . . 80

5.2.3 Wyznaczone struktury modelu . . . 82

6 Walidacja i analiza porównawcza modeli 93

6.1 Walidacja zaproponowanych modeli . . . 93

6.2 Porównanie z modelami znanymi z literatury . . . 98

7 Podsumowanie 103

A Wyznaczanie parametrów elipsy 107

B Uwagi implementacyjne 109

Bibliografia 115

xiii

Spis rysunków

1.1 Pierwszy obraz z sekwencji dla pacjenta w wieku 12 lat, przedsta- wiaj ˛ acy staw kolanowy w pełnym wypro´scie, oraz ostatni obraz z se- kwencji przedstawiaj ˛ acy pełne zgi˛ecie. . . . 2 1.2 Oznaczenie układów współrz˛ednych opisuj ˛ acych konfiguracje ko´sci

na obrazach fluoroskopowych oraz przykładowej chwilowej osi ob- rotu (ICR, z ang. Instantaneous Center of Rotation) stawu kolanowego. . 3 1.3 Graficzna reprezentacja przykładowej struktury podmodelu pozycji,

jako krzywej poziomicowej paraboloidy eliptycznej. . . . 6 2.1 Przykładowe modele kinematyczne stawu kolanowego: (a) model za-

kładaj ˛ acy ruch w 6DOF, (b) model zakładaj ˛ acy ruch w 2DOF. . . 14 3.1 Stanowisko do przeprowadzenia bada ´n: (a) zdj˛ecie aparatu fluoro-

skopowego [108], (b) graficzna reprezentacja przebiegu badania. . . 24 3.2 Zdj˛ecie RTG stawu kolanowego [85]: (a) półrocznego dziecka, (b) czte-

rolatka, (c) dziesi˛eciolatka, (d) pietnastolatka. . . 26 3.3 Przykładowa sekwencja X obrazów RTG odpowiadaj ˛ aca pacjentowi P1,

zgodnie z Tab. 3.1. . . 27 3.4 Wybrane cechy kluczowe ko´sci udowej i piszczelowej. Oryginalne

zdj˛ecie zostało poddane prostym operacjom przetwarzania obrazu w celu wizualizacji. . . 28 3.5 Punkty kluczowe obu ko´sci widocznych na obrazach RTG: (a) punkty

r˛ecznie naniesione na obrazy, (b) referencyjne punkty kluczowe, zgod- nie z opisem zawartym w Tab. 3.2. . . 29 3.6 Obraz RTG z oznaczeniem soczewki fluoroskopu (wi˛ekszy okr ˛ ag) oraz

okr˛egu poszukiwania kraw˛edzi ko´sci (mniejszy okr ˛ ag). . . 32 3.7 Przykładowe histogramy obrazów z tej samej sekwencji: (a) pierwszy

obraz, odpowiadaj ˛ acy pełnemu wyprostowi, (b) ostatni obraz, odpo- wiadaj ˛ acy pełnemu zgi˛eciu. . . 33 3.8 Obraz RTG, któremu odpowiada histogram z Rys. 3.7(b). Obraz został

przetworzony w celu wizualizacji. . . 33 3.9 Poszczególne kroki Algorytmu 1: (a) obraz oryginalny X

, (b) obraz

przefiltrowany r X

, (c) obraz X

(rozja´sniony w celu wizualizacji), (d) obraz wynikowy p X

(przedstawiony w formie negatywu oryginału w celu wizualizacji). . . 34 3.10 Wynik jednego z kroków algorytmu opartego o adaptacyjne progo-

wanie. Obraz ˘ X

uzyskiwany jest przez usuwanie nieistotnych obiek- tów z obrazu przedstawionego na Rys. 3.9(d). Przedstawiony w for- mie negatywu oryginału w celu wizualizacji. . . 36 3.11 Wyniki estymacji orientacji LA dla ko´sci udowej. Opis pacjentów jest

zgodny z Tab. 3.1. . . 38

xiv

3.12 Przykładowe wyniki estymacji pozycji k ˛ atowej LA ko´sci udowej dla poszczególnych ramek sekwencji. Opis pacjentów jest zgodny z Tab. 3.1. 39 3.13 Wynikowe obrazy po poszczególnych krokach algorytmu wykrywa-

nia punktów kluczowych LA: (a) obrazy po adaptacyjnym progowa- niu, (b) obrazy po operacjach morfologicznych. . . 40 3.14 Przyj˛eta struktura splotowej sieci neuronowej ze splotowym blokiem

B

oraz blokiem g˛esto poł ˛ aczonym B

. . . 42 3.15 Operacja wycinania okna o rozmiarze 178px178px z obrazu RTG o

oryginalnym rozmiarze 536px536px. Obraz RTG został przetwo- rzony w celu wizualizacji. . . 46 3.16 Proces normalizacji danych wej´sciowych [55]. . . . 46 3.17 Warto´sci wyj´scia funkcji aktywacji: (a) ReLU, (b) Leaky ReLU. . . 51 3.18 Przykładowe wyniki graficzne estymatora pozycji punktów kluczo-

wych PS: białe znaczniki - r˛eczne oznaczenie, czerwone znaczniki - wynik estymacji sieci. Obrazy zostały przetworzone w celu wizuali- zacji. . . 52 3.19 Warto´sci funkcji straty dla zbioru walidacyjnego [26]: (a) histogram

funkcji L wraz z rozkładem normalnym, (b) rozkłady normalne funk- cji L dla poszczególnych punktów kluczowych PS. . . 53 3.20 Przykładowe wyniki graficzne sieci neuronowej z du ˙z ˛ a warto´sci ˛ a funk-

cji straty: białe znaczniki - r˛eczne oznaczenie, czerwone znaczniki - estymacja sieci. Obrazy zostały przetworzone w celu wizualizacji. . . . 54 4.1 Oznaczenie układu współrz˛ednych okre´slaj ˛ acego konfiguracj˛e ko´sci

udowej. Obraz został przetworzony w celach wizualizacji. . . 56 4.2 Oznaczenie układu współrz˛ednych okre´slaj ˛ acego konfiguracj˛e ko´sci

piszczelowej. . . 57 4.3 Bł˛edy estymacji ka ˙zdej zmiennej konfiguracyjnej ko´sci udowej na ob-

razach ze zbioru walidacji, histogram wraz z rozkładem normalnym. . 60 4.4 Graficzna reprezentacja tworzenia dodatkowych punktów kluczowych

ko´sci udowej zebranych w macierzy K

. Obraz został przetworzony w celu wizualizacji. . . 63 4.5 Bezwzgl˛edna warto´s´c bł˛edu k ˛ atowego orientacji ko´sci udowej w funk-

cji numeru obrazu z sekwencji: (a) wyniki dla pacjenta P2, (b) wyniki dla pacjenta P13. . . 65 4.6 Odległo´s´c Φ: (a) wyniki dla pacjenta P2, (b) wyniki dla pacjenta P13. . 65 5.1 Graficzna reprezentacja zadania identyfikacji podmodelu pozycji, dla

stałej g

. . . 67 5.2 Porównanie przybli ˙zenia okr˛egiem i elips ˛ a zbioru punktów [ _rx

y _r

]

dla pacjenta P1. . . 70 5.3 Graficzna wizualizacja zaproponowanego podmodelu orientacji wraz

z oznaczonymi konfiguracjami stawu kolanowego dla pacjenta P1. . . 72 5.4 Graficzna reprezentacja wej´s´c i wyj´s´c poszczególnych podmodeli wy-

ra ˙zonych za pomoc ˛ a sieci neuronowych: (a) podmodel pozycji, (b) pod- model orientacji. . . 76 5.5 Graficzna wizualizacja zadania interpolacji funkcji. Czarnymi znacz-

nikami oznaczono w˛ezły, tj. warto´sci odczytane z danych wej´scio-

wych, natomiast czerwonymi znacznikami warto´sci mi˛edzy-w˛ezłowe. 77

xv

5.6 Przykładowe grafy wraz z oznaczeniem struktury danego w˛ezła, skom- plikowania i ´sredniego bł˛edu testowego oraz kosztu przej´scia do s ˛ a- siadów: (a) modele analityczne, (b) struktury sieci neuronowej dla przykładowo przyj˛etego ξ = 62. Przedstawione warto´sci ¯J nie odpo- wiadaj ˛ a rzeczywistym warto´sciom. . . 80 5.7 Przykładowe krzywe opisuj ˛ ace podmodel pozycji, odpowiadaj ˛ ace pa-

cjentowi P4, z oznaczon ˛ a konfiguracj ˛ a Q

. . . 83 5.8 Przykładowe struktury podmodelu orientacji dla pacjenta P4: (a) kra-

jobraz dla struktury S

( 4, 3 ) , (b) krajobraz dla struktury S

( 1, 2 ) . . . . 85 5.9 Histogram ró ˙znicy pomi˛edzy oczekiwanym a rzeczywistym wyj´sciem

z sieci dla podmodelu pozycji: (a) porównanie wyniku dla S

( 10, 11, 0, 0 ) oraz S

( 25, 31, 8, 0 ) , (b) wynik dla S

( 25, 31, 8, 0 ) z podziałem na trzy zbiory ucz ˛ ace. . . 87 5.10 Warto´sci MSE ka ˙zdego ze zbiorów ucz ˛ acych podmodelu pozycji w

funkcji numeru epoki uczenia. . . . 88 5.11 Przykładowa sekwencja odpowiadaj ˛ aca pacjentowi P2: oczekiwane i

rzeczywiste wyj´scie sieci. . . 88 5.12 Histogram ró ˙znicy pomi˛edzy rzeczywistym a oczekiwanym wyj´sciem

z sieci dla podmodelu orientacji: (a) wynik porównawczy, (b) wynik z podziałem na zbiory ucz ˛ ace dla S

( 17, 11, 4, 4 ) . . . 90 5.13 Warto´sci MSE ka ˙zdego ze zbiorów ucz ˛ acych podmodelu orientacji

S

( 17, 11, 4, 4 ) w funkcji numeru epoki uczenia. . . 90 5.14 Przykładowa sekwencja odpowiadaj ˛ aca pacjentowi P2: oczekiwane

i rzeczywiste wyj´scie sieci dla podmodelu orientacji oraz krajobraz wyj´scia sieci. . . 91 6.1 Warto´sci RMSE dla ka ˙zdej z analizowanych struktur modelu kinema-

tycznego stawu kolanowego, z podziałem na dwie współrz˛edne mo- delu. . . 94 6.2 Porównanie obu współrz˛ednych modelu M

podczas ruchów wyko-

nanych przed oraz po przeprowadzonej operacji: (a) pacjent P10 pod- czas ruchu wyprostu, (b) pacjent P8 podczas ruchu zgi˛ecia. . . 97 6.3 Porównanie obu współrz˛ednych modelu M

dla pacjenta P5 podczas

dwóch typów ruchu wykonanych przed przeprowadzon ˛ a operacj ˛ a. . 98 6.4 RMSE dla wszystkich analizowanych kinematycznych modeli z po-

działem na dwie współrz˛edne modelu. . . 100 6.5 Wyznaczone ´srodki okr˛egów modelu 4bar, dla ró ˙znych punktów klu-

czowych ko´sci, dla pacjenta P4. . . 101 6.6 Porównanie obu współrz˛ednych modeli dla pacjenta: (a) P1 (najmłod-

szy z grupy), (b) P15 (najstarszy z grupy). . . 101 B.1 J ˛ adrowy estymator g˛esto´sci prawdopodobie ´nstwa jednowymiarowej

zmiennej [62]. . . 110 B.2 Przykładowe rozkłady prawdopodobie ´nstwa G i Z z oznaczeniem

wylosowanej grupy nowych kandydatów. . . . 110

xvii

Spis tablic

3.1 Dane pacjentów, dla których zebrano sekwencje obrazów RTG. Na czerwono oznaczono sekwencje, które zarejestrowano po przeprowa- dzeniu operacji. . . 25 3.2 Wybrane cechy kluczowe ko´sci, odpowiadaj ˛ ace im r˛ecznie oznaczone

punkty oraz referencyjne punkty kluczowe. . . 28 3.3 Rozmiary okien adaptacyjnego progowania dla poszczególnych pa-

cjentów. Opis pacjentów jest zgodny z Tab. 3.1. . . 39 3.4 Liczebno´s´c zbiorów ucz ˛ acych przygotowanych dla CNN. . . 44 3.5 Hiperparametry CNN oraz ich rozwa ˙zane warto´sci. . . 47 4.1 Wykorzystanie ró ˙znych typów konfiguracji stawu kolanowego, jako

danych wej´sciowych do poszczególnych zada ´n. . . 60 5.1 Warto´sci bł˛edu MSE dla wszystkich badanych struktur podmodelu

pozycji dla identyfikacji z przyj˛et ˛ a struktur ˛ a modelu. . . 69 5.2 Dobrane warto´sci podmodelu pozycji dla identyfikacji parametrycz-

nej. . . 71 5.3 Warto´sci MSE podmodelu orientacji dla identyfikacji parametrycznej. 73 5.4 Wyznaczone warto´sci parametrów podmodelu orientacji dla identy-

fikacji parametrycznej. Dane wej´sciowe zostały uprzednio skalowane. 73 5.5 Liczebno´s´c zbiorów ucz ˛ acych przygotowanych dla sieci neuronowej. . 79 5.6 Porównanie struktur podmodelu pozycji dla wszystkich rozwa ˙zanych

zestawów konfiguracji. . . 83 5.7 Porównanie przykładowych struktur podmodelu orientacji dla wszyst-

kich rozwa ˙zanych zestawów konfiguracji. . . 84 5.8 Przykładowe struktury sieci opisuj ˛ acej podmodel pozycji dla identy-

fikacji z wyznaczaniem struktury. . . . 86 5.9 Przykładowe struktury sieci opisuj ˛ acej podmodel orientacji dla iden-

tyfikacji z wyznaczaniem struktury. . . 89 6.1 Warto´sci skomplikowania, ´srednie bł˛edy zbioru testowego oraz war-

to´sci przyj˛etych statystyk dla omawianych modeli. . . 95 6.2 Warto´sci opisuj ˛ ace prymitywy w predykcyjnych modelach znanych z

literatury [52,

B.1 Parametry maszyn pracuj ˛ acych w klastrze, wykorzystanych w proce-

sie poszukiwania optymalnej struktury CNN. . . 111

xix

Lista skrótów

LA Long Axis - o´s ´srodkowo-trzonowa

PS Patellar Surface - wci˛ecie mi˛edzykłykciowe GP Growth Plate - chrz ˛

astka wzrostowa

ICR Instantaneous Center of Rotation - chwilowa o´s obrotu IHA Instantaneous Helical Axis - ´srubowa o´s obrotu

MRI Magnetic Resonance Imaging - obrazowanie metod ˛

a rezonansu magnetycznego

- zdj˛ecie RenTGenowskie

CT Computed Tomography - tomografia komputerowa

SIFT Scale-Invariant Feature Transform - skalo-niezmiennicze przekształcenie cech DOF Degree Of Freedom - stopie ´n swobody

LS Least Squares - metoda najmniejszych kwadratów

MLE Maximum-Likelihood Estimation - metoda estymacji najwi˛ekszej wiarygodno´sci CNN Convolutional Neural Network - splotowa sie´c neuronowa

TPE Tree-structured Parzen Estimator - drzewo estymatorów Parzen PCA Principal Component Analysis - analiza głównych składowych ICP Iterative Closest Point - iteracyjny najbli ˙zszy punkt

MSE Mean Squared Error - bł ˛

ad ´sredniokwadratowy

RMSE Root Mean Square Error - pierwiastek bł˛edu ´sredniokwadratowego BN Batch Normalization - normalizacja partii

DO DropOut - wygaszanie neuronów

Funkcje aktywacji sieci neuronowych:

ReLU Rectified Linear Unit - funkcja liniowa zwracaj ˛

aca 0 dla ujemnych warto´sci

LR Leaky ReLU - modyfikacja ReLU z niezerowymi wyj´sciami dla ujemnych warto´sci SELU Scaled Exponential Linear Unit - wykładnicza modyfikacja ReLU

SIG SIGmoid - funkcja logistyczna

HS Hard Sigmoid - funkcja liniowa z nasyceniem

xxi

Lista symboli

K

macierz punktów kluczowych powi ˛ azana z i-tym obrazem sekwencji K macierz punktów kluczowych odpowiadaj ˛ aca całej sekwencji

k

j-ty punkt kluczowy

d

konfiguracja ko´sci piszczelowej na i-tym obrazie sekwencji, zło ˙zona z:

θ_di

orientacji ko´sci piszczelowej d ¯

pozycji ko´sci piszczelowej

g

konfiguracja ko´sci udowej na i-tym obrazie sekwencji, zło ˙zona z:

θg_i

orientacji ko´sci udowej

¯g

pozycji ko´sci udowej

φ_d

funkcja przekształcaj ˛ aca K

na konfiguracj˛e ko´sci piszczelowej

funkcja przekształcaj ˛ aca K

na konfiguracj˛e ko´sci udowej q

konfiguracja stawu kolanowego na i-tym obrazie, zło ˙zona z:

θ_qi

orientacji stawu kolanowego

x

współrz˛ednej x pozycji stawu kolanowego y

współrz˛ednej y pozycji stawu kolanowego f

( _ ) podmodel pozycji

f

( _ ) podmodel orientacji

p

optymalne parametry podmodelu pozycji p

optymalne parametry podmodelu orientacji J ( p

) bł ˛ ad identyfikacji podmodelu pozycji J ( p

) bł ˛ ad identyfikacji podmodelu orientacji M kinematyczny model stawu kolanowego

ˆ

x

aproksymacja x

wyznaczona przez podmodel pozycji ˆ

y

aproksymacja y

wyznaczona przez podmodel pozycji θ ˆ

aproksymacja θ

wyznaczona przez podmodel orientacji

m ´srodek geometryczny wci˛ecia mi˛edzykłykciowego ko´sci udowej

σ

odchylenie standardowe

K

rozszerzona macierz punktów kluczowych ko´sci udowej obrazu i-tego k

dodatkowe punkty kluczowe, równomiernie rozło ˙zone na wci˛eciu

mi˛edzykłykciowym ko´sci udowej

k

dodatkowe punkty kluczowe, równomiernie rozło ˙zone na osi ko´sci udowej

||  || norma euklidesowa

Konfiguracje stawu kolanowego dla całej sekwencji obrazów

Q ogólna konfiguracja stawu kolanowego, zło ˙zona z:

θ

orientacji stawu kolanowego x

, y

pozycji stawu kolanowego

Q

konfiguracja wyznaczona na podstawie r˛ecznie oznaczonych punktów, zło ˙zona z:

θ

orientacji stawu kolanowego

x

, y

pozycji stawu kolanowego

xxii

Q

konfiguracja wyznaczona na podstawie estymowanych punktów, zło-

˙zona z:

θ

orientacji stawu kolanowego x

y

pozycji stawu kolanowego

Q

konfiguracja wyznaczona na podstawie r˛ecznie oznaczonych punktów z dodatkowym szumem, zło ˙zona z:

θ

orientacji stawu kolanowego x

y

pozycji stawu kolanowego

Q r konfiguracja wyznaczona dla sekwencji obrazów z r˛ecznie nało ˙zon ˛ a pozycj ˛ a ko´sci udowej, zło ˙zona z:

θ r orientacji stawu kolanowego rx r y pozycji stawu kolanowego

Q r

konfiguracja wyznaczona przez dodanie szumu do r Q, zło ˙zona z:

θ r

orientacji stawu kolanowego rx

y r

pozycji stawu kolanowego

X sekwencja obrazów

S liczba obrazów w sekwencji X

i-ty obraz sekwencji

s rozmiar okna adaptacyjnego progowania

dodatkowy próg adaptacyjnego progowania X r

wynik filtracji obrazu X

filtrem u´sredniaj ˛ acym X

obraz pomocniczy adaptacyjnego progowania X p

obraz wyj´sciowy adaptacyjnego progowania

X ˘

obraz uzyskiwany przez usuwanie nieistotnych obiektów z obrazu p X

liczba wykrytych punktów na kraw˛edziach ko´sci, parametr algorytmu

adaptacyjnego progowania

c warto´s´c kontrastu wykorzystana w procesie rozszerzania zbiorów SAT funkcja saturacji

I

obraz wyj´sciowy, po zmianie kontrastu I

obraz wej´sciowy, przed zmian ˛ a kontrastu T zbiór macierzy transformacji dla całej sekwencji

` operator nakładania transformacji na macierz punktów kluczowych Φ kryterium odległo´sci dla zadania nakładania obrazów

T

( _ ) macierz transformacji i-tego zdj˛ecia zło ˙zona z rotacji i translacji o na- st˛epuj ˛ acych warto´sciach parametrów:

θ_t

k ˛ at obrotu

x

translacja w poziomie y

translacja w pionie

B

j-ty blok sieci neuronowej

liczba bloków sieci neuronowej D g˛esto poł ˛ aczona warstwa sieci Conv splotowa warstwa sieci N

liczba warstw w j-tym bloku

liczba bloków splotowych

p

ϕ estymowane wyj´scie sieci neuronowej

ϕ oczekiwane wyj´scie sieci neuronowej

xxiii

L funkcja straty

X partia w procesie uczenia sieci

n

liczba neuronów w warstwie g˛esto poł ˛ aczonej n liczba neuronów w warstwach D

η_j

liczba filtrów w warstwie splotowej j-tego bloku w

rozmiar okna splotu j-tego bloku

p

rozmiar okna redukcji j-tego bloku

max ( _ ) redukcja przez wybór warto´sci maksymalnej mean ( _ ) redukcja przez wybór warto´sci ´sredniej d

prawdopodobie ´nstwo dropoutu j-tego bloku

$

funkcja aktywacji

α

k ˛ at nachylenia dla ujemnych warto´sci dla funkcji aktywacji Leaky ReLU

ñ operacja spłaszczania danych do wektora kolumnowego

, β

parametry uczone normalizacji partii

γ₁

, γ

warto´sci progowe funkcji straty dla kryterium wczesnego zatrzymania E ˆ ( _ ) warto´s´c oczekiwana zmiennej

Var y ( _ ) wariancja zmiennej

u liczba hiperparametrów sieci

M unikalny zestawem hiperparametrów CNN k iteracja algorytmu optymalizacji struktury CNN

G populacja modeli CNN gwarantuj ˛ aca nisk ˛ a warto´s´c funkcji starty Z populacja modeli CNN z wysok ˛ a warto´sci ˛ a funkcji straty

EI ( _ ) współczynnik spodziewanej poprawy

S

( _ ) struktura podmodelu pozycji S

( _ ) struktura podmodelu orientacji

próg warto´sci funkcji kosztu

r stopie ´n wielomianu, dla wielomianów pełnej reprezentacji r

, r

pot˛ega zmiennych wielomianu w niepełnej reprezentacji N liczba warstw g˛esto poł ˛ aczonych

C ( _ ) skomplikowanie modelu

¯J ( p ) ´srednia warto´s´c J ( p

) dla zbioru testowego C

skomplikowanie j-tego w˛ezła

ζ_i

funkcjonał dla i-tego w˛ezła algorytmu A*

h

funkcja heurystyczna algorytmu A*

U zbiór otwarty algorytmu A*

F zbiór zamkni˛ety algorytmu A*

S

pocz ˛ atkowa struktura podmodelu wykorzystana w algorytmie A*

Z nast˛epnik aktualnej struktury modelu, wykorzystany w algorytmie A*

W

i-ta statystyka uwzgl˛edniaj ˛ aca zło ˙zono´s´c modelu

` parametr skierowany ´scie ˙zki stawu kolanowego, wykorzystany w strukturach opartych o sieci neuronowe

M ilo´s´c zebranych danych wykorzystanych do identyfikacji

F ( _{, } ) rozkład F prawdopodobie ´nstwa zmiennej losowej o stopniach swo- body okre´slonych parametrami wej´sciowymi

ρ

statystyka z rozkładem F, wykorzystana przy testach zło ˙zono´sci mo- delu

ξ_q

długo´s´c wektora wej´s´c sieci neuronowej podmodelu pozycji

długo´s´c wektora wej´s´c sieci neuronowej podmodelu orientacji

xxv

Prac˛e t ˛ a dedykuj˛e mojej Rodzinie (przez wielkie R)

1

Rozdział 1

Wprowadzenie

1.1 Motywacja i przedstawienie problemu

Wraz z post˛epem medycyny oraz wzrostem mo ˙zliwo´sci urz ˛ adze ´n pomiarowych, które mo ˙zna bezpiecznie zastosowa´c do bada ´n na pacjentach, zauwa ˙zalny jest rów- nie ˙z wzrost zainteresowania opisem zjawisk zachodz ˛ acych w ludzkim ciele. Jed- nym z obszarów tego zainteresowania s ˛ a badania nad opracowaniem modelu kine- matycznego poszczególnych stawów. Jako model kinematyczny w niniejszej pracy rozumiane b˛ed ˛ a zachodz ˛ ace w danym stawie zwi ˛ azki geometryczne. Wobec tego model kinematyczny jest zdefiniowany w przestrzeni współrz˛ednych geometrycz- nych, a nie pr˛edko´sci. Opracowanie odpowiednio dokładnego modelu kinema- tycznego stawu jest istotne ze wzgl˛edów u ˙zytkowych jak i poznawczych. Spo´sród tych pierwszych, wymieni´c mo ˙zna przykładowo usprawnienie procesu rehabilitacji (współczesne metody oraz zastosowane w nich uproszczenia mog ˛ a by´c niewystar- czaj ˛ ace), umo ˙zliwienie porównania zmian zachodz ˛ acych dla jednego pacjenta (m.in.

przed i po operacji, w trakcie trwania rehabilitacji), mo ˙zliwo´s´c wyznaczenia ilo´scio- wych miar podobie ´nstwa słu ˙z ˛ acych do klasyfikacji pacjentów.

W niniejszej pracy skupiono si˛e na analizie modeli kinematycznych stawu kola- nowego, jednak zaprezentowana metodyka mo ˙ze by´c z łatwo´sci ˛ a dostosowana do stawów innego typu. Staw kolanowy jest zło ˙zonym stawem ciała ludzkiego, któ- rego reprezentacja jest cz˛esto nadmiernie upraszczana. W literaturze mo ˙zna znale´z´c ró ˙zne opracowania dotycz ˛ ace modeli kinematycznych kolana zarówno w dwóch wymiarach, jak i w trzech. Modele te dotycz ˛ a głównie osób dorosłych, a grupa ba- danych jest zazwyczaj jednolita pod wzgl˛edem wyst˛epuj ˛ acego schorzenia lub mało liczna. Badania rzadko s ˛ a prowadzone u osób nieletnich, ze wzgl˛edu na konieczno´s´c ograniczenia poziomu napromieniowania, co uniemo ˙zliwia zastosowanie wielu me- tod akwizycji obrazu. Identyfikacja kinematycznego modelu stawu kolanowego dzieci, u których proces kostnienia nie jest jeszcze uko ´nczony, wydaje si˛e by´c bardzo interesuj ˛ aca, a wyniki uzupełniałyby te ju ˙z istniej ˛ ace, dla szerszej grupy wiekowej.

Proces identyfikacji modelu b˛edzie składał si˛e z kilku etapów. Ka ˙zdy z nich zo- stanie teraz ogólnie nakre´slony, bez przytaczania pozycji literaturowych wspieraj ˛ a- cych przyj˛eta metodyk˛e. Szczegółowy opis aktualnego stanu wiedzy przedstawiony zostanie w Roz. 2, natomiast precyzyjny opis ka ˙zdego z kroków identyfikacji znaj- duje si˛e w odpowiadaj ˛ acym mu rozdziale.

Pierwszym krokiem identyfikacji modelu jest opracowanie i przygotowanie eks-

perymentu maj ˛ acego na celu zebranie danych pomiarowych. W omawianym przy-

padku, wi ˛ a ˙ze si˛e to z rejestracj ˛ a ró ˙znych pozycji ko ´nczyny dolnej człowieka. Ze

wzgl˛edu na ograniczenia dotycz ˛ ace protokołu przeprowadzanych bada ´n u dzieci,

2 Rozdział 1. Wprowadzenie

wskazane jest, aby do analizy stawu kolanowego posłu ˙zy´c si˛e wynikami bada ´n ko- niecznymi ze wzgl˛edów medycznych. Podczas wielu operacji chirurgicznych leka- rze wspomagaj ˛ a si˛e technikami obrazowania medycznego. Przykładowo, przy re- konstrukcji wi˛ezadeł lekarze wykorzystuj ˛ a zdj˛ecia rentgenowskie (RTG) lub fluoro- skopowe. Jest to metoda mało inwazyjna, umo ˙zliwiaj ˛ aca rejestracj˛e pełnego zakresu ruchu, a aparaty fluoroskopowe oraz aparaty rentgenowskie s ˛ a szeroko dost˛epne w polskich szpitalach. Z analiz ˛ a zebranych zdj˛e´c zwi ˛ azany jest szereg trudno´sci, których wpływu na dalsze kroki identyfikacji nie mo ˙zna podwa ˙zy´c.

Drugi etap dotyczy odpowiednio dokładnego okre´slenia konfiguracji, tj. pozycji i orientacji, ko´sci udowej oraz piszczelowej

na dost˛epnych obrazach. Metoda wy- znaczania konfiguracji ko´sci na obrazach powinna odpowiada´c typowi dost˛epnych danych pomiarowych. Jak ju ˙z wspomniano wcze´sniej, charakter uzyskanych zdj˛e´c medycznych wpływa bezpo´srednio na trudno´sci ich analizy.

Przykładowe zdj˛ecie z sekwencji przedstawiono na Rys. 1.1. Przede wszystkim, dost˛epne dane s ˛ a dwuwymiarowe, przy czym rejestrowany ruch nast˛epuje w trzech płaszczyznach, zatem widoczny obraz ko´sci stanowi jego projekcj˛e na płaszczyzn˛e ekranu fluoroskopowego. Dodatkowo, zauwa ˙zalny jest wysoki poziom uwydatnie- nia cech na tkankach mi˛ekkich (np. błona maziowa oznaczona na Rys. 1.1), poja- wiaj ˛ a si˛e znacz ˛ ace ró ˙znice w poziomie skostnienia stawu kolanowego oraz zakłóce- nia obrazu (implanty, ´sruby, aparat Ilizarowa

). Ze wzgl˛edu na powy ˙zsze, wykry- wanie cech kluczowych standardowymi metodami analizy obrazu nie daje satysfak- cjonuj ˛ acych wyników dla opisywanych danych. Istotne jest zatem zaproponowanie odpornego algorytmu okre´slania konfiguracji ko´sci na obrazie.

Rzepka Kość udowa

Kość piszczelowa Błona maziowa

Kość strzałkowa

RYSUNEK 1.1: Pierwszy obraz z sekwencji dla pacjenta w wieku 12 lat, przedstawiaj ˛acy staw kolanowy w pełnym wypro´scie, oraz

ostatni obraz z sekwencji przedstawiaj ˛acy pełne zgi˛ecie.

Trzeci etap identyfikacji mo ˙zna opisa´c jako szukanie struktury modelu, który z odpowiedni ˛ a dokładno´sci ˛ a (wymagan ˛ a do analizy danych medycznych) odzwier- ciedli ruch ko´sci piszczelowej wzgl˛edem ko´sci udowej, zakładaj ˛ ac, ˙ze mo ˙zliwe jest

1Zmiany konfiguracji rzepki oraz ko´sci strzałkowej, pomimo ich obecno´sci na analizowanych ob- razach, b˛ed ˛a w niniejszych rozwa ˙zaniach pomini˛ete. Co prawda rzepka wpływa na kinematyk˛e i dy- namik˛e ruchów w stawie kolanowym [29], lecz do opisu ruchu wystarczaj ˛aca jest analiza wzajemnego przemieszczania si˛e ko´sci udowej i piszczelowej.

2Niektórzy spo´sród badanych przechodz ˛a terapi˛e wydłu ˙zania ko´sci aparatem Ilizarowa. Metoda polega na przeci˛eciu ko´sci i stopniowym wydłu ˙zaniu odległo´sci pomi˛edzy jej cz˛e´sciami. Urz ˛adze- nie mocuje si˛e bezpo´srednio do ko´sci za pomoc ˛a specjalnych pr˛etów, które ł ˛acz ˛a si˛e z zewn˛etrznymi pier´scieniami.

1.1. Motywacja i przedstawienie problemu 3

okre´slenie konfiguracji ko´sci na ka ˙zdym obrazie. Dla i-tego obrazu mo ˙zna przyj ˛ a´c,

˙ze g

P IR

opisuje konfiguracj˛e ko´sci udowej a d

P IR

okre´sla konfiguracj˛e ko´sci piszczelowej. Układy współrz˛ednych zwi ˛ azane z konfiguracjami ko´sci oznaczono schematycznie na Rys. 1.2.

RYSUNEK1.2: Oznaczenie układów współrz˛ednych opisuj ˛acych kon- figuracje ko´sci na obrazach fluoroskopowych oraz przykładowej chwilowej osi obrotu (ICR, z ang. Instantaneous Center of Rotation)

stawu kolanowego.

Konfiguracje ko´sci na obrazach dotycz ˛ a ich pozycji i orientacji g

= ^

¯g



, d

= ^

d ¯



, (1.1)

gdzie θ

i θ

okre´slaj ˛ a i-t ˛ a orientacj˛e ko´sci udowej i piszczelowej, a przez ¯g

i ¯ d

ozna- czono i-t ˛ a pozycj˛e ko´sci udowej i piszczelowej, zgodnie z oznaczeniami zamieszczo- nymi na Rys. 1.2. Wówczas i-t ˛ a konfiguracj˛e stawu kolanowego q

mo ˙zna rozumie´c jako konfiguracj˛e ko´sci piszczelowej d

wyra ˙zon ˛ a w układzie współrz˛ednych ko´sci udowej g

. Wobec tego mo ˙zna zapisa´c, ˙ze

q

=





x

y



 =





θ_di

θ

 cos (

) sin (

)

sin (

) _cos (

)



[ d ^¯

¯g

]



 . (1.2)

Sformułowanie problemu

Dla przyj˛etego opisu konfiguracji stawu kolanowego, identyfikacj˛e kinematycz- nych modeli mo ˙zna przedstawi´c jako wyznaczenie parametrów minimalizuj ˛ a- cych sum˛e bł˛edów estymacji, dla zbioru sekwencji obrazów ró ˙znych pacjentów, ze znanymi (oznaczonymi r˛ecznie lub automatycznie) konfiguracjami ko´sci na obrazach. Zakłada si˛e, ˙ze struktura modelu b˛edzie stała dla wszystkich pacjen- tów, jedynie parametry modelu b˛ed ˛ a indywidualne dla ka ˙zdej osoby.

Ka ˙zdy z powy ˙zszych aspektów zostanie pokrótce opisany. Jako sekwencj˛e obra-

zów X b˛edzie rozumiany zbiór S dwuwymiarowych obrazów X

w skali szaro´sci,

4 Rozdział 1. Wprowadzenie

odpowiadaj ˛ acy jednemu typowi ruchu jednego pacjenta

X _tX

_u

, X

P IR

, (1.3) gdzie obraz X

jest przedstawiony w postaci macierzy o warto´sciach odpowiadaj ˛ a- cych jasno´sci poszczególnych pikseli. Warto´s´c 536 jest zwi ˛ azana z rozdzielczo´sci ˛ a wykorzystywanych obrazów, nie stanowi to jednak ˙zadnego ograniczenia metody.

Ka ˙zdy obraz sekwencji przedstawia staw kolanowy dla innego k ˛ ata ugi˛ecia. Se- kwencje przedstawiaj ˛ a dwa rodzaje ruchów, tj.:

• zgi˛ecie - od pełnego wyprostu do pełnego ugi˛ecia, przy czym ka ˙zdy nast˛epny obraz przedstawia staw ugi˛ety pod wi˛ekszym k ˛ atem,

• wyprost - od pełnego ugi˛ecia do pełnego wyprostu, przy czym ka ˙zdy nast˛epny obraz przedstawia staw ugi˛ety pod mniejszym k ˛ atem.

Warto podkre´sli´c, ˙ze kolejno´s´c nie jest istotna w rozwa ˙zanym przypadku. Rozwa-

˙zane b˛ed ˛ a zbiory dyskretnych konfiguracji ko´sci, przy czym ich uszeregowanie nie wpływa na wynikowe postaci modelu.

W celu jednoznacznego opisu konfiguracji ko´sci na i-tym obrazie, wprowadza si˛e f punktów kluczowych k

. Z i-tym obrazem uto ˙zsamiony jest zbiór punktów kluczowych przedstawiony za pomoc ˛ a nast˛epuj ˛ acej macierzy

K

 k

k

. . . k



= ^x

^x

^{. . . x}

y

y

. . . y



, K

P IR

, (1.4) gdzie x

oraz y

opisuj ˛ a współrz˛edne j-tego punktu kluczowego k

. Dla całej se- kwencji definiuje si˛e macierz

K = K

K

. . . K



P IR

. (1.5)

Dla i-tego obrazu sekwencji przyjmuje si˛e nast˛epuj ˛ ace zale ˙zno´sci

g

=

( K

) ,

: IR

Ñ IR

, oraz (1.6) d

=

( K

) ,

: IR

Ñ IR

, (1.7) gdzie zale ˙zno´sci φ

oraz φ

przekształcaj ˛ a macierz punktów kluczowych (1.4) do konfiguracji ko´sci piszczelowej i udowej, danych równaniem (1.1). Zakłada si˛e, ˙ze przekształcenia φ

oraz φ

b˛ed ˛ a takie same dla wszystkich obrazów.

Dla rozwa ˙zanej sekwencji zdefiniowano zbiór chwilowych konfiguracji stawu kolanowego, który mo ˙zna przedstawi´c w postaci macierzowej przez

Q =



 θ

x

y



 ,

θ

= ^

. . . θ

 , x

= x

x

. . . x

 , y

= y

y

. . . y

 .

(1.8) Dla tak przyj˛etej notacji zdefiniowano problem identyfikacji modeli kinematycz- nych stawu kolanowego z podziałem na dwa podproblemy: problem identyfikacji podmodelu orientacji i problem identyfikacji podmodelu pozycji.

Podmodel pozycji mo ˙zna wyrazi´c dla znanych: sekwencji obrazów X , macierzy

punktów kluczowych K , oraz przekształce ´n φ

i φ

. Identyfikacja polega na odnale-

zieniu funkcji aproksymuj ˛ acej f

, która odwzorowuje z przyj˛et ˛ a dokładno´sci ˛ a zbiór

punktów x

y

 i składa si˛e z nast˛epuj ˛ acych kroków:

1.1. Motywacja i przedstawienie problemu 5

1. Wyboru struktury podmodelu pozycji f

( p

) , gdzie przez p

=  p

p

. . . p

 oznaczono wektor parametrów podmodelu pozycji, przy czym p

P IR . Jako C oznaczono liczb˛e parametrów modelu. Ze wzgl˛edu na rozpatrywane struk- tury podmodelu pozycji, zapis ten jest jednoznaczny, gdy ˙z dla danego zestawu parametrów mo ˙zna opisa´c jedynie jedn ˛ a struktur˛e podmodelu.

2. Identyfikacji najlepszych ze wzgl˛edu na przyj˛ete kryterium warto´sci parame- trów p

przez rozwi ˛ azanie zestawu S równa ´n postaci

f

( p

, x

, y

) = _{0, dla i} = _{1 . . . S, (1.9)} które b˛ed ˛ a szczegółowo opisane w Roz. 5, przykładowo równaniem (5.1). Opty- malne warto´sci parametrów zapewniaj ˛ a spełnienie warunku

p

= arg min

pq

J ( p

) , J ( p

) 

¸

f

( p

, x

, y

)

. (1.10)

3. Dla przyj˛etej struktury podmodelu oraz parametrów uzyskanych w poprzed- nim punkcie nale ˙zy wyznaczy´c punkty [ x ˆ

y ˆ

] , które s ˛ a rzutem ortogonalnym punktów [ x

y

] na krzyw ˛ a f

( p

, ˆx

, ˆy

) = 0. Krok ten mo ˙zna przedstawi´c za pomoc ˛ a nast˛epuj ˛ acego schematu

[ x

y

] _Ñ f

( p

, ˆ x

, ˆ y

) = 0 Ñ [ x ˆ

y ˆ

] .

Ò Ò Ò

dyskretny opis ci ˛agły dyskretny

zbiór zbiór

(1.11)

Po zako ´nczonym procesie identyfikacji podmodel pozycji z ju ˙z dobranymi struktur ˛ a i warto´sciami parametrów, przybli ˙za warto´sci pozycji stawu kolano- wego, do wyznaczonej krzywej. Warto podkre´sli´c, ˙ze warto´sci pozycji mog ˛ a by´c spoza zbioru rzeczywistych odczytów z obrazów wej´sciowych. Podmo- del pozycji jest funkcj ˛ a skalarn ˛ a, która zwraca jednowymiarow ˛ a krzyw ˛ a tj.

f

( p

, ,  ) : IR

Ñ IR spełniaj ˛ac ˛a ograniczenie f

( p

, ˆ x

, ˆ y

) = 0.

Podmodel pozycji mo ˙ze by´c rozpatrywany jako krzywa poziomicowa b˛ed ˛ aca wynikiem przeci˛ecia dwóch powierzchni. Przykładowo, gdyby elipsa była wybrana jako struktura podmodelu, wówczas mo ˙zna j ˛ a rozpatrywa´c jako krzyw ˛ a poziomi- cow ˛ a o poziomie równym 0 paraboloidy eliptycznej, jak pokazano na Rys. 1.3. Wów- czas podmodel pozycji b˛edzie wynikiem spełnienia dwóch ogranicze ´n:

# z = 0,

z =

+

_1. ^(1.12) Przedstawiony przykład dotyczy krzywej zamkni˛etej, jednak podmodel pozycji nie jest ograniczony do krzywych tego typu.

Warto zauwa ˙zy´c, ˙ze wybrany funkcjonał J ( p

) jest równy zeru, wtedy i tylko

wtedy, gdy wszystkie punkty [ x

y

] _{le ˙z ˛} a na tej wybranej krzywej uwikłanej. Im

wi˛eksza warto´s´c funkcjonału, tym bardziej oddalone s ˛ a punkty od krzywej opisu-

j ˛ acej model. Krzywa poziomicowa oznaczona kolorem niebieskim na Rys. 1.3 jest

krzyw ˛ a zamkni˛et ˛ a, z kolei dla podmodelu pozycji wa ˙zny jest jedynie pewien dys-

kretny zbiór punktów na tej krzywej.

6 Rozdział 1. Wprowadzenie

RYSUNEK1.3: Graficzna reprezentacja przykładowej struktury pod- modelu pozycji, jako krzywej poziomicowej paraboloidy eliptycznej.

Analogiczne rozumowanie mo ˙zna przyj ˛ a´c dla podmodelu orientacji. W odró ˙znie- niu od poprzedniego podmodelu, zakłada si˛e znajomo´s´c wektora orientacji stawu kolanowego θ

dla całej sekwencji. Identyfikacja polega na odnalezieniu funkcji aproksymuj ˛ acej f

, która odwzorowuje najwierniej θ

na podstawie zbioru punk- tów [ x

, y

] i składa si˛e z nast˛epuj ˛ acych kroków:

1. Wyboru struktury podmodelu orientacji f

( p

) , gdzie wektor parametrów pod- modelu oznaczono przez p

=  p

p

. . . p

, przy czym p

P IR. Ze wzgl˛edu na rozpatrywane struktury podmodelu, zapis ten jest jednoznaczny, gdy ˙z dla danego zestawu parametrów mo ˙zna opisa´c jedynie jedn ˛ a struktur˛e podmo- delu.

2. Identyfikacji optymalnych warto´sci parametrów p

przez rozwi ˛ azanie zestawu równa ´n postaci

f

( p

, x

, y

)
_θ

= 0, (1.13) które b˛ed ˛ a szczegółowo opisane w Roz. 5, przykładowo równanie (5.3). Za- le ˙zno´s´c (1.13) okre´sla macierz S równa ´n, które mo ˙zna przedstawi´c w formie uwikłanej oraz jawnej. Optymalne warto´sci parametrów zapewniaj ˛ a spełnie- nie warunku

p

= arg min

p_θ

J ( p

) , J ( p

) 

¸

} f

( p

, x

, y

)
_θ

_}

. (1.14)

3. Dla przyj˛etej struktury podmodelu oraz parametrów uzyskanych w poprzed- nim punkcie nale ˙zy wyznaczy´c ˆ θ

odpowiadaj ˛ ace warto´sciom [ x

y

] . Krok ten mo ˙zna przedstawi´c za pomoc ˛ a nast˛epuj ˛ acego schematu:

[ x

y

] _Ñ f

( p

, x

, y

) _Ñ θ ^ˆ

.

Ò Ò Ò

dyskretny opis ci ˛agły dyskretny

zbiór zbiór

(1.15)

Po zako ´nczonym procesie identyfikacji podmodel orientacji z ju ˙z dobranymi

struktur ˛ a i warto´sciami parametrów, wyznacza warto´sci orientacji stawu ko-

lanowego, na podstawie wyznaczonej zale ˙zno´sci. Podmodel orientacji jest

1.1. Motywacja i przedstawienie problemu 7

funkcj ˛ a skalarn ˛ a, która zwraca warto´s´c nale ˙z ˛ ac ˛ a do zbioru przeliczalnego tj.

f

( p

, ,  ) : IR

Ñ S spełniaj ˛ac ˛a ograniczenie f

( p

, x

, y

) = ˆθ

.

Zaproponowany model kinematyczny stawu kolanowego jest par ˛ a dwóch pod- modeli, pozycji i orientacji, st ˛ ad mo ˙zna zapisa´c

M = f

( _ ) , f

( _ )
_. (1.16) Po etapie identyfikacji nast˛epuje walidacja modelu M , która odbywa si˛e dla ju˙z wyznaczonych parametrów podmodeli. W niniejszej pracy poj˛ecie walidacja b˛edzie oznaczało weryfikacj˛e utworzonego modelu pod wzgl˛edem warto´sci bł˛edów wyj-

´sciowych oraz pod wzgl˛edem zło ˙zono´sci modelu. Zakłada si˛e, ˙ze podmodele b˛ed ˛ a posiadały indywidualny zestaw parametrów dla ka ˙zdego pacjenta. Model wyko- rzystuje wyznaczony podmodel pozycji do okre´slenia estymowanej pozycji stawu kolanowego. Nast˛epnie na bazie tej estymowanej pozycji wyznacza estymowan ˛ a orientacj˛e, przy u ˙zyciu podmodelu orientacji. Mo ˙zna zatem zapisa´c, ˙ze

x

y



Ñ f

( _ ) = 0 Ñ  ˆx

y ˆ



Ñ f

( _ ) _{Ñ ˆθ}

.

(1.17)

Zało˙zenie

W przedstawionym rozumowaniu przyj˛eto zało ˙zenie, ˙ze x

oraz y

odwzoro- wuje wzajemnie jednoznacznie θ

, wi˛ec k ˛ at ugi˛ecia kolana mo ˙zna okre´sli´c jako funkcj˛e pozycji x

y

, tj. θ

= f

( p

, x

, y

) .

Uwzgl˛edniono tu kilka istotnych aspektów dotycz ˛ acych stawu kolanowego. Po pierwsze, gdy w ruchu zauwa ˙zalna jest wył ˛ acznie rotacja, chwilowa o´s obrotu (ICR, z ang. Instantaneous Center of Rotation) jest ulokowana w obszarze ograniczonym ko ´ncem dalszym ko´sci udowej, jak pokazano na Rys. 1.2. Z kolei ¯ d jest powi ˛ azane z konfiguracj ˛ a ko´sci piszczelowej. Zatem, ICR nigdy nie pokryje si˛e z pocz ˛ atkiem układu współrz˛ednych okre´slaj ˛ acego konfiguracj˛e piszczeli.

Po drugie, podczas ruchu zginania/prostowania staw kolanowy podlega kombi- nacji rotacji i translacji, przy czym nie s ˛ a zauwa ˙zalne ruchy zło ˙zone z samego prze- suni˛ecia [66]. Zatem nast˛epuj ˛ ace relacje s ˛ a prawdziwe

x

= x

( θ

) , y

= y

( θ

) , (1.18) tzn. θ

jest jedynym argumentem funkcji opisuj ˛ acych pozycj˛e.

Wynikowe modele kinematyczne stawu kolanowego, w celach porównawczych, b˛ed ˛ a przyjmowały postaci dwóch współrz˛ednych danych równaniem (1.18). Aby zilustrowa´c zaproponowany i zdekomponowany na dwie współrz˛edne zapis, dla uproszczenia rozwa ˙zmy staw zawiasowy, zło ˙zony ze zwykłego zł ˛ acza obrotowego.

Mo ˙zna wówczas zapisa´c nast˛epuj ˛ ac ˛ a posta´c modelu kinematycznego stawu x

= _{|| ¯} d
¯g|| cos θ

,

y

= _{|| ¯} d
¯g|| sin θ

, (1.19)

gdzie przez ||  || oznaczono norm˛e euklidesow ˛a.

8 Rozdział 1. Wprowadzenie

1.2 Opis istniej ˛ acych rozwi ˛ aza ´n

Prace naukowe znane autorce wskazuj ˛ a, i ˙z badania nad opracowaniem modeli kine- matycznych stawu kolanowego dzieci prowadzone s ˛ a od stosunkowo niedawna [5].

Wi˛ekszo´s´c prac dotycz ˛ acych analizy kinematyki stawów dzieci, zazwyczaj ograni- cza si˛e do zaledwie kilku badanych pacjentów[2], jest dostosowana do opisu w ˛ askiej grupy schorze ´n lub jest oparta na teoretycznych zało ˙zeniach, a nie na rzeczywistych pomiarach. Zdecydowana wi˛ekszo´s´c opracowa ´n dotyczy danych zebranych dla do- rosłych osób.

W przypadku analizy stawów kolanowych dorosłych, mo ˙zna wymieni´c wiele rozwi ˛ aza ´n. Cz˛esto identyfikacja kinematycznych modeli stawu kolanowego opiera si˛e na analizie przebiegu ´scie ˙zki punktu kontaktu pomi˛edzy powierzchniami sta- wowymi [4,

uniemo ˙zliwia zastosowanie tej grupy modeli do analizy ruchu w stawie kolanowym małych pacjentów.

Warto nadmieni´c, ˙ze zastosowanie okre´slonej metody identyfikacji modelu jest silnie uzale ˙znione od typu uzyskanych danych pomiarowych. Mo ˙zna wymieni´c wiele technik akwizycji danych medycznych, zastosowanych do analizy kinema- tycznych modeli stawów, np. fluoroskopi˛e jedno- i dwupłaszczyznow ˛ a, RTG, ste- reofotogrametri˛e rentgenowsk ˛ a, systemy analizy chodu (np. Optitrack, Vicon), ob- razowanie metod ˛ a rezonansu magnetycznego (MRI, z ang. Magnetic Resonanse Ima- ging). W kolejno´sci zostan ˛ a przedstawione przykładowe istniej ˛ ace rozwi ˛ azania z podziałem na typ danych wej´sciowych.

Pierwsze badania nad kinematyk ˛ a stawu kolanowego wykorzystywały zdj˛ecia RTG [95]. Prace te głównie opierały si˛e na r˛ecznym oznaczaniu znaczników na zdj˛eciach ko´sci. W kolejno´sci zaproponowano u ˙zycie zdj˛e´c fluoroskopowych [69], które równie ˙z wymagaj ˛ a promieniowania rentgenowskiego, jednak dawka jest na du ˙zo ni ˙zszym poziomie

. Niew ˛ atpliw ˛ a zalet ˛ a metody jest fakt, ˙ze badania te cz˛esto wykonuje si˛e podczas operacji wspieranych obrazowaniem medycznym, wówczas promieniowanie jest nieuniknione. Wykonane zdj˛ecia s ˛ a dwuwymiarowe, w prze- ciwie ´nstwie do wszystkich kolejno wymienionych technik.

Aby zniwelowa´c główn ˛ a wad˛e zdj˛e´c fluoroskopowych, zaproponowano wyko- rzystanie fluoroskopii dwupłaszczyznowej [68,

acej dwóch aparatów fluoroskopowych umiejscowionych po obu stronach badanego. Metoda po zasto- sowaniu algorytmu nakładania zdj˛e´c umo ˙zliwia analiz˛e w trzech wymiarach. Nie- stety, wymaga zastosowania specjalistycznego sprz˛etu, rzadko dost˛epnego w pol- skich szpitalach, bada ´n nie mo ˙zna przeprowadza´c w trakcie operacji, a dawka pro- mieniowania jest wi˛eksza ni ˙z w przypadku standardowej fluoroskopii.

Stereofotogrametria rentgenowska jest równie ˙z u ˙zywana do analizy kinematyki stawu kolanowego [56]. Metoda ta cz˛esto wymaga zamocowania znaczników bez- po´srednio do ko´sci, zatem nie mo ˙ze by´c postrzegana jako mało inwazyjna i nie jest wskazana dla dzieci.

W ostatnim czasie du ˙z ˛ a popularno´sci ˛ a cieszy si˛e zastosowanie zdj˛e´c MRI do ana- lizy ruchu stawów [67]. Metoda ta jest nieinwazyjna i nie powoduje promieniowa- nia. Natomiast jej zastosowanie nie jest odpowiednie dla wszystkich aplikacji ze wzgl˛edu na ograniczon ˛ a przestrze ´n, w której wykonanie pełnego zakresu ruchów

3W niniejszej pracy do analizy wykorzystano zdj˛ecia pozyskane metod ˛a fluoroskopow ˛a, ze wzgl˛edu na mniejsz ˛a dawk˛e promieniowania w porównaniu do standardowych zdj˛e´c RTG. Poj˛ecia obrazy RTG oraz obrazy fluoroskopowe b˛ed ˛a stosowane zamiennie w odniesieniu do wej´sciowych obra- zów wykorzystanych w pracy.

1.3. Teza, cel i zakres pracy 9

jest niemo ˙zliwe. Dodatkowo, analiza kinematyki ruchu opiera si˛e na opisie krzy- wizn ko´sci, co jak wspomniano wcze´sniej, nie jest mo ˙zliwe dla osób, u których pro- ces kostnienia nie jest uko ´nczony.

Inn ˛ a interesuj ˛ ac ˛ a metod ˛ a akwizycji medycznych danych pomiarowych jest sko- rzystanie z systemu analizy chodu, wykorzystuj ˛ acego zestaw kamer oraz znacz- ników montowanych do skóry [17]. Jak pokazuj ˛ a wcze´sniejsze prace, pomiary s ˛ a obarczone znacznym bł˛edem ze wzgl˛edu na ruchy skóry wzgl˛edem ko´sci [64], co ogranicza praktyczne zastosowanie tej metody.

W opisywanej pracy analizowane b˛ed ˛ a sekwencje obrazów fluoroskopowych stawów kolanowych dzieci w wieku od 5 do 19 lat. Wszystkie kolana s ˛ a patolo- giczne, jednak nie ma nało ˙zonych ˙zadnych ogranicze ´n na typ schorzenia. Opisy- wane prace wpasowuj ˛ a si˛e w nurt prac z ostatnich lat, w których proponuje si˛e mo- dele stawów dla coraz szerszej grupy wiekowej oraz dla coraz szerszego zakresu schorze ´n.

1.3 Teza, cel i zakres pracy

Zgodnie z przedstawionymi rozwa ˙zaniami okre´slono nast˛epuj ˛ ace zało ˙zenia wst˛epne niniejszej pracy:

1. Ze wzgl˛edu na charakter danych pomiarowych, modele obejmuj ˛ a jedynie ki- nematyk˛e planarn ˛ a

2. Mo ˙zliwe jest wyznaczenie zbioru punktów kluczowych obrazów fluoroskopo- wych jako wyniku estymacji splotowej sieci neuronowej (CNN, z ang. Convo- lutional Neural Network)

3. Istnieje przekształcenie pozwalaj ˛ ace na wyznaczenie konfiguracji ko´sci pisz- czelowej i udowej na podstawie obrazów fluoroskopowych, wobec tego mo ˙z- liwe jest wyznaczenie funkcji φ

oraz φ

4. Ko´s´c piszczelowa nie przemieszcza si˛e wzgl˛edem ko´sci udowej bez zauwa ˙zal- nej rotacji, pomi˛edzy dwoma nast˛epuj ˛ acymi po sobie obrazami z sekwencji, zatem model mo ˙ze by´c zdefiniowany w funkcji k ˛ ata ugi˛ecia stawu kolano- wego, tj. x

= x

( θ

) oraz y

= y

( θ

)

5. Mo ˙zliwe jest wyznaczenie kinematycznego modelu, który wiernie odzwier- ciedla wła´sciwo´sci stawu, przy uwzgl˛ednieniu ró ˙znic anatomicznych, które s ˛ a zauwa ˙zalne dla pacjentów w ró ˙znym wieku i przy ró ˙znych schorzeniach Postawiona na podstawie powy ˙zszych zało ˙ze ´n teza brzmi:

Teza niniejszej pracy

Wła´sciwe wst˛epne przetwarzanie danych pomiarowych pochodz ˛ acych z se- kwencji zdj˛e´c fluoroskopowych, z wykorzystaniem wiedzy eksperckiej wspartej działaniem zoptymalizowanej splotowej sieci neuronowej, pozwala na zauto- matyzowan ˛ a identyfikacj˛e uniwersalnej struktury modelu stawu kolanowego uwzgl˛edniaj ˛ acego ró ˙znice anatomiczne pacjentów.

Planuje si˛e, ˙ze rezultatem prac opisywanych w niniejszej rozprawie b˛edzie al-

gorytm, który na podstawie sekwencji zdj˛e´c fluoroskopowych wyznaczy parametry