Rekbepaling in de linker hartkamer op basis van MORIo-tagged beelden. Mascha Mamhout WFW-Rapport no. 94.141
Stage verslag, Oktober 1994. Begeleider : Dr Ir P.H.M. Bovendeerd. Vakgroep fundamentele Werktuigkunde. Faculteit Werktuigbouwkunde. Technische Universiteit Eindhoven.
INHOUDSOPGAVE
1
Inleiding
2
Rekbepaling m.b.v. markers 2.1
Rekbepaling uit markerverplaatsingen
2.2
Implementatie en verificatie van methode 2.2.1
Keuze van het testprobleem
2.2.2
Uitwerking van het testprobleem
2.2.3
Vergelijking van de experimenteel en analytisch bepaalde rekken
2.2.4
pag. 13
Conclusies
pag. 20
Rekbepaling in het hart
pag. 21
3.1
Bewerking M.R.1.-tagged beelden
pag. 21
3.2
Toepassing van methode op het hart
pag. 25
Discussie
pag. 32
4.1
Discussiepunten methode rekbepaling
pag. 32
4.2
Discussiepunten rekbepaling in het hart
pag. 35
Conclusies en aanbevelingen
pag. 38
5.1
Conclusies
pag. 38
5.2
Aanbevelingen
pag. 39
Referenties
pag. 40
Bijlagen
pag. 41
3
4
5
1 INLEIDING.
Het hart zal tijdens het vemllem vam zijn functie, namelijk het rondpompen van het bloed door het lichaam, deformeren. Door middel van de mechanische modelvorming van de linkerventrikel van het hart, waarmee Bovendeerd (1990) zich tijdens zijn promotieonderzoek heeft bezig gehouden, is getracht inzicht te krijgen in de factoren die van invloed zijn op het mechanische gedrag van de linkerventrikel tot op het niveau van de individuele spiercel tijdens de ejectie- of uitdrijvingsfase. De vraag is of elke spiercel een gelijke bijdrage levert aan de pomparbeid van het gehele hart. In de ejectiefase stuwen de ventrikels een gedeelte van het bloed in de arteriële ostia, de uitstroomopeningen. In het geval van de linkerventrikel betreft dit de aorta. De totale duur van de ejectiefase is ongeveer 300 ms. De verworven kennis betreffende het mechanische gedrag van de linkerventrikel kan vervolgens worden gebruikt voor herkenning van veranderingen m.b.t de arbeidsvverdeiing die optreden bij bijvoorbeeld ischemie. Een uitgebreide bespreking van het model van Bovendeerd zal hier achterwege gelaten worden, hiervoor wordt verwezen naar referentie [6] en [8]. Een korte bespreking van de resultaten van het onderzoek is echter wel van belang. Uit het model is gebleken dat de vezeloriëntatie van de spiervezels die zich in de wand van de linkerventrikelbevinden, invloed heeft op het mechanische gedrag van de ventrikel. Door middel van al, de hoek tussen vezelrichting en de omtreksrichtingvan de ventrikel, worden de verschillende vezeloriëntaties beschreven. Door een groot aantal onderzoekers zijn experimenten uitgevoerd om te achterhalen hoe de vezelhoek precies verloopt van het endocard, de binnenkant van de ventrikelspier,naar het epicard, de buitenkant van de spier. Bovendeerd heeft twee verschillende vezeloriëntaties doorgerekend die beide binnen de spreiding van de experimenten liggen. In figuur 1.1 wordt voor beide configuraties aangegeven welke verschillen in actieve spanning in de spiercellen gevonden worden. In model 2 blijken de spanningen over de dikwandige afgeknotte ellipsoïde, waarmee de linkerventrikel in goede benadering kan worden gemodelleerd, overal hetzelfde te zijn. Dit laatste lijkt ook het meest waarschijnlijk, omdat dan alle spiervezel op dezelfde wijze belast worden gedurende een hartcyclus. Om nu te verifiëren of de spanningen over de dwarsdoorsnede van de linkerventrikel inderdaad overal ongeveer hetzelfde zijn, zou men deze spanningen willen meten.
+90 al
.......................
["I
...................
-90
bu
Fig. 1.1
-
bi
bu
bi
bu
bi
bu
bi
bu
Het effect van vezeloriëntatie op de spanningen in de spiercellen. -(Bovendeerd, 1990)
Rechtstreekse meting van deze spanningen m.b.v. bijvoorbeeld rekstrookjes is echter onmogelijk wegens de onherroepelijke schade die het hartspierweefsel zou oplopen. Een alternatief is om uit markerverplaatsingen lokale rekken te berekenen volgens de schattingsprocedure zoals beschreven door Peters (1987). Indien experimenteel bepaalde en m.b.v. het model berekende rekken met elkaar overeen komen, wordt aangenomen dat de m.b.v. het model berekende spanningen een betrouwbaar beeld geven van de spanningen die in werkelijkheid in de hartwand heersen. Voor de rekberekening m.b.v. de schattingsprocedure volgens Peters zijn twee verschillende toestanden van een deformerend lichaam op opeenvolgende tijdstippen t,, en t, nodig. Hierbij moeten markers op dat lichaam worden aangebracht waaruit verplaatsingen van tijdstip $, naar tijdstip t, kunnen worden bepaald.
e
Het aanbrengen van de markers in een doorsnede van een kloppend hart van een gezond persoon is gedaan door middel van de beeldvormende techniek M.R.I. (MagneticResonance Imaging). Met deze techniek is het mogelijk om m.b.v. het zgn. tagging een regelmatig rooster van markerpunten in een plak weefsel aan te brengen. Deze markerpunten blijven enige tijd aanwezig, waardoor m.b.v. M.R.I. opeenvolgende beelden van de doorsnede van de linkerventrikel tijdens ejectiefase kunnen worden gemaakt, waarin het markerrooster berkenbaar aznwezig is. Op de techrdek clie wordt toegepast bij deze M.W.1.-tagged beelden wordt verder niet ingegaan, hiervoor wordt verwezen naar referentie [l].
In hoofdstuk 2 paragraaf 2.1 zal de schattingsprocdure volgens Peters besproken worden. Vervolgens zal in hoofdstuk 2 paragraaf 2.2 een testprobleem worden besproken wat we hebben gebruikt om het programma REKKEN, waarin de algoritmenvoor de rekberekening zijn geïmplementeerd,te verifiëren. In hoofdstuk 3 paragraaf 3.1 zal beschreven worden hoe de M.R.1.-tagged beelden, zoals ze verkregen zijn door het A.Z.M. (Academisch Ziekenhuis Maastricht), bewerkt worden om tot een dz?a-file te kcmen met de voor onze schattingsprocedure volgens Peters bruikbare markercoördinaten. Hierna wordt het programma RHART, de aangepaste versie van het programma REKKEN, toegepast op de data-file van de linkerventrikel tijdens ejectie. Deze toepassing en de daaruit voortvloeiende resultaten zullen besproken worden in Hoofdstuk 3 paragraaf 3.2. Vervolgens worden in Hoofdstuk 4 verschillende aspecten van de toepassing van de programma’s op het testprobleem en de linkervenrtikel kritisch besproken. Tenslotte volgen in Hoofdstuk 5 de conclusies en aanbevelingen.
3
2 REKBEPALING MJ3.V. MARKERS.
Bij de theoretische modelvorming van vervorming van een lichaam wordt uitgegaan van twee naburige punten, waarvan de verschiivectoren a0, in de referentietoestand (to) en&. in de momentane toestand (ti) infinitesimaal zijn, i.e. naderen tot O. De deformatietensorlr wordt dan als volgt gedefinieerd:
Deze deformatietensor vormt weer de basis voor de definitie van de verschillende rektensoren. Wanneer we nu, zoals hier, op basis van meetgegevens deformatiegrootheden willen schatten wordt dezelfde weg gevolgd. Echter het verschil is dat er nu sprake is van eindige verschilvectoren. Mede hierdoor en door het feit dat meetwaarden altijd behept zijn met meetfouten, zijn de deformatiegrootheden niet exact te bepalen. Hoogstens is een zo goed mogelijke schatting van de gewenste grootheden mogelijk. We zullen hier dan ook de schattingsprocedure hanteren, zoals beschreven door Peters (1987). Hiertoe worden twee naburige markers P en Q op eindige afstand van elkaar beschouwd. De vectoren die P en Q in de referentie-en momentane toestand verbindenworden geschreven als respectievelijk AX(, en AZl, i.e. eindige verschilvectoren ( zie figuur 2.1).
Fig. 2.1 Lichaam M met twee markerpunten P en Q. ~
De modelfout f’ kan dan geschreven worden als:
1 i
i
Deze vergelijking vormt de basis voor de schatting van deformatiegrootheden uit de gemeten markercoördinaten. We beschouwen nu een groep van n markers, die dicht bij elkaar liggen, een z.g.n. rekgroep. Wanneer nu de coördinaten in twee toestanden (referentie en momentaan) bekend zijn, kunnen de vectoren AZ& en AZli berekend worden ( zie figuur 2.2).
Fig. 2.2 De deformatietensor wordt geschat uit de verschilvectoren van de referentie- en momentane toestand. Dit zijn de vectoren die de centrale marker P met de rest van de markers verbinden en welke als stochastischegroothedenworden beschouwd. Dit vanwege de meetfouten waarmee ze behept zijn. De volgende notatie wordt dan gebruikt:
6 = Axli- F- AZoi
; i=l,...,n
Als de posities van de markers in een rekgroep willekeurig gekozen worden, wordt de daarvan afhankelijke vector ook als stochastische grootheid beschouwd. Deze wordt gedefinieerd als:
x
. n
Met f’, de gemiddelde vector van de vectoren
8 en si de willekeurige meetfout volgt dan: -
De som van de willekeurige meetfout wordt gedefinieerd als zi :
sien het willekeurige deel van de modelfoutö&
5
Er kunnen nu schatters voor de modelfout 7 en voor de deformatietensor F S worden bepaald m.b.v. de kleinste kwadraten methode, i.e. door de volgende functie J te minimaliseren :
Deze sch2tters volgen ui: de eiseii :
Op de uitwerking zal verder niet worden ingegaan, hiervoor wordt verwezen naar referentie [4].Hier wordt alleen het resultaat gegeven :
waarbij Ago en AZl de gemiddelde verschilvectoren zijn en & en hl de tensoren die de distributie van de verschilvectoren vastleggen:
2&=
n
AXoi AXoi- A$ A% i=l
6
Vervolgens is het mogelijk om schattersvoor de verschillende rektensoren te bepalen. Voor de Green-Lagrange rektensor geldt bijvoorbeeld:
De hierboven genoemde algoritmen voor het schatten van deformatiegrootheden, zijn nu geïmplementeerd in het programma REKKEN, dat is verkregen door het aanvullen en aanpassen van het door van Ratingen verkregen programma voor rekberekening. Hiermee is het mogelijk lokale Green-Lagrange rektensoren te bepalen van een willekeurig deformerend lichaam. A l s invoer voor het programma moet gebruik worden gemaakt van data-files, waarin de coördinaten van een eindig aantal punten van dat lichaam zijn opgeslagen, en wel in twee of eventueel meerdere toestanden. Het programma REKMEN is zo geprogrammeerd dat de rekgroepenvan het deformerendelichaam uit een op te geven aantal punten n bestaan. De verschilvectoren AZ& en AZli worden vervolgens berekend t.o.v. 'net rriiddeipunt, i.e. zwaartepunt, van de veïschillenûe punten vzfi de ïekgïsey 'u? respectievelijk referentie- en momentane toestand. Dit is een belangrijk verschil met het oorspronkelijke programma, dat de verschilvectoren berekende t.o.v. één bepaalde marker maar alvorens een rekgroep omheen was bepaald. Hier wordt dus aangenomen dat het middelpunt in referentie- en in momentane toestand hetzelfde materiële middelpunt van de rekgroep is.Deze aanname is gerechtvaardigd wanneer we te maken hebben met homogene deformaties. De listing van het programma REKKEN is opgenomen in bijlage Al.
2.2 Implementatie en verificatie van methode.
Om te verifiëren in hoeverre de met de programmatuur bepaalde rekken overeen komen met de werkelijk heersende rekiren in een deformerenci ‘lichaam, Kan uitgegaan worden van een bekende situatie. Als deformerend lichaam is hier gekozen voor een incompressibele schijf met binnenstraal RI = 2 cm en buitenstraal R, = 5 cm ( zie figuur 2.3).
.-
Fig. 2.3 Geometrie van de incompressibele schijf.
-
Deze keus is gemaakt wegens analogie met de geometrie van de linker hartkamer, weliswaar sterk vereenvoudigd. Door nu aan deze schijf een bekende deformatie op te leggen, zijn de experimenteel bepaalde rekken ( met het programma REKKEN) te vergelijken met de analytisch berekende werkelijke rekken. Alvorens in te gaan op de uitwerking van de opgelegde deformaties, dient te worden opgemerkt dat wegens de gekozen geometrie van het deformerend lichaam, het eenvoudiger is met cylindercoördinaten te werken, dan met cartesische coördinaten. Er zal dan ook gesproken worden over radiele rek ( En),tangentiele rek ( E**) en afschuifrek ( E*). Dit vergt enige toevoegingen in de programmatuur.
8
In eerste instantie wordt de lokale Green-Lagrange rektensor bepaald t.o.v. een cartesische basis { ëx Zy } : EM
=
Em êxêxx+ En êyêy+Ev ëxZy+EP ZyZx
Waarbij 0 gedefinieerd is als arccos(êr-ëx) (zie figuur 2.3). In matrixrepresentatie wordt Em*%- geschreven als:
Door nu de z.g.n. rotatietensor R toe te passen op EM., volgt de Green-Lagrange rektensor t.o.v. een cylindrische basis {ë' ë+}, EqL. De rotatietensor R wordt hier gedefinieerd als:
In matrix representatie:
E=
cos4
z sin0
-sin4 cos0
1
De Green-hgrange rektensor Evt. volgt uit de matrix representatie:
Bovenstaande algoritmen zijn reeds toegevoegd aan het programma REKKEN, waarvan de listing zoals eerder vermeld is opgenomen in bijlage Al. Het programma REKKEN bepaalt dus naast E, En en E,(=EF), ook E,, E,, en E,(=E,,), respectievelijk de radiële-, tangentiële- en afschuifkek.
9
2.2.2 Uitwerking van het testprobleem.
Uitgaande van de incompressibele schijf, zal het testprobleem worden uitgewerkt, waarbij we uit- gaan van de referentietoestand zsds getoond in figuur 2.4.
Fig. 2.4 Referentietoestand van de schijf met een willekeurig punt P. De binnenstraal van de schijf wordt veranderd met een verlengingsfactor Al. Daarnaast is ook een moment aanwezig dat de buitenrand roteert over een hoek a(R2-RJ. Voor de straal van een punt op de schijf in momentane toestand geldt nu :
Verder geldt voor de hoek:
9 = 90
+
a('0-q
=
90 A 4 +
Wanneer we om beperken tot een twee-dimensionale situatie, geldt voor de positievector van een punt op de schijf in referentie toestand :
In momentane toestand wordt deze positievedor : X=rër
Voor ël in momentane toestand geldt wegens de hoekverdraaiing : ê40 ;
ër = cos(A4) êlo + &(A+)
A+= a(ro- R,)
Door substitutie van ël in de positievector volgt dan : -3
=
r . cos(A4) ëlo+P. sin(A4) ë'o
Uit de definitie van de deformatie tensor :
P= ( vo3 ) " is nu op eenvoudige wijze af te leiden dat in deze situatie voor F geldt : acos(A4) ) i?,Oët0 + (1cos(A4)) ê+oë+o r0 ar0
Op grond van de incompressibiliteitseis : det ( F )
=
ro dro
r dr
1
moet gelden : =
Door het oplossen van deze differentiaal vergelijking en substitutie hierin van :
r volgt dat voor X(r$ moet gelden :
=
A@& ro
Voor r geldt dus :
r
=
-/
Uit de definitie van de Green-Lagrange rektensor :
’
E = - (Fe- F - I ) 2
volgt door substitutie van de deformatie tensor F en de hoekverdraaiing A+ dat voor E moet gelden :
1 ar2 1 r 2 (-) Z4Zr0 + - ((-) - 1 ) Z40z+0
+-
2
ro
‘0
Door substitutievan r in deze vergelijking voor de Green-Lagrange rektensor vinden we de vergelijkingen waaraan Em, en E,,+ moeten voldoen :
E4+
=
i
(A2- 1)
1 E,,+ = - ( a 2
*
A2
*
ro )
met:
2.2.3 Vergelijking van de analytisch- en experimenteel bepaalde rekken.
De werkelijk heersende rekken afhankelijk van de straal (ro) in referentie toestand zijn nu bekend, afgezien van een te kiezen constante A, en a!. Deze rekken moeten nu vergeleken worden met de experimenteel bepaalde rekken m.b.v. het programma REKKEN. Hiervoor zijn data-files nodig met x- en y-coördinaten in referentie en momentane toestand van de defûlmerende schijf. Deze datii-files zijn gsgenereerb m.b.v. een speciaal daarvoor geschreven programma, n.1. het programma SCHWF. De listing van dit programma is opgenomen in bijlage B. Aan de hand van de te kiezen waarden voor Ai en a! zijn twee verschillende situaties geanalyseerd. In beide situaties is uitgegaan van een verlengingsfactor A, met de waarde 0,8. In de eerste situatie is voor a! de waarde O genomen, in deze situatie is dus geen hoekverdraaiing aanwezig. In de tweede situatie, waar wel sprake is van hoekverdraaiing, is voor a! de waarde 0,02 genomen. Deze waarden zijn willekeurig gekozen. Tevens zijn voor beide situaties data-files gegenereerd met twee verschillende markerdichtheden. De markerdichtheid waarbij de kleinste afstand tussen twee markerpunten 1cm is, wordt hier gedefinieerd ä l m ë i n e markeraicnthëi-d.Ee-n-gro te-marker d i c h t h e i c h i c r o r d t g e d e f i ~ ~ ~ d ~ ~die dichtheid waarbij de kleinste afstand tussen twee markers 0,5 cm is (zie figuur 2.5).
Schijf met grote markerdichtheid
Schijfmet kleiie markerdichtheid
n 0 0 0 0 0 0 0
oooooûûûûûûu.
0 0 0 0 0 0
Fig. 2.5 Schijf in referentietoestand met de verschillende markerdichtheden.
-
Door de verschillende data-files in het programma REKKEN in te laden worden nu voor bijbehorende situaties de rekken experimenteel bepaald. Hierbij wordt uitgegaan van vijf markers plus het zwaartepunt per rekgroep (n=5). Bij deze waarde van n bleken de experimentele rekken de analytische het best te benaderen. Dit is in te zien doordat bij de keuze van n=5 bij een bepaalde marker vier dichtstbijzijnde andere markers gezocht worden, die wegens het regelmatige markerpatroon netjes om deze marker heen liggen. De lokale rekken in het zwaartepunt, dat samenvalt met de middelste marker, worden nu berekend uit de verplaatsingen t.o.v. dit zwaartepunt in alle richtingen. In het programma REKKEN zijn plot-routines geïmplementeerd. Hiermee kan o.a. geïllustreerd worden welke verplaatsingen de verschillend markers van referentie naar momentane toestand ondergaan. Om een indruk te geven zijn in figuur 2.6 en 2.7 deze verplaatsingen weergegeven. De eerste situatie is zonder hoekverdraaiing en de tweede met, beide met een grote markerdichtheid.
14
Fig. 2.6 Verplaatsingen van markers op de schijf zorider hoekverdraaiing. (* geeft markers a m iïì ïefeieiitietoesta7d,-e In iz~mcntaaetoest22d)
Fig. 2.7
Verplaatsingen van markers op de schijf met hoekverdraaiing. (* geeft markers in referentie toestand aan, o in momentane toestand)
15
Het is ook mogelijk de vergelijking tussen experimenteel en analytisch bepaalde rekken te illustreren als functie van de straal r van het zwaartepunt. Voor de schijf zonder en met hoekverdraaiing (a! = O en a! = 0.02) zijn de experimenteel en analytisch bepaalde rekken (radiële-, tangentiële- en afschuifrek) te zien in respectievelijk de figuren 2.8 en 2.9. In beide gevallen is gebruik gemaakt van een grote markerdichtheid. dgerneen blijkt uit deze i!!ustratie5 dat de expeïhenteel bepaalde rekken de werkelijk heersende analytische rekken met een absolute fout kleiner dan 0,004 % benaderen. Voor de radiële- en tangentiële rek komt dit neer op een relatieve fout van respectievelijk kleiner dan 7 % en kleiner dan 4 %. Een uitzondering hierop vormen de rekken berekend in punten kleiner dan één markerafstandverwijderd van de randen. De algemene absolute fout die in deze gebieden gemaakt wordt is kleiner dan 0,07 %, wat dus ruim tien keer zo groot is als de fout buiten deze randgebieden.
De experimenteel en analytisch bepaalde rekken van de schijf zonder hoekverdraaiing (a! = O) en kleine markerdichtheid zijn te vinden in figuur 2.10. Uit deze illustraties blijkt dat de fout waarmee de werkelijke rekken bij deze kleine markerdichtheid benaderd worden kleiner is dan 60 %. Deze fout is dus aanzienlijk groter dan bij het gebruik van de grote markerdichtheid.
Radieie. rek -0.02 O
L
2
2.5
3
3.5
4
s m a l [ml
straal ícml
Fig. 2.8 Radiële- tangentiëleen afschuifrek als functie van de straal, met grote markerdichtheid, zonder afschuiving.
3-
1%
4.5
5
I 5.5
Radiele rek
0.3
0.2
t: w 0.15
1 -0.02
t
-0.06
-
-0.0s
-
-
3 w
-
I
-0.1-
*,"' I&* A.
-0.12 0.1 -
0.05
I .
:o
-0.16-
-
-0.1s
o'
r'
-0.14-
*
2
25
3
3.5
4
4.5
5
-
5.5
,J': ;
2
2.5
3
straal [cml
3.5
4
4.5
5
55
straal [cml
* afschuifrek
Fig. 2.9 Radiële- tangentiëleen afschuifrek als functie van de straal, met grote markerdichtheid, met afschuiving.
0.1
0.09
-
0.0s
-
0.07
-
0.060.05
-
0.04 -
0.03 I . 2
25
3
3.5
4
straal [cml
45
5
5
Radielerek
tangentiele rek
-0.12
:*
,' =:
-0.14-
-0.16-
- :
-0.18
1 2.5
2
3
straat icml
3.5
4
stmal [cml
afschuifrek
xi05
I
Fig 2.10 Radiële- tangentiëleen afschuifrek als functie van de straal, met kleine markerdichtheid, zonder afschuiving.
* *
*
s
*
4.5
5
5.5
2.2.4 Conclusies.
In het algemeen is op de binnenstraal van de schijf een grotere spreiding om de werkelijke rekwaarden te zien. Bij de afschuifrek is dit ook het geval op de buiten straal, hoewel hier rekening moet worden gehouden met de kleinere schaalverdeling van de y-as. Deze spreiding is toe te schrijven aan de randeffecten. Aan de randen van de schijf kan n.1. geen rekgroep gevoïmd woïden waarin de markers symmetrisch om het zwaartepunt gerangschikt zijn, wat het gevolg is van het gebruik van eindige verschilvectoren. Specifiek valt nog op dat bij de schijf met hoekverdraaiing er ook een duidelijk grotere spreiding van de werkelijke radiële- en tangentiële rekwaarden te constateren is op de buitenstraal. Dit is weer toe te schrijven aan de randeffecten die in dit geval groter zijn wegens grotere verplaatsingen aan de buitenstraal veroorzaakt door de hoekverdraaiing. Wegens al eerder genoemde meetfouten, gebruik van eindige verschilvectoren en de aanname dat het middelpunt in referentie toestand hetzelfde materiële punt is als het middelpunt in momentane toestand, zullen we altijd rekenjng moeten bEjvm hoiider, =et afijkkgen ra;; de werkelijke waarden. Deze fouten zullen kleiner worden naarmate er gewerkt wordt met grotere markerdichtheden en dus kleinere verschilvectoren.
20
3 REKBEPALING IN HET HART.
Vanaf ongeveer het begin van de ejectie fase van de hartcyclus zijn om de f 20 ms M.R.1.-beelden gemaakt van de doorsnede van de linker ventrikel. Voor deze beelden is gebruik gemaakt van gezonde proefpersonen. In totaal zijn de eerste tien beelden relevant bevonden voor verdere bewerking. Deze tien beelden beslaan dus een tijdsbestek van ongeveer 200 ms. We kunnen er dan ook vrijwel zeker van zijn dat deze beelden ook werkelijk de linker ventrikel tijdens de ejectiefase weergeven. Dit wegens het feit dat de ejectiefase ongeveer 300 ms duurt. In figuur 3.1 zijn enkele van de M.R.1.-tagged beelden geïllustreerd zoals ze verkregen zijn van het Academisch Ziekenhuis Maastricht.
~
Fig. 3.1 M.R.1.-tagged beelden van een transversale doorsnede van een thorax met het hart van een gezond proefpersoon. (frame 1,5 en 10) Om te komen tot voor ons bruikbare markercoördinaten, moeten de beelden eerst bewerkt worden. De informatie over deze beeldverwerking is verkregen via F.Aelen, BíoQsica, R.E.
21
De conventionele M.R.1-tagged beelden, waarin de markers aanwezig zijn, bevatten veel randinformatie die niet nodig is. Daarom wordt ingezoomd op dat gedeelte van het beeld waar de linker hartkamer zich bevindt. Dit ingezoomde beeld wordt gefilterd, zodat de markers beter te detecteren zijn. Dit filteren gebeurt door het beeld te transformeren naar het spatiële frequentie domein. Hierdoor ontstaat een 2-dimensionaal frequentie spectrum (zie figuur 3.2).
Fig. 3.2 2-D frequentie domein van M.R.1.-tagged beeld. -
Wanneer we nu alleen die frequenties doorlaten die essentieel zijn voor het rooster, welke te zien zijn in figuur 3.3, ontstaat een beeld waarin de markers duidelijker zichtbaar Zijn.
22
Fig. 3.3 Gefilterd frequentie domein van M.R.1.-tagged beeld.
-
Na het filteren moeten de middelpunten van de markers worden bepaald, om zo de posities van de markers te achterhalen. Dit gebeurt m.b.v een drempelwaarde. Dit is een bepaalde in te stellen grijstint. Pixels met een grijstint boven deze drempelwaarde krijgen de waarde 1(wit) en pixels met een grijsiint beneden deze drempelwaarde krijgen de waarde O (zwart). Hierdoor ontstaat een binair beeld van de markerverdeling, i.e. wel of geen marker. In figuur 3.4 is een serie van beelden te zien waarin het uiterst linkse beeld een origineel M.R.1.-beeld is, het middelste beeld een gefilterd beeld, en het rechtse beeld een 'gedrempeld', binair beeld.
23
Fig. 2.4
Serie met orgineel (links), gefilterd (midden) en gedrempeld (rechts) beeld.
De dddeipunteri vam de markers hmefi 1iu bepadcl wordeii door de gemiddelde Ä= e=.ywaarde per marker te bepalen. Het aantal pixels waaruit een marker mag bestaan wordt hiertoe eerst gelimiteerd (bijv. minimaal 5 pixels en maximaal 30 pixels). Van elk frame worden op deze manier de markerposities bepaald. Tenslotte moet nog bepaald worden welke markers van de verschillende frames bij elkaar horen. Dit gebeurt met een speciaal daarvoor ontwikkeld programma van AMuytjens en J. Roos, Medische Informatica en Statistiek, R.L. Het uiteindelijke resultaat is een data-file met markercoördinaten in tien verschillende frames, i.e. toestanden, vanaf ongeveer het begin van de ejectie fase.
24
3.2 Toepassing van methode op het hart
Aangenomen dat het programma REKKEN aan de hand van markerverplaatsingen een goede benadering geeft van de werkelijk heersende rekken, zal het programma nu worden toegepast op het hart. Hiertoe is het programma REKKEN aangepast, waaruit het -programma RHART is ontstaan. De listing van dit programma is te vinden in bijlage A2. De essentie van het programma, n.1. de rekberekening uit verplaatsingen is uiteraard hetzelfde gebleven. De veranderingen betreffen vooral de verschillende plotroutines voor het presenteren van gegevens en resultaten. Van de data-file van het hart, zoals die van F. Aelen verkregen is, Zijn in figuur 3.5 de markerbanen van frame 1 t/m frame 10 te zien.
*
Fig. 3.5 Markerbanen van oorspronkelijke data-file. (* geeft markers in referentietoestand weer) Wegens verwachtingen betreffende de geometrie van de doorsnede van de linkerventrikel en het feit dat de gradiënten in het verplaatsingsveld niet te groot mogen worden, is uit figuur 3.5 geconcludeerd dat enkele van de ge'ïllustreerde markers geen deel uitmaken van de linkerventrikelwand. Vermoedelijk djn deze markers verbonden met omringend weefsel zoals bijvoorbeeld het borstbeen. De oorspronkelijke data-file is dan ook aangepast door deze punten uit te sluiten.
25
De overgebleven markers met weer hun banen tijdens de tien doorlopen toestanden, zijn te zien in figuur 3.6. In deze figuur is ook de middelpuntsverplaatsing, berekend uit de gemiddelde x- en y-coördinaten van alle markers per frame, ge-ïllustreerd.
y
f
$ 4
f
-
I?//- * 'y
L
0 1
t k E \ \
\
fi
L
\-*4,
b Y
h
k
ss\~L %L.\
Fig. 3.6 Markerbanen en middelpuntsverplaatsing van de aangepaste data-file. (* geeft markers in referentietoestand weer, + geeft middelpunt in ïefireiitietoeJtmd weer) Van de baan van dit middelpunt hebben we aangenomen dat deze de starre verplaatsing van de linker ventrikel weergeeft. Door deze middelpuntsverplaatsing van alle individuele verplaatsingen af te trekken, houden we de markerbanen over zoals te zien in figuur 3.7. Het in deze figuw weergegeven plus-teken geeft het vastgezette middelpunt weer ten opzichte waarvan de markers verplaatsen.
Fig. 3.7 Markerbanen min de starre translatie van de linker ventrikel. (* geeft markers in referentietoestand weer, + geeft vastgezette middelpunt weer)
26
Vervolgens hebben we de berekende gemiddelde straal van de linker ventrikel uitgezet tegen de framenummers, i.e. de achtereenvolgens doorlopen toestanden, (zie figuur 3.8).
gemiddelde straal van het hart per frame i
n 4 v)
o>
.ï
a
U 4
(d
e
Y
v!
38l
38
O
2
4
6
8
10
frame nummer
Fig. 3.8 Gemiddelde straal van de linkerventrikel t.o.v. de framenummers. Uit deze figuur blijkt dat de straal nagenoeg lineair daalt als functie van de tijd, i.e. framenummer. Vanwege deze continue daling van de gemiddelde straal kan worden aangenomen dat de tijdsperiode vanaf frame 1t/m 10 inderdaad binnen de tijdsperiode van de ejectiefase ligt. We kunnen nu m.b.v. het programma RHART de rekken experimenteel bepalen, waarbij we gebruik maken van de aangepaste data-file. Als referentietoestand kiezen we frame 1en als momentane toestand nemen we frame 10. Verder gaan we weer uit van vijf markers per rekgroep, omdat dit bij de toepassing van het programma op het testprobleem de beste keus bleek.
27
In figuur 3.9 zijn de hoofdrekken geïllustreerd zoals ze berekend zijn door het programma RHART. Hierbij geeft de getrokken lijn de hoofdrek weer in de richting waar de grootste verlenging optreedt en de gestreepte lijn de hoofdrek in de richting van de grootste verkorting.
Fig. 3.9 Rekveld met de hoofdrekken in de linkerventrikelwand tijdens ejectie.
.
Uit deze figuur blijkt dat de heersend hoofdrekken enigszins in radiële en tangentiële richting zijn, wat ook logisch is wegens het afnemen van de gemiddelde straal van de linkerventrikel tijdens de ejectiefase. Bij het testprobleem hadden we te maken met een volledig cylindersymmetrische situatie. We gaan er nu van uit dat de geometrie van de doorsnede van de linker ventrikelwand ook cylindersymmetrisch is en zullen de rekken dan ook weer presenteren als radiële- tangentiële- en afschuifrek als functie van de straal van het zwaartepunt van de verschillende rekgroepen, gemeten vanaf het berekende middelpunt van de linker ventrikel.
28
Wanneer figuur 3.5 en figuur 1.1met elkaar vergeleken worden, kan vastgesteld worden dat het gedeelte van de ventrikelwand links in figuur 3.5 overeenkomt met het septum, de scheiding tussen linker- en rechterventrikel. Het gedeelte rechts in deze figuur representeert de vrij wand van de linkewentrkel. Om een globale indruk te krijgen van de te verwachten onderlinge verschillen in de rekken in deze verschillende gedeelten van de linkerventrikel is een kwadrantindeling toegepast zoals geïllustreerd in figuur 3.10.
\
Y
I
/*
/
* * * * * *\
Fig. 3.10 Kwadrantindeling van de linkerventrikel met de markers in referentietoestand. De radiële- tangentiële- en afschuifrek zijn te zien in figuur 3.11. De verschillende tekens in deze figuur geven aan in welk kwadrant een bepaalde marker, waarvan de rek wordt geillustreerd ligt (* geeft het eerste kwadrant aan, + het tweede kwadrant, o kwadrant drie en x kwadrant vier).
29
tangentiede rek per k w ~ ( ~ l , + = 2 , 0 = 3 , x 4 )
*
O*
t:
0.4
-
8
O
PI
o.
0.2
-
O
oxxx
-0.05
-
-0.1
-
0-
-
*
O
30
+
-0.2
-
-0.25
-
0
+
++ox
Lc
*o +++ o
-
0
+
+ xxx
x
* X x x X X
-0.3
O
40
O
+
-0.15 -
O
-0.2
+ O
+
+ +=x + + c %x” +t ++ xd(+
o
50
60
straal ipixekl
30
40
50
straal [pix&]
Fig. 3.11 Radiële- tangentiëleen afschuifrek als functie van de straal van de linkerventrikel (De tekens geven aan in welk kwadrant de rekken berekend werden)
straal [pixels1
60
Uit deze figuur is over het verband van de verschillende rekken met de straal weinig te zeggen. Hooguit is uit de puntenwolk een gemiddelde waarde te schatten die een zeer algemene indruk geeft van de heersende radiële-, tangentiële- en afschuifrek in de ventrikelwand. Wanneer we de rekken per kwadrant bekijken valt vooral bij de radiële rek op dat de berekende rekken vooral in het tweede en vierde kwadrant een relatief geringe spreiding vertonen. Als we nog eens naar figuur 3.10 kijken blijkt dit ook te verklaren uit het feit dat er in deze kwsrdrmten weral twee markers =i; de doorsiiede va= de wand liggen, waardoor de markers per rekgroep meer symmetrisch om het zwaartepunt heen liggen en de berekende rekken dus meer met elkaar overeen komen. Verder lijkt het niet reëel om uit de illustraties van de rekken iets zinnigs te zeggen over het verband van de rek met de straal die loopt van binnen- naar buitenstraal, wegens het geringe aantal markers op de doorsnede van de ventrikelwand bij een bepaalde hoek @ (gemiddeld namelijk twee). De zwaartepunten ten opzichte waarvan de rekken berekend worden liggen dus ongeveer op dezelfde straal. Concluderend is er dus alleen iets zinnigs te zeggen over de grote van de rekken in het algemeen en eventueel per kwadrant. Vrijwel met zekerheid -hmm we stellen dat we te maken hebben met een positieve radiële rek en een negatieve tangentiëleen afschuifrek.
i
,
4 DISCUSSIE.
4.1 Discussiepunten methode rekbepaling.
Bij het testprobleem zoals beschreven in hoofdstuk 2 zijn een aantal keuzes gemaakt betreffende de ~ j z van e de reitberekening. In hoeverre deze Keuzes correct zijn vait te bezien. Het gaat hier om de aanname betreffende het zwaartepunt van een rekgroep en de keuze van vijf markers per rekgroep. De keuze om de lokale rekken te berekenen t.o.v. het zwaartepunt van een rekgroep is gemaakt omdat de experimentele rekken op de randen van de schijf bij deze keuze minder afwijkingen van de analytische rekken vertoonden. Dit is in te zien door te kijken wat er gebeurt bij de rekgroepsbepalingvan een marker die precies op de buitenrand ligt. Wanneer de ïsk bepaald zou xorden t.o.+'. deze marker z m d m de mmksrs imamit de verplaatsingsvectorenbepaald worden allen op een Meinere afstand van het middelpunt van de schijf liggen dan deze marker, waardoor de experimentele rek iets hoger is dan de werkelijk heersende rek. Door nu alvorens de rek uit de verplaatsingsvectoren te berekenen, het zwaartepunt c.q. middelpunt van een rekgroep te bepalen en t.o.v. dit punt de rekken te berekenen zullen de experimentele rekken op de randen de analytische rekken beter benaderen, doordat de berekende gemiddelde rek van een rekgroep nu gekoppeld is aan de gemiddelde positie in een rekgroep, n.l. die van het zwaartepunt. Daarnaast is het zo dat het gebruik van rekgroepen met een berekend zwaartepunt alleen dan te verantwoorden is wanneer er sprake is van homogene deformaties, dus als elk infinitesimaal klein blokje materie dezelfde deformatie ondergaat. Het materiële blokje dat in referentie toestand zwaartepunt is moet dat in momentane toestand dus ook nog zijn. Bij het testprobleem dat wij hebben gebruikt is dit echter niet helemaal het geval. Het blijkt dat als we twee punten nemen, één op de binnen- en één op de buitenrand van de schijf bij enkel radiële deformatie ( A, =0,s en a! = O ) de berekende positie van het zwaartepunt in momentane toestand van deze twee punten 0,06 cm afwijkt van de positie van het werkelijke materiële punt dat in referentie toestand verbonden was met het zwaartepunt. We hebben hier dus te maken met een fout van ongeveer 2 %. Wanneer er ook nog afschuiving aanwezig is zal deze fout nog iets hoger liggen.
We moeten echter niet vergeten dat deze fout van f 2 % voortkomt uit een situatie waarbij een zeer kleine markerdichtheid gehanteerd wordfn.1. de kleinste afstand tussen twee markerpunten is 3 cm.Bij de markerdichtheden die wij hebben gebruikt zal deze fout dus iets kleiner zijn dan 2 %. Bij de keuze van wel of geen toepassing vapp een maartepunt in een rekgroep is dus de afweging gemaakt tussen of het verlagen van de relatief grote randeffecten en een kleine algemene fout bij de rekberekening of het elimineren van deze algemene fout met als consequentie de grotere randeffecten. Zoals al eerder vermeld is voor het eerste gekozen, n.1. het toepassen van een zwaartepunt in een rekgroep, omdat in de gebruikte geometrie van de schijf en het hart, er relatief veel randen aan deze lichamen zitten, en de markerdichtheid klein is. Zoals in Hoofdstuk 2 al vermeld werd is gekozen voor vijf markers per rekgroep. Wanneer namelijk werd uitgegaan van een grote markerdichtheid zoals gedefinieerd in Hoofdstuk 2 (kleinste afstand tussen twee markers 0,s cm), bleken bij deze keuze de experimenteel bepaalde rekken de analytische werkelijke rekicen het best te benaderen. Dit is te veriiiaren doordat bij deze keuze in de meeste rekgroepen de verdeling van markers om het zwaartepunt symmetrisch is. Tevens valt het zwaartepunt samen met de middelste marker, zie figuur 4.1.
I
I
*
Fig. 4.1 Rekgroep zods die gemiddeld gekozen za7 worden. (* geeft markers weer, o geeft het berekende zwaartepunt weer)
33
Door de symmetrische verdeling van de markers om het zwaartepunt worden de rekken in dit punt het best benaderd, omdat deze nu berekend worden uit verschilvectoren die symmetrisch om het zwaartepunt heen liggen. Bij een keuze van negen markers in een rekgroep (n=9) liggen de markers sok 5ymnetrisch om het zwaartepunt heen. Uit de rekberekening met het programma REKKEN bleek echter dat bij de keuze van n=9 de experimentelerekken aan de randen aanzienlijk minder goed de werkelijk heersende rekken benaderden. Dit is te begrijpen wanneer men bedenkt dat bij de keuze van n=9 het inhomogene rekveld over een groter gebied gemiddeld wordt.Tevens komt het zwaartepunt van een grotere rekgroep aan bijvoorbeeld de buitenrand van de schijf, dichter bij het middelpunt van de schijf komt te liggen. Hierdoor wordt het verschil tussen de buitenstraal en de grootste straal van de schijf waar nog lokale rekken worden berekend groter. Men zou zich kunnen afvragen of de genoemde effecten niet onderdrukt kunnen worden door de keuze van een veel grotere markerdichtheid. Door het werken met grotere markerdichtheden nemen de eindige verscbihxtûïen îìaïììelijl~in gïmtte af, waardoor de deformatietensor beter geschat zal worden en er dus ook een betere benadering van de lokale Green-Lagrange rektensor bepaald kan worden. De reden waarom in het testprobleem toch voor een relatief lage markerdichtheid is gekozen ligt in het feit dat m.b.v. M.R.1.-tagged beelden er maar een beperkte markerdichtheid gerealiseerd kan worden.
34
4 2 Discussiepunten rekbepaling in het hart.
Bij de toepassing van het programma RHART op het hart kunnen een aantal punten betreffende de betrouwbaarheid van de resultaten ter discussie worden gebracht. Over de M.R.1.-tagged beelden dient te worden opgemerkt dat deze verkregen zijn uit een plakje weefsel van 5 8 mm dik. Dat betekent dus dat de bewegingen van ae markers gemiddeld worden over deze plak, wat fouten met zich meebrengt. Daarnaast is het zo dat de positie van de plak vastligt. Bij beweging van het hart in longitudinale richting in de mens wil dat zeggen dat elk beeld in een serie een ander plakje weefsel weergeeft. Er is daarentegen aangenomen dat in de hier gebruikte serie van tien frames, elk frame hetzelfde plakje weefsel weergeeft. Ook dit leidt dus tot fouten. Over de grootte orde van deze fouten is echter niets bekend. r T.
uit Îiguur 3.í in hoofdstuk 3 hebben g e zzíîgeíì~mefidat edele narterp~~teri piei verbonden zijn met het weefsel van de linkerventrikel. Vervolgens hebben we deze markers uitgesloten bij de rekberekening. In hoeverre wij de juiste markers hebben meegenomen in de rekberekening (de markers die dus wel met het linkerventrikel weefsel waren verbonden) is niet met zekerheid te zeggen. Om dit wel met zekerheid te kunnen zeggen, zou men al tijdens de beeldbewerking van de conventionele M.R.1.-beelden die markers moeten uitsluiten die niet verbonden zijn met het linkerventrikel weefsel. In de conventionele M.R.1.-beelden kunnen de verschillende weefsels met de daarop aangebrachte markers namelijk nog goed herkend worden.
Wanneer we figuur 3.4 nog eens goed bekijken valt op dat rechtsonderaan in deze figuur twee markers in de linkerventrikelholte lijken te liggen. Het is aannemelijk om te veronderstellen dat deze markers verbonden zijn met een papillairspier. Deze twee afwijkende markers zorgen voor fouten in de rekberekening, temeer wegens de zeer kleine markerdichtheid die is toegepast, waar we later nog terug komen. Ook deze markers zouden tijdens de M.R.1.-beeldbewerking al herkend moeten zijn en moeten zijn uitgesloten van de data-file.
35
In hoofdstuk 3 hebben we aangenomen dat de geometrie van de dwarsdoorsnede van de linkerventrikel cylindersymmetrischis. Het is nu natuurlijk de vraag in hoeverre dit correct is. Bekijken we hiertoe figuur 4.2, welke de dwarsdoorsnede van het hart ter hoogte van de equator weergeeft tijdens systole en diastole, bG&t toch dat de Enkerventrikel, rechts gelegen in de figuur, geometrisch bij goede benadering gezien kan worden als een cylindrische schijf. Of de linkerventrikel qua materiaal, bijvoorbeeld de vezeloriëntatie, cykdersymmetrisch is vait te betwijfeien.
Fig. 4.2 Anatomische doorsnede van hat hart ten hoogte van de equator. (A. tijdens systole. B. tijdens diastole)
36
Tenslotte is het de vraag in hoeverre de bepaalde rekken de werkelijke rekken benaderen. Uit Hoofdstuk 2 bleek dat wanneer we een markerdichtheid gebruikten waarbij zeven markers op de doorsnede van binnenstraal neer buitenstraal waren aangebracht, de experimentele rekken met een acceptabele spreiding OIII de werkelijke waarden bepaald werden. Als er op deze doorsnede maar vier markers waren gebruikt, bleek de spreiding al aanzienlijk groter te worden. Misschien kunnen we hier zelfs al beter over afwijkingen spreken dan over spreiding. -Wanneer we nu de situatie in net geval van de linkerventrikel bekijken, zien we dat hier gemiddeld slechts twee markers op de doorsnede van binnennaar buitenstraal aantreffen. Het is dus terecht om de verkregen resultaten met enig voorbehoud te bezien.
5 CONCLUSIES EN AANBEVELINGEN.
5.1 Conclusies.
Naar aanleiding van de resultaten bij de toepassing van het programma REKKEN op het tesrpr~bleemPE het progra~maREMRT op de Eïìkervenîrillei ‘kunnen we de volgende conclusies trekken:
1. Voor een betrouwbare rek bepaling is het vereist om een zo groot mogelijke markerdichtheid te gebruiken, waardoor fouten met als oorzaak het gebruik van eindige verschilvectoren kleiner worden.
2. Tevens is gebleken dat wanneer we de markers in een rekgroep symmetrisch om het punt waar men de ïokaïe rek wil weten legt, deze lokale rek het best de werkelijk heersende rek benadert. 3. In de doorsnede van de linkerventrikelwand wordt de vormverandering gedurende ejectie gekarakteriseerd door een positieve radiële rek en een negatieve tangentiëleen afschuifrek.
Opgemerkt dient nog te worden dat de keuze van een zo groot mogelijke markerdichtheid ook van belang is wanneer we iets willen zeggen over het verloop van de rekken over de doorsnede van de linkerventrikelwand. De reden hiervoor is logischerwijs dat we zoveel mogelijk punten willen hebben over deze doorsnede bij een bepaalde hoek @, waarna van elk van die punten de lokale rek bepaald kan worden.
38
5.2 Aanbevelingen.
Als aanbeveling voor vervolgstudie kan derhalve in overweging genomen worden dat m.b.v. de M.R.1.-tagged beelden de heersende rekken goed benaderd kainnen worden mb.v het programma RHART. Als voorwaarde hiervoor geldt dan wel dat er een markerdichtheid gebruikt moet worden die ongeveer vier keer zo groot is als de markerdichtheid die toegepast is in de door ons gebruikte M.R.I.-tagged beeiden. Een aantai markers van minimaal acht over de doorsnede van de linkerventrikel is dus een vereiste voor betrouwbare lokale rekberekening. Bij deze minimale markerdichtheid wordt dan een aantal markers per rekgroep van vijf aanbevolen (n=5). In hoeverre men een minimale markerdichtheid van acht markers over de linkerventrikelwand m.b.v. de M.R.1.-tagging techniek kan verwezenlijken is onbekend, hiervoor zou het Academisch Ziekenhuis Maastricht geïnformeerd moeten worden.
39
~
6 REFERENTIES
B.R. Friedman et al., Principles of M.R.T., Me Craw-Hill Inc.1989. J.A. Bernards en L.N. Bouman, Fysiologie van de mens, vijfde druk, Houten/Antwerpen, 1988. Barry J. Anson, Morris’ Human Anatomy, twaalfde druk, 1966. G.W.M. Peters, Tools for the measurement of stress and strain fields in soft tissue, Proefschrift, Rijksuniversiteit Limburg, Maastricht, 1987. Jery L. Prince en Elliot R. Me Veigh, Motion estimation from tagged MR image sequences, IEE transaction on medical imaging, vol.11, 110.2, juni 1992. P.H.M. Bovendeerd, Dependence of local left ventricular wall mechanics on myocardial fiber orientation : a model study, J. Biomechanics, vo1.25, no. 10, pag. 1129 -1140, 1992. 1 I
T. Arts, Description of the deformation of the left ventricle by a kinematic model, J. Biomechanics, vo1.25, no.lO, pag. 1119 -1127, 1992. P.H.M. Bovendeerd, Influence of endocardial-epicardial crossover fibers on left ventricular wall mechanics, J biomechanics, vo1.27, no.7, pag. 941-951, 1994. M.R. van Ratingen, Mechanical identification of inhomogeneous solids, Proefschrift Technische Universiteit Eindhoven, mei 1994.
I
l
iI I
I
Ii
i I
I
I
I
40
7 BIJLAGEN In deze bijlage zijn de listings van de in matlab geschreven programma’s bijgevoegd. De programa’s in bijlage A zijn verkregen door het aanvullen en aanpassen van het door M.van Ratingen verkregen programma voor rekberekening. Het programma in bijlage B is geschreven voor het genereren van data-files behorende bij de op verschillende manieren deformerende schijf. Bijlage A l
programma REKKEN.
Bijlage A2
programma RHART.
Bijlage B
programma SCHIJF.
BIJLAGE A l :
REKKEN.M
hold off echo off clear clg axis ( 'square') % % % % % %
Deze M-file is bedoeld om data files te laden. Hierna worden rekken bepaald aan de hand van verplaatsingen en rekken berekend aan de hand van bekende verlengingsfactor(1abda) en eventueel bekende afschuiffactor(a1fa).
fiienaml = input(;Geef naam van file : : f ' s ' j ; filenam2 = input('Geef extensie van file : ' , 'S') ;
filename=[filenaml,'.',fiïenam2];
disp(/Laden data-file') eval(['load
',filename]);
disp(/Data-file geladen') mrks = eval( [filenaml, ( :,1)'3 ) ; nmrk = max(mrks) ; coorxx = eval([filenaml,'(l:nmrk,2)']); cooryy = eval([filenaml,'(l:nmrk,3)']); npix = eval([filenam2]); coorx = (lü.*coorxx)./npix; coory = (10.*cooryy) /npix; coordi = [ coorx coory];
.
= mean (coordi); xg = 0(1,1); Yg = 0(1,2); coorxn = coorx-xg; cooryn = coory-yg; fi = atan2 (cooryn,coorxn);
O
%
%
x en y waarde van middelpunt.
% %
x en y waarden ten opzichte van middelpunt.
.....................................
dist = zeros(nmrk,nmrk); r = [I;
mpix = 12; Vpix = [O mpix
0
mpix];
disp('Bepa1en markerafstanden') for i = 1:nmrk xl = coorx(i); y1 = coory(i) ; xro = coorxn(i) ; yro = cooryn(i) ;
% afstand matrix (dist) vullen.
sqrt (xroA2+yroA2) ; [r rol;
ro
=
r
=
j
= 1;
while j <= i distance = sqrt( (xl-coorx(j) ) ^2+(yl-coory( j) ) ^2); dist(i,j) = distance; dist (j,i) = distance; j = j+l; end; end ; m i n = min(coorx) ; xmax = max(coorx) ; ymin = min(cosrYj ymax = max(coory) ;
hulp1 = sort(dist( :) ) ; dissel = hulpl(nmrk+l);
%
Alle afstanden op volgorde van groote
axis ('square') cont=l; while cont==l Epsl Eps2 Eps3 Exx EYY EXY Veigl Veig2 frame rekgr Hoofdrekken Eg1 Err Efifi Erfi Cylrekken rE cxm cym
=
= = = = =
= = = = = =
= = = = =
= =
[I;
El; El; El;
[I; [I; [I; [I; [I; [I; [I; [I; [I; [I; El; [I; [I; [I; [I;
%
matrixen nodig voor rek berekening.
nstr = input('Aanta1 markers in rekgroep is : '1; case = input('We1k frame moet vergeleken worden met ref. frame : dicp('referentie frame wordt vergeleken met momentane frame')
deformxx deformyy
= =
eval([f~lenaml,~((case-l)*nmrk+l:case*nmrk,2)~]); eval([filenaml,~((case-l)*nmrk+l:case*nmrk,3)~])~
deformx deformy
=
(ïû.*deformxx)./npix; (lO.*deformyy)./npix;
=
axis (Vpix);
'1;
polymark(coorx,coory,'g*') polymark(deformx,deformy,'ro')
for i = 1:nmrk polyline ( [ coorx ( i) deformx (i)3 , [ coory (i) deformy ( i)3 ) end X = [O mpix mpix 0 O]; Y = [O O mpix mpix O]; polyline(X,Y) title(['Markers voor ',filename,' in ref. en moment. frame ( * and o)']) pause; eval(['meta ',filenaml,'nl']); eval(['!gpp ',filenaml,'nl ldps']); clc ClCJ
disp('rekken deform
= [
strgrp
=
for i
=
worden berekend') deformx deformy];
zeros(nmrk,nstr) ;
1:nmrk re rEE rekgr
% %
% coordinaten van momentane frame.
frune
% loopt over alle markers. = r(i); = [rEE
re]; = [rekgr i] ; = [ f r a z e case] ;
[distl, distl-i]
=
%
matrixen die nodig zijn voor
sort(dist(:,i)); % % % % %
strgrp(i, : ) = distl-i(1:nstr) ';
% afstand tussen marker n en andere markers op volrorde van groote. matrix met markers die in de rekgroep van elke observatie marker horen.
= coorx (strgrp(i,:) ,:) ; coorxrkg cooryrkg = coory(strgrp(i,:),:); deformxrkg = deformx(strgrp(i,:),:); deformyrkg = deformy (strgrp(i, :) , :) ; coorxm = mean (coorxrkg); coorym = mean (cooryrkg); cxm = [cxm coorxm] ; CYm = [cym coorym]; deformxm = mean(deformxrkg); deformym = mean(deformyrkg); rrg = sqrt( (coorxm-xg)^2+(coorym-yg)^ 2 ) ; rE = [rE rrg];
point = [ coorxm coorym]; % x en y coordinaten van mark dpoint = [ deformxm deformym]; % ref. frame en moment. frame pointar = ones(nmrk,l)*point; % 80 bij 2 matrixen met 80 keer x dpointar= ones(nmrk,l)*dpoint; % waarde van marker n van de beid R
% %
fik
=
= [
cos(fik) -sin(fik); sin(fik) cos(fik) 3 ;
atan2((coorym-yg),(coorxm-xg)); %
rotatie matrix.
berekend dx0mean dvr = (coordi(strgrp(i,:),:) dx0mean = mean (dvr); berekend dxlmean
-
pointar(strgrp(i,:),:));
dvd = (deform(strgrp(i,:),:) - dpointar(strgrp(i,:),:)); dxlmean = mean (dvd); berekend XOO and XO1 xOO= dvrr*dvr./nstr-dxOmeanr*dxOmean; xO1= dvrr*dvd./nstr-dxOmeanf*dxlmean;
%
F E Ecyl
=
Exx EYY Exy Err Efifi Erfi
=
= =
xO1r*inv(xOO); 0.5*(Fr*F-eye(2)); R/*E*R;
[Exx E(l,l)]; = [EYY E(2,2)1: = rExY E(lf211; = [Err ECyl(i,ij]; = [Efifi Ecy1(2,2)]; = [Erfi Ecyl(l,2)];
[Veig,epsilon]
=
eig(E) ;
% % %
% %
berekend deformatie tensor. berekend Green Lagrange
% % % %
rek ‘censor (E). van alle observatie markers worden de Green Lagrange rekken in rijen gezet.
de eigenvectoren (Veig) en eigenwaarden (epsilon) van E worden in een 2 bij 2 matrix gezet.
if epsilon (1) >= epsilon (2) Epsl = [Epsl epsilon(l,l)]; Eps2 = [ Eps2 epsilon (2,2)3 ; Eps3 = CEps3 0.5*(epsilon(l,l)-eps~lon(2,2))]; Veigl = [veigi Veig( :,i) 3 ; Veig2 = [Veiga ~eig(:,2)]; % de hoofdrekken worden else % van alle markers in een ri Epsl = [Epsl epsilon(2,2)]; Eps2 = [ Eps2 epsilon (1,1)3 ; Eps3 = [ Eps3 O. 5* (epsilon( 2 ,2 ) -epsilon ( 1,1)) 3 ; Veigi = [veigi veig(:,2)]; Veig2 = [Veigz veig(:,î)]; end ; next i
%
...
end;
clc
xx
=
input(’0m de rekgroepen per marker te zien, type [l] :
if xx
==
1
for i = 1:nmrk clg polymark(coorx(strgrp(i,:)),coory(strgrp(i,:)),~gx~)
hold on polymark(cxm(i) ,cym(i) ,‘wor),pause hold off end end %
epslmax = ceil(max(Eps1) *250)/250; epslmin = -ceil(max(-Eps2)*250) /250; V = [O epslmax epslmin O]; axis (V);
clg plot(Epsl,Eps2,'go') grid title(['Rek-domein voor ',filename]) xlabel('Eps -1') ;ylabel('Eps -2') pause strl
str2 fac
= =
=
Veigl'.*[Epsl' Veig2'.*[Eps2'
Epsl']; Eps2'1;
1.4*dissel/max(max(abs([strl str21)));
coordim = [cxm' cynf]; M11 = coordim(l:nmrk,:)-0.5*fac*strl; M12 = coordim(l:nmrk,:)+0.5*fac*strl; M21 = coordim(l:nmrk,:)-0.5*fac*str2; M22 = coordim(l:nmrk,:)+0.5*fac*str2; xll = Mll(:,l); x12 = M12 (:, 1);
x21 = M21(:,1); x22 = M22(:, 1); y11 = Mll(:,2); y12 = M12(:,2); y21 = 6121(:,2); y22 = M22(:,2); c w axis([xmin xmax ymin ymax]); for i = 1:nmrk if Epsl(i) > O polyline([xîl(i) else polyline( [xll(i) end if Eps2(i) > O polyline ( [x21(i) else polyline ( Ex21 (i) end end
x12(i)],
[yïi(i) yï2(i)],fwr);
xi2 (i)1 ,
[yïî(i) y12 (i)I ,
:r
)
;
x22 (i)], [y21(i) y22 (i)3 ,'w') ; x22 (i)3, [y21(i) y22 (i)3, : ' ) ; f
% % %
title(['Rek-veld voor ',filename]); xlabel([rHoofdrekken met ',num2str(nstr),' markers in rekgroep']) pause %eval(['meta f,filenaml,rn2']); %eval(['!gpp ',filenaml,'nL ldps']); clc
disp(' De analytische rekken worden nu vergeleken met de experimentele')
echo off clg
Erran = Efifian = Erfian =
[I; [I; [I;
rmnr deformxm deformym rm rmnm fim [rv,r-l] mmxr fimn fimx
= =
= = = =
= =
min(r); deformx-xg; deformy-yg; sqrt((deformxm.*czformxm)+(d f rmym. *deformym) ) ; min(rm) ; atan2(deformym,deformxm); sort (r); max(rv-1) ;
= fijrrirrixrj; =
fim(mmxr);
rmnm./rmnr; rmnr ; max (r); Ri-(Ri./lO); Ro+(Ro./lO); (fimx-fimn) / (Ro-Ri); (Ro-Ri)./loo; [Ri:int:Ro]; length (strl);
% %
labda, binnen- en buiten-straal worden bepaald.
.
z z = input(’A1s er sprake is van afschuiving type[2], zo niet type[l]:
if z z == 1; for i=i:q; str = strl(i) ; Lstr = sqrt(l+(Ri^2./str.^2)*(LA2-1)) ; Erra = .5*((l./(Lstr).^2)-1); Efifia = .5* ( (Lstr) ^2-1) ; Erfia = O; % rekken worden analytisch Erran = [Erran Erra]; % berekend Efifian = [Efifian Efifia]; Erfian = [Erfian Erfia]; end; end;
.
if z z == 2; for i=i:q; str = strl(i) ; Lstr = sqrt(l+(Ri^2./str.^2)*(LA2-l)) ; Erra = .5*( (l./(Lstr) .^2)+( (alf*Lstr*str). ^ 2 ) - ï ) ; Efifia = .5*( (Lstr).^2-1); Erfia = .5*(alf*LstrA2*str.^2)./(str); % rekken worden analytisch Erran = [Erran Erra] ; % berekend. Efifian = [Efifian Efifia]; Erfian = [Erfian Erfia]; end ; end p’ Emar Emer2 Emerl
max (Erran); max(Err) ; = min(Err) ; = =
if Emar>=Emer2, Emx2r = Emar;
else Emx2r end
=
Emer2;
Emxlr = Emerl; Emr2
Emx2r+Emx2r./lO;
=
if Emxlr>=O, Emrl = O; else Emrl = Emxlr+(Emxlr./lO); end
Emafi = min (Efifian); Emefil = min(Efifi); Emefi2 = max (Efifi) ; if Emafi>=Emefil; Emxlfi = Emefil; else Emxlfi = Emafi; end; Emfil Emx2fi
Emxlfi+(Emxlfi./lO); Emefi2;
= =
if Emx2fi>=O, Emfi2 = ErnefiL+(Emefi2./10); else Emfi2 = O; end B
=
[Rir Ror Emfil Emfi23;
Emarfil = max(Erfian); Emarfi2 = min(Erfian); Emerfi1 = max (Erfi) ; Emerfi2 = min (Erfi) ; if Emarfil>=Emerfil, Emxrfi = Emarfil; else Emxrfi = Emerfil; end if Emarfi2>=Emerfi2, Emnxrfi = Emerfi2; else Emnxrfi = Emarfi2 end Emrfi
=
Emxrfi+(Emxrfi./lO);
if Emnxrfi>=O Emnrfi = Emnxrfi-(Emnxrfi./lO); else
% % %
analytische en experimentele rekken als functie van de straal worden in een grafiek gezet.
Emnrfi = Emnxrfi+(Emnxrfi./lO); end C
=
[Rir Ror Emnrfi Emrfi];
axis(A) ; plot(rE,Err,'r*',strl,Erranfrg--t 1 title(['Radiele rekr]); xlabel([/straal [~m]~]); ylabel(['Err']); pause eval(['meta ',filenaml,'rk']); evai(['igpp ',fiienaml,'rk i,,--','1 ) ; clc clg 23-
axis (B) ; plot(rE,Efifi,/r*',strl,Efififian,'g--r 1
title(['tangentiele rek']); xlabel(['straal [cmIr]); ylabel( [rEfi€i']) ; pause eval(['meta ',filenaml,rrk']); eval([[!gpp fl,filenaml,rrkldps']); clc clg
axis ( C ) ; plot(rE,Erfi,/r*r,strl,Erfianf'g--' 1 title( [ 'afschuifrek']) ; xlabel(['straal [cm]']); ylabel(['Erfi']); pause eval(['meta eval(['!gpp
',filenaml,rrkr]); /,fi1enamltrrkldps']);
CIC cont=input('voor end; end;
bepalen rekken andere situatie, type[l]:
I);
BIJLAGE A2 :
FU3ART.M
hold off echo off clear Cu3 axis ( square' ) % % % %
Deze M-file is bedoeld om data files te laden van het hart. Hierna worden rekken bepaald aan de hand van verplaatsingen.
filenaml = input('Geef filenam2 = input('Geef
naam van file F extensie van file : f
' S f )
;
', ' S f ) ;
filename=[filenaml,f.f,filenam2];
disp('Laden
data-file')
eval(['load
',filename]);
disp ( Data-file geladen' ) = eval( [filenaml, ( : ,2) 3) ; Yc = eval([filenamî,' (:,3)']); LX = length (Xc); mrks = eval([filenaml,f(:,l)f]); nmrk = max(mrks) ; nfr = LX. /nmrk; npnt = (nfr.*nmrk)-nmrk;
xc
coorx = eval([filenaml,f(l:nmrk,2)f]); coory = eval([filenamlff(l:nmrk,3)f]); coordi = [ coorx coory]; O
=
xg Yg coorxn cooryn fi
%
mean (coordi); 0(1,1); coorx-xg; coory-yg; atan2 (cooryn,coorxn) ;
dist = zeros(nmrk nmrk) ; r = [I; = = =
eval([filenam2]); npix+2; [O mpix û mpix];
XX = zeros(nfr,l); YY = zeros(nfr,l) ; Xg = zeros(nfr,l); Yg = zeros(nfr,l);
k = 1;
x en y waarde van middelpunt.
0(1,2);
**** * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
npix mpix Vpix
%
% x en y waarden ten opzichte van % middelpunt.
while k < nfr for i = l:nmrk:npnt+l Xg(k) = mean(Xc(i:i+nmrk-1)); Yg(k) = mean(Yc(i:i+nmrk-1)); k = k+l; end end minx = min(Xg)-1; maxx = max(Xg)+l; miny = min(Yg)-l; maxy = max(Yg)+l;
axis([minx maxx miny maxy]); % baan van de starre translatie wordt geplot polymark(Xg(1) ,Yg(l) I /g+/) polyline(Xg,Yg) A = [minx maxx maxx minx minx]; B = [miny miny maxy maxy miny]; polyline(A,B) title(['baan.van de starre translatie voor /,filename]) pause %eval([lmeta J,filenami,fl/]); %eval([/!gpp r,filenaml,'l ldps']); clc clg axis (Vpix); polymark(coorx,coory,rr*') polymark(Xg(1) ,Yg(l) ,rg+') polyline(Xg,Yg) for i = 1:nmrk k = 1; while k < nfr; for j = 0:nmrk:npnt XX(k) = Xc(i+j); % markerbanen worden geplot. YY(k) = Yc(i+j); k = k+l; end end polyline(XX,YY) end X = [ O mpix mpix 0 O]; Y = [ O 0 mpix mpix O]; polyline(X,Y) title(['markerbanen voor /,filename]) pause %eval(['meta /,filenami,f2/]); %eval([/!gpp /,filenarnl,/2/dps/]); clc clg % markerbanen min starre translatie % worden geplot.
dltXg = zeros(nfr,1); dltYg = zeros (nfr,i) ; dltXg(1) = o; dltYg(1) = O; for i = 2:nfr dltXg (i) = Xg (i-1)-Xg (1); dltYg (i) = Yg (i-1)-Yg (1);
end XXc = zeros(nfr*nmrk,1); YYc = zeros(nfr*nmrk,1); for ifr = 1:nfr for jmrk = 1:nmrk XXc(jmrk+(ifr-1) *nmrk)=Xc(jmrk+(ifr-1) *nmrk)-dltXg(ifr) ; YYc(jmrk+(ifr-1) *nmrk)=Yc(jmrk+(ifr-1) *nmrk)-dltYg(ifr) ; end end axis (Vpix); =^lym^rk!CoQ~x;coory~/r*/) polymark(Xg(1) ,Yg(ij ,;g+;j for i = 1:nmrk k = 1; while k < nfr; for j = 0:nmrk:npnt XX(k) = XXc(i+j); YY(k) = YYc(i+j); k = k+l; end end polyline(XX,YY) end x = [ O mpix =pix o o ] ; Y = [ O 0 mpix mpix O]; polyline(X,Y) title([/markerbanen min de starre translatie voor ',filename]) pause %eval(['meta f,filenam1,f3/]); %eval([/!gpp /,filenaml,'3 /dps/]);
= [I; rf rfs = [I; nframe = [I;
Rcx RCY
= =
xc-xg; Yc-xg;
for j=O:nmrk:npnt for i=1: nmrk Rcx2 = (Rcx(i+j) ) * (Rcx(i+j)) ; Rcy2 = (Rcy(i+j)) * (Rcy(i+j)) ; rfr = sqrt(Rcx2+ Rcy2); rf = [rf rfr]; end srfl = rf ((j+l) : (j+nmrk)); srf = sum(srf i) ; rfgem = srf./nmrk; rfs = [rfs rfgem]; end
fX
mrfsl mrfs2 ryl
= = = =
nfr+nfr./lO; min(rfs); max(rfs) ; mrfsl-mrfsl./lO;
% gemiddelde straal van het hart wordt % bepaalt en geplot als functie % van het frame no.
ry2
=
mrfs2+mrfs2./10;
axis([O fx ryl ry2]); plot (rfs ) title('gemidde1de straal van het hart per frame') xlabel('frame nummer') ylabel('gemid.straa1 [pixels]') pause eval(['meta ',filenaml,'4']); eval(['!gpp ',filenaml,'4 /dpsf]); clc clg cont=l; while cont==l caseref
=
input('we1k
frame kies je als ref. frame :
'1;
coorx = eval([f~lenaml,'((caseref-l)*nmrk+l:caseref*nmrk,2)~]); coory = eval([f~lenaml,'((caseref-l)*nmrk+l:caseref*nmrk,3)~]); disp('Bepa1en markerafstanden')
for i
1:nmrk xl = coorx(i) ; y1 = coory(i) ; % afstand matrix (dist) vullen. xro = coorxn(i) ; yro = cooryn(i) ; ro = sqrt(xroA2+yroA2); r = [r rol; j = 1; while j c= i distance = sqrt((xl-coorx(j))A2+(yl-coory(j))A2); dist (i,j) = distance; dist (j,i) = distance; j = j+l; end; end; =
xmin = min(coorx) ; xmax = max(coorx) ; ymin = min(coory) ; ymax = max(coory) ; hulp1 = sort(dist( :) ) ; dissel = hulpl(nmrk+l);
%
Alle afstanden op volgorde van groote
axis ( 'square')
Epsl Eps2 Eps3 Exx EYY
=
= =
= =
[I; 61; [I; [I; [I;
% matrixen nodig voor rek berekening.
EXY Veigl Veig2 frame rekgr Hoofdrekken Eg1 Err Ef ifi Erfi Cylrekken rE
= =
= = = = = =
= = = = =
CWm
cym
= =
flZW
case nstr
= =
[I; [I; [I; [I; [I; [I; [I; [I; Cl; [I; [I; [I; [I;
c1; ij;
input('We1k frame moet vergeleken worden met ref. frame : input('Aanta1 markers in rekgroep is : '1;
'1;
disp('referentie frame wordt vergeleken met momentane frame')
deformx = eval([filenaml,~((case-l)*nmrk+l:case*nmrk,2)~]); deformy = e v a l ( ~ f i l e n a m l , ~ ~ ~ c a s e - l ~ * n m r k + l : c a s e * n m r k , 3 ~ ~ ~ ~ ~
axis (Vpix); polymark(coorx,coory,'go') polymark(deformx,deformy,'r*')
for i = 1:nmrk polyline ( [ coorx (i) deformx (i)3 , [ coory (i) deformy (i)3 ) end X = [O mpix mpix 0 O]; Y = [O 0 mpix mpix O]; polyline(X,Y) title(['Markers voor ',filename,' in ref. frame ',num2str(caseref),' en mom. f pause; %eval(['meta ',filenaml,'5']); %eval(['!gpp ',filenaml,'5 ldps']); clc clg disp('rekken deform
= [
worden berekend')
deformx deformy];
% coordinaten van momentane frame.
strgrp = zeros(nmrk,nstr); for i = 1:nmrk rekgr frame
% loopt over alle markers. = =
%
[rekgr i]; [frame case];
[distl, distl-i]
=
matrixen die nodig zijn voor
% afstand tussen marker n en % andere markers op volrorde van % groote. % matrix met markers die in de
sort(dist(:,i));
strgrp(i,:) = distl-i(1:nstr)';
% rekgroep van elke observatie % marker horen. = coorx (strgrp(i,:) , :) ; coorxrkg cooryrkg = coory (strgrp(i, :) , :) ; deformxrkg = deformx (strgrp(i, : ) , :) ; deformyrkg = deformy(strgrp(i,:),:); coorxm = mean (coorxrkg); coorym = mean (cooryrkg); cxm = [cxm coorxm] ; CYm = [cym coorym] ; deformxm = mean (deformxrkg) ; deformym = mean (deformyrkg) ; rrg = sqrt ( ( coorxm-xg) 2+ ( coorym-yg) 2 ) ; rE = [FE: rrgj;
point = [ coorxm coorym]; % x en y coordinaten van mark dpoint = [ deformxm deformym]; % ref. frame en moment. frame pointar = ones(nmrk,l)*point; % 80 bij 2 matrixen met 8 0 keer x dpointar= ones(nmrk,l)*dpoint; % waarde van marker n van de beid fik fizw R
% % %
=
=
atan2((coorym-yg),(coorxm-xg));
[fizw fik]; cos(fik) -sin(fik); s i a a ( f i k ) cos(fik1 1;
= [
% rotatie matrix.
berekend dx0mean dvr = (coordi(strgrp(i, :) , :) - pointar (strgrp(i,:) ,:) ) ; dx0mean = mean (dvr); berekend dxlmean dvd = (deform(strgrp(i, : ) , : ) - dpointar(strgrp(i, : ) , : ) ) ; dxlmean = mean (dvd); berekend XOO and XO1 xOO= dvrr*dvr./nstr-dxOmeanr*dxOmean; xO1= dvr'*dvd./nstr-dxOmean'*dxlrnean;
F E Ecyl
=
Exx EYY EXY Err Efifi Erfi
= = = = =
= =
xûlf*inv(xûû); 0.5*(Fr*F-eye(2)); Rf*E*R;
[Exx E(l,l)]; CEYY E(2,2)1; [EXY E(1,2)1; [Err Ecyl(l,l)]; [Efifi Ecy1(2,2)]; = [Erfi Ecyl(l,2)];
[Veig,epsilon] = eig(E) ;
% berekend deformatie tensor. % berekend Green Lagrange
% % % %
rek tensor (E). van alle observatie markers worden de Green Lagrange rekken in rijen gezet.
% de eigenvectoren (Veig) en eigen% waarden (epsilon) van E worden % in een 2 bij 2 matrix gezet.
if epsilon(1) >= epsilon(2) Epsl = [Epsl epsilon(1,l)l; Eps2 = [Eps2 epsilon (2,2)3 ; Eps3 = [Eps3 0.5*(epsilon(l,l)-epsilon(2,2))1; Veigï = [Veigï Veig(:,ï)]; Veig2 = [Veig~Veig(:,2)]; % de hoofdrekken worden else % van alle markers in een ri Epsl = [Epsl epsilon(2,2)];
Eps2 = [Eps2 epsilon(1,l) 3 ; Eps3 = [Eps3 0.5*(epsilon(2,2)-epsilon(l,l))]; Veigï = [Veigl Veig(:,2)]; Veig2 = [Veigz Veig(:,ï)]; end; next i
%
...
end; clc xx = input('0m
if
de rekgroepen per marker te zien, type [l] :
'1;
XX == 3.
for i = 1:nmrk clg polymark (coorx(strgrp(if:) ) coory (strgrp(iI :) ) 'gxr) hold on polymark(cxm(i) ,cym(i) , 'wo') ,pause hold off end end %
epslmax eps2min
= =
ceil(max(Eps1) *250) /250; -ceil(max(-Eps2)*250)/250;
V = [O epsîmax eps2min O]; axis(V); clg plot(Epsl,Ep~2,~go') grid title(['Rek-domein voor ',filename]) xlabel('Eps -1') ;ylabel('Eps -2') pause strl = Veigl'.*[Epsl' str2 = Veig2'.*[Eps2/ fac
=
Epsl']; Eps2'1;
1.4*dissel/max(max(abs([strl str21)));
coordim = [cxm' cym']; M11 = coordim(l:nmrk,:)-0.5*fac*strl; M12 = coordim(l:nmrk,:)+0.5*fac*strl; M21 = coordim(l:nmrk,:)-0.5*fac*str2; M22 = coordim(l:nmrk,:)+0.5*fac*str2; xll x12 x21 x22 y11 y12 y21 y22
Mll(:,l); M12 (: ,i) ; M21(:,1); = M22(:,1); = Mll(:r2); = M12(:,2); = M21(:,2); = M22(:,2); = = =
clg axis([xmin xmax ymin ymax]);
% hoofdrekken worden geplot.
for i = 1:nmrk if Epsl(i) > O polyline([xll(i) else polyline ( [xlï(i) end if Eps2(i) > O polyline ( [ x21 (i) else polyline ( [x21(i) end end
x12(i)],
[yiï(i) yï2(i)],/w/);
x12 (i) 3 , [y11(i) y12 (i)3 ,
I--'
1;
x22 (i)3 , [ y21 (i) y22 (i)3 I /w' ) ;
x22 (i)], [y21(i) y22 (i)],I - - '
1;
% % %
title(['Rek-veld voor /,filename,/frame r,num2str(case),~t.o.v. frame /,num2 xlabel([/Hoofdrekken met f,num2str(nstr),t markers in rekgroep']); pause eval(['meta f,filenaml,r6/]); eval(['!gpp rlfilenaml,f6/dps/]); clc clg % De bepaalde rekken worden nu geillustreerd.
rmnr Ri Ro Rir Ror
= = = = =
min(r) ; rmnr; max(r); Ri-(Ri./lO); Ro+(Ro./lO);
% labda, binnen- en buiten-straal % worden bepaald.
Emer2 = max(Err) ; Emerl = min(Err) ; Emr2
=
Emer2+Emer2./10;
if Emerl==O Emrl = -0.01; else if Emerl>O Emrl = Emerl-Emerl./lO; else Emrl = Emerl+Emerl./lO; end end A
=
[Rir Ror Emrl Emr21;
axis (A) ; plot(rE,Err, title([/Radiele rek ( * ) voor xlabel1 [ straal/]) ;
,filenameIrframe ',num2~tr(case),~t.o.v. fram
ylabel(['Err pause % eval(['meta % eval(['!gpp clc clg
met f,num2str(nstr),f markers in rekgroep']); ',filenam1,'7f]); ',filenaml,'7 ldps']);
Emefi1 = min (Efifi) ; Emefi2 = max (Efifi) ; Emfi2
=
% analytische en experimentele % rekken als functie van de straal % worden in een grafiek gezet.
Emefi2+EmefiL./iO;
if Emefil==O Emfil = -0.01; else
if Emefil>O Emfil = Emefil-Emefil./lO; else Emfil = Emefil+Emefil./lO; end end
B
=
[Rir Ror Emfil Emfi21;
cixis(B) ;
plot(rE,Efifi,'r*') title(['tangentiele rek ( * ) voor ',filename,' frame ',nurn2~tr(case),~t.o.v. ~label([~straal']); ylabel(['Efifi met ',num2~tr(nctr),~markers in rekgroep']); pause % eval(['meta f,filenaml,'8f]); % eval(['!gpp ',filenarnl,'8 /dpsf]); clc clg Emerfi2 = max (Erfi) ; Emerfi1 = min (Erfi) ; if Emerfi2>=O Emrfi2 = Emerfi2+Emerfi2./10; else Emrfi2 = Emerfi2-Emrfi2./10; end if Emerfil>=O Emrfil = Emerfil-Emerfil./lO; else Emrfil = Emerfil+Emerfil./lO; end C
=
[Rir Ror Emrfil Emrfi21;
axis (C); plot(rE,Erfi,{r*') title(['afschuifrek ( * ) voor ',filename,' frame /,num2str(case),' t.o.v. fram ~label([~straal']); ylabel(['Erfi met ',num2str(nstr),' markers in rekgroep']); pause % eval([fmeta f,filenaml,'9f]); % eval(['!gpp ',filenaml,'9 /dpsf]);
%cont=input('voor bepalen rekken andere situatie, type[l] anders Enter: %end; %end; cont2 =l; while cont2 == 1
% % % % %
iiskken worden nü per kwâdrârìt bepaalt erì a p a r t
in een grafiek geplot m.b.v. tekens.
hoek = input('Geef kwad = input('Geef rEk Errk Efifik Erfik
= = = =
de hoek t.o.v.de x-as, waar het eerste kwadrant begint: het te illustreren kwadrant ( r e k ilustratie): I ) ;
[I; [I; [I; [I;
if kwad == 1 hoekl = hoek; hoek2 = hoek+(0.5*pi);
k = 1; for i = 1:nmrk if fizw (i)>hoekl&fizw (i)
= 1;
for i = 1:nmrk if fizw (i)>hoek1 I fizw (i)
k = 1; for i = 1:nmrk if fizw (i)>hoekl&fizw (i)
hoek2 rEk (k) = rE(i); = Err(i); Errk(k) Efifik(k) = Efifi(i); Erfik(k) = Erfi(i); k = k+l; else k = kf end end end
hoekkl = hoek+(kwad-l)*0.5*pi; hoekk2 = hoekkl+0.5*pi; min(r) ; rmnr; max(r); = Ri-(Ri./lO); = Ro+(Ro./lO) ;
rmnr Ri Ro Rir Ror
= = =
Emer2 Emerl
= =
Emr2
=
% %
labda, binnen- en buiten-straal worden bepaald.
max(Errk) ; min(Errk);
Emer2+Emer2./10;
if Emerl==O Emrl = -0.01; else if Emerl>O Emrl = Emerl-Emerl./lO; else Emrl = Emerl+Emerl./lO; end end A
=
[Rir Ror Ernar1 ErnrL];
axis ( A ) ; plot(rEk,Errk, Ir*') title(['Radiele rek ( * ) voor ',filename,' frame ',num2~tr(case),~ t.o.v. fram xlabel(['straal (kwadrant van ',n~m2str(hoekkl),~ tot f,num2str(hoekk2),' rad
ylabel(['Err pause % eval(['meta % eval(['!gpp clc clg
met ',num2str(nstr),' markers in rekgroep']); ',filenaml,'kl']); ',filenaml,'kl /dps']);
Emefi1 Emefi2
= =
min (Efifik) ; max (Efifik) ;
Emfi2
=
Emefi2+~mefi2./iO;
analytische en experimentele rekken als functie van de straal % worden in een grafiek gezet. % %
if Emefil==O Emfil = -0.01; else if Emefil>O Emfil = Emefil-Emefil./lO; else Emfil = Emefil+Emefil./lO; end end B
=
[Rir Ror Emfil Emfi21;
axic(B) ; plot(rEk,Efifik,'r*') title(['tangentiele rek ( * ) voor ',filename,' frame ',num2str(case),' t.o.v. xlabel(['straal (kwadrant van ',num2str(hoekkl),' tot ',num2str(hoekk2),' rad ylabel(['Efifi met ',num2str(nstr),' markers in rekgroep']); pause % eval(['meta ',filenaml,'k2']); % eval(['!gpp ',filenaml,'kL /dps']); clc clg
Emerfi2 = max (Erfik) ; Emerfif = min (Erfik) ; if Emerfi2>=0 Emrfi2 = Ernerfi2+Emerfi2./10; else Emrfi2 = EmerfiL-Emrfi2./10; end if Emerf il>=O Emrfil = Emerfil-Emerfil./lO; else Emrfil = Emerfil+Emerfil./lO; end C
=
[Rir Ror Emrfil Emrfi21;
axis ( c ) ; plot (rEk,Erfik,'r* ) title(['afschuifrek ( * ) voor ',filename,' frame ',num2str(case),' t.o.v. fram xlabel([lstraal (kwadrant van ',num2str(hoekkl),~tot ',num2str(hoekk2) ,' rad ylabel(['Erfi met ',nurn2str(nstr),' markers in rekgroep']); pause % eval(['meta ',filenaml,'k3']); % eval(['!gpp ',filenaml,'k3 ldps']);
~ont2=input(~voor ilustreren rek ander kwadrant type[l], anders Enter: end % % % % % % %
');
De rekken per kwadrant worden tenslotte allen geplot in een grafiek, waarin de tekens aangeven in welk kwadrant een bepaalde marker waarvan de rek bepaald werd ligt.
disp('Er pause; rEkl Errkl Efifikl Erfikl rEk2 Errk2 Efifik2 Erfik2 rEk3 Errk3 Efifik3 Erfik3 rEk4 Errk4 Efifik4 Erfik4
volgt nu een gecombineerde illustratie van de rekken per kwadrant.') =
= = =
= =
= =
= = = =
= = =
=
for hoekl hoek2
[I; [I; [I; [I; [I; [I; [I; [I; [I; [I; [I; [I; [I; [I; [I; [I;
= =
hoek; hoek+(0.5*pi);
k = 1; for i = 1:nmrk if fizw (i)>hoekl&fizw (i)hoek1 I fizw (i)
Efifik2(k) = Efifi(i); Erfik2(k) = Erfi(i); k = k+l; else k = k; end end end for hoekl = hoek-pi; hoek2 = hoek- (O. 5*pi) ;
k = 1; for i = 1:nmrk if fizw (i)>hoekl&fizw (i)hoek2 = rE(i); rEk4 (k) Errk4(k) = Err(i); Efifik4(k) = Efifi(i); Erfik4 (k) = Erfi (i); k = k+l; else k = k; end end end
Errk Efifik Erfik Errkmn Efifikmn Erfikmn
%hoekkl %hoekk2 rmnr Ri Ro Rir Ror Emer2 Emerl
= = = = = =
= =
[max(Errkl),max(Errk2),max(Errk3),max(Errk4)]; [max(Ef~f~kl),max(Ef~f~k2),max(Ef~f~k3),max(Ef~f~k4)]; [max(Erf~kl),max(Erf~k2),max(Erfik3),max(Erf~k4)]; [min(Errkl),min(Errk2),min(Errk3),min(Errk4)]; [min(Efifikl),min(Efifik2),min(Efifik3),min(Efifik4)]; [min(Erfikl),min(Erfik2),min(Erfik3),min(Erfik4)];
hoek+(kwad-l)*0.5*pi; hoekkl+0.5*pi;
min(r) ; rmnr; = max(r) ; = Ri-(Ri./lO); = Ro+(Ro./lO); =
=
= =
max(Errk) ; min (Errkmn);
% %
labda, binnen- en buiten-straal worden bepaald.
Emr2
=
Emer2+Emer2./lO;
if Emerl==O Emrl = -0.01; else if Emerl>O Emrl = Emerl-Emerl./lO; else Emrl = Emerl+Emerl./lO; end end
A
=
[Rir Ror Emrl Emr21;
axis(A) ; plot(rEkl,Errkl,/r*/,rEk2,Errk2,",rEk3/Errk3//bo',rEk4/Errk4//wx/) title(['Radiele rek per kwadrant(*=l,+=2,0=3,~=4)~]);
xlabel(['straal [pixels] )/I); ylabel(['Err']); pause eval([/meta ',filenarnl,'kl']); eval(['!gpp ',filenaml,'kl /dps/]); clc clg Emefil = min(Efifikmn); Emefi2 = max (Efifik) ;
% analytische en experimentele % rekken als functie van de straal
% worden in een grafiek gezet.
Emfi2
=
Emefi2+Emefi2./10;
if Emefil==O Emfil = -0.01; else if Emefil>O Emfil = Emefil-Emefil./lO; else Emfil = Emefil+Emefil./lO; end end B
=
[Rir Ror Emfil Emfi21;
axis (B); plot(rEkl,Efifikl,'r*',rEk2,Efifik2,'g+',rEk3/Efifik3/'bo//rEk4/Efifik4//~~/) title(['tangentiele rek per kwadrant(*=l,+=2,0=3,~=4)']);
xlabel(['straal [pixels]']); ylabel(['Efifi']); pause eval(['meta ',filenaml,'k2/]); eval(['!gpp ',filenamlrfk2/dps/]); clc clg Emerfi2 = max (Erfik) ; Emerfi1 = min (Erfikmn) ; if Emerfi2>=O Emrfi2 = Emerfi2+Emerfi2./10; else Emrfi2 = Emerfi2-Emrfi2./10;
end if Emerfil>=O Emrfil = Emerfil-Emerfil./lO; else Emrfil = Emerfil+Emerfil./lO; end C
= [Rir Ror Emrfil Emrfi21;
axis (C); plot!rEkl,Efifikl,'r*',rEk2'Efifik2.'g+"rEk3~Efifik3,'bo',rEk4,Efifik4,'~~') title([iafscnuifrek per kwaàrant(*=i,t=2,0=3,x=4j'j); xlabel i[ 'straal [pixels] 3 ) ; ylabel(['Erfi']); pause eval(['meta ',filenaml,'k3']); eval(['!gpp ',filenaml,{k3 ldps']); clc clg cont=input('voor end end
bepalen rekken andere situatie, type[l] anders Enter:
');
cont=input('voor end end
bepalen rekken andere situatie, type[l] anders Enter:
');
~~
BIJLAGE B :
SCHIJF.M
hold off echo off clear clg clc axis ('square') % % % % % % %
Deze M-file is bedoeld om data-files te genereren van een deformerende schijf. Die deformatie wordt bepaald door een op te geven verlengingsfactor (labda) en rotatiehoek (alfa).
npix
=
input('geef het aantal pixels:
x = o; y = o; xref = [I; yref = [I; mrk = [I; xsr = [I; Ysr = [I; coordiref =
[I;
disp ('rasterpunten worden bepaald') for j
Y
0:npix j+l; for i = 0:npix x = i+l; y = j+l; xref = [xref x]; yref = [yref y]; end =
=
end xgm Ygm Ri Ro mr
mean (xref); mean(yref) ; = npix./5; = npix./2; = (npix+l)^2; =
=
for k = 1:mr xs = xref (k); ys = yref (k); xsch = xgm-xs; ysch = ygm-ys; rs = sqrt (xschA2+ysch^2) ; if rs<=Ro & rs>=Ri xsr = [xsr xs]; ysr = [ysr ys]; end end nmr = length(xsr) ; mrk = (l:nmr]; coordiref = [mrk' xsr' ysr'];
% Hoe groter het aantal pixels % hoe groter de markerdichtheid.
I);
disp('coordinaten
markers in de referentie toestand zijn bepaald')
xref = coordiref(l:nmr,2); yref = coordiref(l:nmr,3); coorref = [xref yref]; L alfa
= =
input ( 'geef labda (Ri): input('geef alfa: '1;
xmom ymom npnt
=
[I; [I; [I;
=
xgem ygem
= = mean (coorref); = Qs(l!l); = ûs(I,Lj;
for i
=
Os
xr yr ro Lro
'1;
1:nmr; xref(i); = yref(i); = sqrt( (xr-xgem)"2+(yr-ygem) "2) ; =
= (sqrt(l+(RiA2./roA2)*(L"2-l))); = Lro*ro; fi = atan2 ( (yr-ygem), (xr-xgem)) ; fia = fi+alfa* (ro-Ri); xm = xgem+r* (cos(fia)) ; ym = ygem+r* (sin(fia) ) ; xmom = [xmom;xm]; ymom = [ymom,ym];
r
npnt
=
[npnt'i];
end % de coordinaten in momentane toestand % zijn bepaald.
coordiref = [npnt' xref yref]; coordimom = [npnt' xmom' ymom']; schffi = [coordiref;coordimom]; save schijfd.out schffi /asdi /double disp('momentane coordinaten(schijfd.out) zijn bepaald') % % %
De data-file is gegenereerd.
I
