DE60201923T2

DE60201923T2 - Materialverbiegungssimulationsgerät inbesondere für weiche Körpergewebe

Info

Publication number: DE60201923T2
Application number: DE60201923T
Authority: DE
Inventors: Nicholas Ayache; Delingette; Guillaume Picinbono
Original assignee: Institut National de Recherche en Informatique et en Automatique INRIA
Current assignee: Institut National de Recherche en Informatique et en Automatique INRIA
Priority date: 2001-02-08
Filing date: 2002-02-04
Publication date: 2005-11-24
Anticipated expiration: 2022-02-05
Also published as: EP1231568B1; DE60201923D1; US20020183992A1; FR2820534B1; ATE282862T1; FR2820534A1; US7239992B2; EP1231568A1

Description

Die Erfindung betrifft die Verformungssimulation von Materialien bzw. Stoffen, insbesondere von körperlichen Weichteilgeweben.
Üblicherweise ist eine Visualisierungsschnittstelle vorgesehen, um wenigstens ein Objekt entsprechend einer vorher festgelegten Netzform darzustellen, mit einem gewählten Netzmuster. Eine Benutzerschnittstelle mit einem mechanischen Erfassungsorgan ermöglicht, eine oder mehrere global auf das Objekt ausgeübte Kräfte zu simulieren.
Eine Speicherzone ist auch vorgesehen, zur Speicherung von:

– mechanischen Parametern des Materials bzw. Stoffs, aus dem das betreffende Objekt in dreidimensionaler Form besteht,
– Positionsdaten des Objekts, gemessen an den Spitzen eines gewählten Netzwerks, und
– Kraftdaten, die Stärke und Position einer auf das Objekt auszuübenden Belastung darstellen.

Ein Rechner, fähig mit der Speicherzone zu kooperieren, bewertet neue Positionen der Spitzen in Abhängigkeit von der ausgeübten Belastung und der mechanischen Parameter des Stoffs. Um die Verformung und/oder die Rückkehr des sich bewegenden Objekts in eine Ruhestellung sichtbar zu machen, bestimmt der Rechner eine zeitliche Entwicklung der Positionen der Spitzen in Abhängigkeit von dieser Belastung. Bei einer Anwendung zur Simulation der Verformung (oder des Durchschneidens) von Körpergewebe, kann eine qualifizierte Person (zum Beispiel ein Chirurg) die Simulation einer Verformung und/oder einer Rückkehr in eine Ruhestellung verfolgen.
Jedoch sind die simulierten Verformungen mit großen Verschiebungen üblicherweise nicht zuverlässig.
Die vorliegende Erfindung hat diese Situation verbessert.
Die Erfindung ist definiert in dem beigefügten Vorrichtungs-Anspruch 1.
In der bevorzugten Ausführungsform umfasst das gewählte Netzwerkmuster wenigstens eine Dreiecksbasis.
Wenn das Objekt also die Form einer hohlen dreidimensionalen Hülle hat, ist das gewählte Netzwerk vorzugsweise dreieckförmig, so dass der Rechner fähig ist, eine Zusammensetzung der Kräfte in jeder Dreieckspitze in Abhängigkeit von dem Unterschied zwischen der laufenden Länge jeder Kante des Dreiecks und der Länge dieser Kante im Ruhezustand zu bestimmen.
Wenn das Objekt eine volle dreidimensionale Form aufweist, ist das gewählte Netzwerk vorzugsweise tetraederförmig, so dass der Rechner fähig ist, in jeder Spitze des Tetraeders eine Kräftezusammensetzung zu ermitteln, in Abhängigkeit von dem Unterschied zwischen der laufenden Länge jeder Kante des Tetraeders und der Länge dieser Kante im Ruhezustand.
Vorzugsweise ist die Speicherzone so konzipiert, dass sie in Verbindung mit jeder Masche mechanische Parameter des Materials der Masche speichert, die wenigstens partiell lokal definiert sind, insbesondere in Höhe der Masche oder von Elementen von dieser.
In der veröffentlichten Anmeldung FR-A-2771202 des Anmelders basiert die Simulation der Verformungen des Objekts bei der beschriebenen Realisierung auf einem linearen Modell. Obwohl die Anwendung der in FR-A-2771202 beschriebenen Vorrichtung sich auf dem Gebiet der chirurgischen Assistenz bei der Simulation des Durchschneidens von Weichteilgeweben (Leber, Haut usw.) als sehr vielversprechend erwiesen hat, war die Simulation von großen, möglicherweise nicht-linearen Verschiebungen nicht komplett natürlich. Zudem erlaubte diese Vorrichtung insbesondere bei großen Verschiebungen nicht, die Reaktion auf eine auf irgend einen verformbaren Materialtyp ausgeübte Kraft vorauszusehen.
Bei einer bevorzugten Realisierung der vorliegenden Erfindung ist der Rechner fähig, den Unterschied zwischen der laufenden Länge und der Ruhezustandslänge zu bestimmen, um eine Ableitung der potentiellen Verformungsenergie jeder Masche zu schätzen, ausgedrückt in Abhängigkeit von einem Green-St.Venant-Tensor und von für das Material in dieser Masche typischen Koeffizienten.
Die Ableitung der so ausgedrückten potentiellen Energie in Bezug auf die Position der Spitze ermöglicht, die hier ausgeübte Kraft (oder die Zusammensetzung der ausgeübten Kräfte) zu erhalten.
Bei der Schätzung der also in Abhängigkeit von einem Green-St.Venant-Tensor ausgedrückten Ableitungen der potentiellen Energie an den Spitzen, ist der Rechner fähig, die Unterschiede zwischen den Quadraten der laufenden Längen und der Ruhezustandslängen jeder Kante genauer zu bestimmen, um die Zusammensetzung der vorerwähnten Kräfte zu bestimmen.
Die vorerwähnten mechanischen Parameter umfassen vorteilhafterweise die Lamé-Koeffizienten des Materials in jeder betroffenen Masche.
Die erfindungsgemäße Vorrichtung erlaubt also, große Verschiebungen und/oder nichtlineare Elastizität von möglicherweise heterogenen Objekten zu simulieren, die aus unterschiedlichen Materialien bzw. Stoffen sind.
Es ist also möglich, die Verformung eines möglicherweise in eine hohle Haut gehüllten vollen dreidimensionalen Objekts zu simulieren. Es empfiehlt sich dann, für das Objekt ein tetraederförmigen und für die Haut ein dreieckförmiges Netzwerk vorzusehen und jedem Netzwerktyp adäquate mechanische Parameter zuzuteilen. Man kann die Verschiebungen des weichen Objekts (zum Beispiel Muskeln) auch in Anwesenheit von steifen und festen Objekten (zum Beispiel Knochen) vorsehen, was weiter unten beschrieben wird.
Die vorliegende Erfindung offenbart auch ein Verfahren zur Simulation von Verformungen von Objekten, insbesondere von körperlichen Weichteilgeweben, das die folgenden Schritte ganz oder teilweise umfasst:

a) Vorsehen einer bestimmten Form des Objekts im Ruhezustand,
b) Vernetzen der genannten Form auf der Basis eines gewählten Musters, vorzugsweise auf dreieckiger Basis (Dreieck oder Tetraeder),
c) Speichern der Positionsdaten der Maschenspitzen im Ruhezustand,
d) Simulieren von einer oder mehreren global auf das Objekt ausgeübten Kräften (oder von Verschiebungen des Objekts),
e) Schätzen einer resultierenden Kraft in jeder Spitze jeder Masche in Abhängigkeit von den mechanischen Eigenschaften des Materials und den Positionen der vorhergehenden Spitzen,
f) Bestimmen der Verschiebung in Abhängigkeit von der auf jede Spitze ausgeübten Kraft und aufgrund dessen der Verformung und/oder der Rückkehr der Maschenspitzen in eine Gleichgewichtsposition.

Der Schritt c) durchläuft außerdem im Rahmen der Ausarbeitung einer Datenstruktur eine Zuteilung von Indizes für die Maschenelemente und möglicherweise eine Vereinfachung der Struktur des Objekts, wie man weiter unten sehen wird. Die Schritte e) und f) stützen sich auf ein durch den Anmelder entwickeltes Modell, das weiter unten im Detail beschrieben wird.
Weitere Merkmale und Vorteile der Erfindung gehen aus der nachfolgenden detaillierten Beschreibung und den beigefügten Figuren hervor:
die 1 stellt schematisch eine Vorrichtung gemäß einer Realisienangsform der Erfindung dar,
die 2A und 2B stellen schematisch dreieckige Maschen dar, ursprünglich bzw. verformt, mit den Bezeichnungen der Spitzen des in der Folge zur Schätzung der Zusammensetzung der Kräfte in jeder Spitze des benutzten Dreiecks,
die 3A ist ein Tetraederschema, das hilft, die für einen Punkt Pj des Tetraeders und den Normalvektor nj der entgegengesetzten Seite benutzten Bezeichnungen zu begreifen,
die 3B bis 3D sind Tetraederschemata, die zusammen mit der vorerwähnten 3A helfen zu begreifen, wie die Summen bezüglich der Kanten und Kantenpaare eines Tetraeders berechnet werden, um eine Kraft zu bewerten, der eine Spitze Pp des Tetraeders ausgesetzt ist,
die 4A und 4B sind Schemata, die helfen, die Zuteilung der Indizes bei verschiedenen dreieckigen Maschen zu verstehen,
die 5A und 5B sind jeweils Schemata von Dreiecksmustern, die eine laufende Spitze (dargestellt durch ein Quadrat) und eine virtuelle Spitze (dargestellt durch eine Scheibe) vor der Vereinfachung und die gleiche vernetzte Form nach der Vereinfachung umfassen,
die 6 zeigt schematisch ein Beispiel einer nicht-konformen Dreiecksvernetzung, und
die 7 zeigt ein schematisches Ablaufdiagramm der Berechnungen der Kräfte und der Positionen der das Objekt darstellenden Spitzen.
Die vorliegende Erfindung bezieht sich auf Formeln und Gleichungen, die als Anhang beigefügt sind.
Die Zeichnungen sowie der Anhang enthalten im Wesentlichen Elemente von feststehender Eigenart. Sie können also nicht nur dem besseren Verständnis der Beschreibung dienen, sondern gegebenenfalls auch zur Definition der Erfindung beitragen.
Zunächst wird auf die 1 Bezug genommen, in der die Vorrichtung sich in der beschriebenen Realisierung dargestellt ist, in Form eines Computers des PC-Typs. Der Computer enthält einen nichtflüchtigen Speicher MEM (zum Beispiel des Typs ROM), in den die Programme zur Berechnung der Maschenspitzenpositionen und die resultierenden Kräfte in diesen Spitzen eingespeichert sind, sowie einen Prozessor μP, fähig mit einem Arbeitsspeicher MT (zum Beispiel des Typs RAM) zu kooperieren, um diese Berechnungen durchzuführen.
Die weiter oben beschrieben Module des Rechners können gegebenenfalls in dem Speicher MEM abgespeicherte oder von einem nicht dargestellten Peripheriegerät, zum Beispiel einem Disketten- oder CD-ROM-Lesegerät, in den Arbeitsspeicher MT geladene Programme sein, oder auch mittels MODEM empfangen werden. Selbstverständlich sind solche Programme besonders wichtige Mittel zur Anwendung der Erfindung. Daher betrifft die vorliegende auch solche Programme – möglicherweise auf einem Träger des vorerwähnten Typs (Diskette, CD-ROM usw.) – von denen weiter unten detaillierte Realisierungsbeispiele angegeben werden.
Der Computer umfasst außerdem einen Monitor VISU, verbunden mit dem Prozessor μP durch eine graphische Schnittstelle IG, wobei die Einheit aus VISU und IG (oben "Visualisierungsschnittstelle" genannt) ermöglicht, die Form eines Objekts sichtbar zu machen, dessen Verformung simuliert wird. Der Computer umfasst noch eine Benutzerschnittstelle IU, verbunden mit einem mechanischen Erfassungsorgan, zum Beispiel einer Tastatur CLA oder einer Maus MO, die einem Benutzer ermöglichen, eine global auf ein Objekt ausgeübte Kraft zu simulieren, oder auch eine Verschiebung, mit Anfangsverformung, dieses Objekts.
Die Form des zu vernetzenden Objekts kann eine in drei Dimensionen dargestellte Fläche sein. In diesem Fall ist die Vernetzung dreieckförmig (2A). Die Form des zu vernetzenden Objekts kann dreidimensional und voll sein, und in diesem Fall ist die Vernetzung tetraederförmig (2B). Möglicherweise ist ein Vernetzungsmodul in dem Speicher MEM abgespeichert und kann durch den Prozessor μP angewendet werden, indem er mit der graphischen Schnittstelle IG kooperiert, um eine Form des vernetzten Objekts auf dem Monitor VISU sichtbar zu machen.
In der Folge wird ein durch den Anwender entwickeltes und bei einer Verformungssimulationsvorrichtung angewandtes Modell beschrieben. Bei diesem Modell wird eine Ruhezustandsform betrachtet, ohne Argumentation in Bezug auf eine Ruheposition (Verschiebungsargumentation, die üblicherweise die Basis der klassischen Modelle ist).
Zunächst wird eine in drei Dimensionen dargestellte, entsprechend einem Dreiecksmuster vernetzte Fläche betrachtet. Ein Dreieck T0 = (P0, P1, P2) mit der Ruhezustandsform, mit Kanten (I₀, I₁, I₂). Dieses Dreieck (2A) erfährt eine Verformung Φ, so dass es zum Dreieck T = (Q₀, Q₁, Q₂) wird (2B).
Die potentielle elastische Energie dieser Verformung wird hier durch den Verformungstensor C von Cauchy-Green (Formel F1 im Anhang) ausgedrückt. Bei Betrachtung des finiten bzw. endlichen Elements, erlangt indem man diesem Dreieck eine Gruppe linearer Interpolationsfunktionen zuordnet, drückt sich der Gradient der Verformung in jedem Punkt X des Dreiecks durch die Formel 2 aus, wo A die Grundfläche des Dreiecks T0 ist, nj der Einheitsvektor ist, senkrecht zum Segment (Pj + 1, Pj + 2) und die αj die Formvektoren des Dreiecks sind.
Der Verformungstensor C wird dann durch den Ausdruck der Formel F3 geliefert. Vorteilhafterweise ist dieser Verformungstensor invariabel durch Translation und hängt folglich nicht von der Wahl des Nullpunkts des Bezugssystems ab. Vorzugsweise befindet sich der Nullpunkt des Bezugssystems im Zentrum des Umfangskreises des verformten Dreiecks T mit dem Radius r (2A), was ermöglicht, den Tensor C in der in der Formel 4 angegebenen Form auszudrücken, wo Lp die Länge der verformten Kante ist, welche die Spitzen Qj und Qk verbindet.
Der Ausdruck der potentiellen Verformungsenergie (geliefert durch die Formel 5 des Anhangs) benutzt die des Verlaufs des Tensors C (geliefert durch die Formel F6) und des Verlaufs des Tensors C² (der weiter unten beschrieben wird).
Ausgehend von den klassischen Relationen bei dem Dreieck T0 (Kosinusgesetz), ausgedrückt durch die Formel F7, erhält man die Skalarprodukte zwischen den Richtungsvektoren, geliefert durch die Relation R1 des Anhangs. Die jeweiligen Verläufe der Tensoren C und E schreiben sich dann entsprechend den Formeln 8 und 9 des Anhangs. Dieser Term hängt nur von den Längenveränderungen der Kanten des Dreiecks und von seiner Form im Ruhezustand ab (dargestellt durch die Größen Ki).
Von dem Ausdruck von tr C² (Formel F10) leitet man trE² ab und schließlich wird die potentielle elastische Verformungsenergie des Dreiecks durch die Relation R2 des Anhangs geliefert.
Die durch das Dreieck auf eine seiner Spitzen ausgeübte elastische Kraft erhält man, wenn man die elastische Energie in Bezug auf die Position dieser Spitze ableitet. In dem in der 2A dargestellten Beispiel hat die Spitze P0 als Verformungsposition Q0. Wenn man die Relation R2 in Bezug auf Q0 ableitet, beginnt man vorzugsweise damit, die Terme L²p abzuleiten, um den Ausdruck der Kraft entsprechend der Formel F11 zu erhalten, welche die Gleichung der durch das Dreieck auf jede seiner Spitzen ausgeübte elastischen Kraft ist. Diese Kraft drückt sich sehr wohl nur in Abhängigkeit von den Längenänderungen der Kanten des Dreiecks (L_i ² – I_i ²), der mechanischen Parameter des Materials (λ und μ) und der Form des Dreiecks im Ruhezustand aus (geliefert durch die Terme Ki).
Vorteilhafterweise enthält diese Kraft implizit eine Eigenschaft der Invarianz durch Rotation, da alle ihre Terme null bleiben, solange die Längen der Kanten des Dreiecks sich nicht verändern, das heißt solange das Dreieck sich nicht verformt. Die Verschiebungsberechnungen können sich also aufgrund der Tatsache vereinfachen, dass es nicht mehr notwendig ist, die starren Transformationen zu kompensieren, indem die Invarianz durch starre Transformation bei der Formulierung des verformbaren Modells intrinsisch ist.
Indem man das Vorhergehende ausdehnt auf eine volle dreidimensionale Form (Fall eines Volumens), vernetzt mit Tetraedem, ermöglicht die Formulierung der Elastizität nach Saint Venant-Kirchhoff, den Deformiertengradienten auszudrücken (Formel F12) sowie den Verformungstensor von Cauchy-Green (Formel F13), wo Aj die Grundfläche der der Spitze Pj entgegengesetzten Seite ist, nj die Normale zu dieser Seite bzw. Fläche ist und V das Volumen des Tetraeders ist (3A).
Indem man den Nullpunkt des Bezugssystems im Zentrum der den Tetraeder T umschließenden Kugel anordnet, kann man das Skalarprodukt der Positionsvektoren von zwei Spitzen des Tetraeder in Abhängigkeit von der Länge der sie verbindenden Kante und dem Radius r der Kugel ausdrücken, was ermöglicht, den Verformungstensor gemäß der Formel F14 zu schreiben.
Die elastische Verformungsenergie drückt sich in Abhängigkeit von dem Verformungstensor von Green-Saint Venant (Referenz E) aus, aber bei einer 3D-Masche (Formel F15). Der Ausdruck des Verlaufs des Tensors C drückt sich wieder durch die Skalarprodukte aus (αj·αk), indem das Kosinusgesetz für die Flächen- bzw. Raumwinkel benutzt wird, und in Abhängigkeit von den Längen der Kanten des Tetraeders (Formeln 16 und 17). Auch der Verlauf von C² schreibt sich entsprechend der Formel 18. Schließlich zeigt der Anmelder, dass die elastische Verformungsenergie sich durch die Formel F19 des Anhangs (ergänzt durch die Formel 20) ausdrückt.
Den Ausdruck der durch das Tetraeder auf eine seiner Spitzen ausgeübte Kraft erhält man, indem man die elastische Energie in Bezug auf die Position dieser Spitze ableitet. Eine leicht anwendbare Formulierung wird durch die Formel F21 des Anhangs ausgedrückt und lässt mit Hilfe der 3A bis 3D interpretieren.
In der Formel F21:

– sind in der zweiten Zeile k und m die komplementären Indizes von j und p,
– ist in der letzten Zeile m der komplementäre Index von j, k und p, und: K_j _kmp = K_jpK_km + K_jmK_kp – K_jkK_mp.

Die Summe bezüglich der Kanten (j, p) bedeutet, dass man die drei in der Spitze Pp (auf die die Kraft ausgeübt wird) zusammenlaufenden Kanten durchlaufen muss, während bei der Summe bezüglich der Kanten (j, k) (von p unterschiedliche Indizes) die drei auf der Gegenseite zur Spitze Pp befindlichen Kanten beteiligt sind.
Mit Bezug auf die 3A bis 3D:

– bei der Summe bezüglich der Kanten (j, p) sind drei Terme beteiligt (3A),
– bei der Summe bezüglich der Kantenpaare ((j, p), (k, p) sind drei Terme beteiligt (3B),
– bei der Summe bezüglich der Kantenpaare ((j, p), (j, k) sind sechs Terme beteiligt (3C) und
– bei der Summe bezüglich der Kantenpaare ((j, k), (m, p) sind drei Terme beteiligt (3D).

Es zeigt sich wieder eine Rotations-Invarianz der vorhergehenden Ausdrücke für eine volle dreidimensionale Form, was vorteilhafterweise ermöglicht, das Auftreten von irrealistischen Distorsionen zu vermeiden, wenn ein Teil des Objekts eine zu großen Rotation erfährt. Diese Invarianz erstreckt sich auch auf die Translation. Das Objekt kann also eine beliebige Translations- und Rotationsbewegung ausführen und dabei eine Energie und interne Kräfte beibehalten, die null sind. Der Anmelder hat verifiziert, dass in diesem Fall die starre Komponente der Verformung null ist.
Zudem läuft die Aufgliederung des Objekts, indem man finite Elemente verwendet, darauf hinaus, die elastische Energie einer Gruppe von Teilbereichen, definiert durch einfache geometrische Formen, zu schätzen. So schlägt das oben vorgeschlagene Modell für jeden Teilbereich einer Energieformulierung vor, die nur von der Formänderung dieses Bereichs abhängt.
Vorteilhafterweise erscheint die Sichtbarmachung einer simulierten Verformung, bei der man eine anfängliche Ruheform (eher als anfängliche Ruhepositionen) betrachtet, natürlicher als die mit bekannten Modellen simulierten Verformungen. Zudem vermeidet die Betrachtung einer Ruheform die Genauigkeitsprobleme, die auftreten, wenn das Modell sich deutlich von seiner Ruheposition entfernt.
In der Folge werden Vorbereitungen werden der Anwendung der vorerwähnten Module des Rechners zur Bestimmung der Kräfte in jeder Spitze des gewählten Netzwerks beschrieben.
Es empfiehlt sich, zunächst eine Datenstruktur zu definieren. Der Anmelder war zunächst mit der Schwierigkeit konfrontiert, den Elementen jeder Masche adäquate Indizes zuzuteilen.
Für Dreiecksmuster besteht die Datenstruktur, die einer Triangulation entspricht, aus einer Dreiecksliste, einer Kantenliste und einer Spitzenliste.:

– für jede Spitze bestimmt man eine angrenzende Kante,
– für jede Kante bestimmt man einen Pointer in den beiden angrenzenden Dreiecken und den beiden die Kante bildenden Spitzen,
– für jedes Dreieck bestimmt man die Pointer in den drei Spitzen, die das Dreieck bilden, die drei Kanten des Dreiecks, die drei an jede Kante angrenzenden Dreiecke sowie drei Werte, die den Nachbarschaftsindex des laufenden Dreiecks für jedes benachbarte Dreieck liefern.

In der Folge wird die Ausrichtung der drei Spitzen in dem Dreieck beschrieben.
Damit die globale Vernetzung korrekt ausgerichtet ist, muss jede Kante in zwei unterschiedlichen Richtungen durchlaufen werden, wenn man jedes der benachbarten Dreiecke durchläuft. Mit Bezug auf die 4A ist die Kante Ai mit Index i die der Spitze Si mit Index i entgegengesetzte Kante. Das benachbarte Dreieck mit Index i ist das an die Kante Ai angrenzende Nachbardreieck TV. Schließlich hat jedes Dreieck immer drei Nachbardreiecke. Um Flächen mit Rändern darzustellen (eine Ebene hat einen Rand, ein Zylinder hat zwei Ränder) benutzt man den Begriff so genannter "virtueller" Dreiecke (4B). Die virtuellen Dreiecke werden nicht angezeigt bzw. sichtbar gemacht und die Position der virtuellen Spitze hat keine Auswirkung auf die Geometrie der Vernetzung.
Es folgt nun ein Beispiel einer Index-Zuteilung in C++-Sprache.
Wenn man das benachbarte Dreieck Nummer i Nachbar[i] nennt, dann wird der Nachbarindex "whichNeighbor[i]" durch eine Instruktion definiert, etwa:
Neighbor[i]->neighbor[whichNeighbor[i]]=this;
Bei tetraederförmigen Vernetzungen besteht die Datenstruktur eines Tetraedervolumens aus einer Tetraederliste, einer Dreiecksliste, einer Kantenliste und einer Spitzenliste.
Für jede Spitze speichert man ihre Position und einen Nachbartetraeder.
Für jede Kante speichert man die Pointer in den benachbarten Spitzen und einen Pointer in einem dieser Kante benachbarten Dreieck.
Für jeden Tetraeder speichert man:

– die Pointer in den vier Spitzen, die den Tetraeder bilden,
– die Pointer in den sechs Kanten des Tetraeders,
– die Pointer in den vier Dreiecken des Tetraeders, und
– die Pointer in den vier benachbarten, jeder der vier Spitzen entgegengesetzten Tetraedern.

Jeder Tetraeder ist ausgerichtet. Ausgehend von der Ausrichtung der vier Spitzen ABCD sucht man die Ausrichtung der vier Dreiecke BCD, ADC, ABD, ACB. Damit die Vernetzung richtig ausgerichtet ist, muss jedes Dreieck in zwei unterschiedlichen Richtungen, indem jeder dieser benachbarten Tetraeder durchlaufen wird. Schließlich ist jeder Tetraeder immer vier Tetraedern benachbart. Man führt die Begriffe der sogenannten "virtuellen Spitzen" und der "virtuellen Tetraeder" ein, um den Rändern der Vernetzung Rechnung zu tragen.
Um die Struktur des verformten Objekts zu vereinfachen, zerschneidet man die Vernetzung, indem man Material wegnimmt (oder genauer eine Spitze, mit der Auflage, dass diese Spitze keine virtuelle Spitze sein kann).
Im Falle einer Dreiecksvernetzung – wenn diese Spitze keiner virtuellen Spitze benachbart ist –, betrachtet man diese Spitze sowie alle diese benachbarten Dreiecke als virtuell. Wenn die Spitze schon einer virtuellen Spitze benachbart ist, entfernt man die Kante, die diesen Punkt mit der virtuellen Spitze verbindet und zwingt alle der virtuellen Spitze benachbarten Dreiecke, virtuell zu sein (5A und 5B).
Bei einem tetraederförmigen Vernetzungsmuster fährt man fort, indem man eine Gruppe miteinander verbundener Tetraeder entfernt. Wenn diese Tetraeder entfernt werden, muss die Vernetzung konform bleiben, wie unten beschrieben. Dazu kann es vorteilhaft sein, mehr Tetraeder als die zu entfernen, welche zu der anfangs erwähnten Gruppe gehören. Man kann auch virtuelle Spitzen hinzufügen oder wegnehmen, wenn man zusammenhängende Komponenten schafft oder wenn man Ränder beseitigt.
Der Anmelder hat sich der Schwierigkeit angenommen, eine Konformität der Vernetzung, wie anschließend definiert, sicherzustellen.
Damit eine Dreiecksnetzwerk konform ist, darf jede Kante zu höchstens zwei Dreiecken gehören und die topologische Nachbarschaft jeder Spitze muss homöomorph bezüglich eines Kreises (oder eines Halbkreises) bleiben, was bedeutet, dass – wenn zwei Dreiecke sich eine selbe Spitze teilen, zwischen ihnen eine Gruppe benachbarter Dreiecke existieren muss, die den Übergang von einem Dreieck zum einem Anderen ermöglichen. Das Dreiecksnetzwerk der 6 ist nicht konform, denn die beiden Dreiecke haben eine Spitze als Schnittpunkt, so dass keine benachbarten Dreiecke existieren, die einen Übergang von einem dieser Dreiecke zum anderen ermöglichen.
Damit ein Tetraedernetzwerk konform ist, darf jedes Dreieck zu höchstens zwei Tetraedem gehören und die topologische Nachbarschaft jeder Spitze muss homöomorph bezüglich einer Kugel (oder einer Halbkugel, wenn sie sich am Rand befindet) sein. Wenn also zwei Tetraeder sich eine selbe Spitze teilen, muss eine Gruppe benachbarter Tetraeder existieren, die den Übergang von einem Tetraeder zu einem anderen ermöglichen.
Die Gleichung, welche die Bewegung des linearen elastischen Modells bestimmt, ist die newtonsche Differentialgleichung E1 des Anhangs. Diese Gleichung ist verknüpft mit den Differentialgleichungen, denen man in der Kontinuum-Mechanik begegnet (Gleichung E2).
Nach der Theorie der finiten Elemente sind die Massenmatrix M und die Dämpfungsmatrix C hohl und hängen von den Eigenschaften jedes Tetraeders ab. Zudem können diese Matrizen M und C hier als Diagonale betrachtet werden, das heißt, dass die Masse und die Dämpfungseffekte in den Spitzen des Netzwerks konzentriert sind. Diese Vereinfachung, masse-lumping genannt, entkoppelt die Bewegungen jedes Knotenpunkts und ermöglicht folglich, die vorhergehende Gleichung E2 wie eine Gruppe unabhängiger Differentialgleichungen des Typs E1 zu schreiben, mit einer Gleichung für jede Spitze.
Die Wahl eines expliziten Integrationsschemas führt dazu, die elastischen Kräfte zum Zeitpunkt t auszudrücken, um mit der Gleichung E3 des Anhangs die neue Position der Spitze zum Zeitpunkt t + 1 zu berechnen. Vorteilhafterweise ist diese Gleichung E3 gültig für jede Art von Kraft, linear oder nicht.
Die expliziten Schemata sind nur bedingt stabil; sie konvergieren folglich weniger schnell als die impliziten Schemata.
Eine weitere Schwierigkeit ist aufgetreten: damit das Schema stabil ist, muss der zur numerischen bzw. digitalen Integration benutzte Zeitschritt ausreichend klein sein. Der kritische Zeitschritt hängt ab von den lokalen Werten der Masse und der Dämpfung. Um den Zeitschritt zu optimieren, wird die Massematrix als gleich der Identität festgelegt und die Optimierung der lokalen Dämpfungswerte – um die Oszillationsrisiken zu begrenzen – hat die Wirkung, die Konditionierung des linearen Systems zu verbessern und folglich dessen Auflösung zu stabilisieren.
Die 7 zeigt ein Ablaufdiagramm, das schematisch und beispielartig die verschiedenen, auf eine aus Dreiecken zusammengesetzten dreidimensionalen Fläche angewandten Verarbeitungsschritte darstellt.
Die Verarbeitung beginnt mit einem ersten Schritt 10, in dem man die oben beschriebene Datenstruktur appliziert, mit einer Zuteilung adäquater Indizes. Man betrachtet anschließend im Schritt 12 jedes Dreieck T_n, für das man die Lamé-Koeffizienten λ_n und μ_n sowie die Grundfläche des Dreiecks A sucht (Schritt 14). Für jede Kante des Dreiecks (j sei eine ganze Zahl von 0 bis 2), berechnet man die Differenz der Quadrate der Kantenlängen zwischen der laufenden Länge l_jund der vorhergehenden Kantenlänge L_j (Schritt 16). Von diesen Differenzen der Kantenlängen im Quadrat leitet man die in jeder Spitze S_j des Dreiecks T_n ausgeübten Kräfte ab (Kraft F_j ⁿ, berechnet im Schritt 18 des Ablaufdiagramms). Indem man zu den nächsten Dreiecken übergeht, integriert man die Kräfte aller Dreiecke mittels einer Index-n-Schleife (Schritt 20).
In den Schritten 30 und 32 bestimmt man die neuen Positionen Q'_n der Spitzen S_j der verformten Dreiecke, ausgehend von ihrer vorhergehenden Position Q_n und den Kräften Fⁿ. Mittels einer Index-n-Schleife (Schritt 36) realisiert man eine Integration aller verformten Objekte. Die Schätzung der neuen Positionen der Spitzen S_j erfolgt vorzugsweise, indem man die vorerwähnten Gleichungen E1, E2 und E3 benutzt.
In dem Ablaufdiagramm der 7 sind die beiden jeweiligen Integrationen getrennt dargestellt (Kräfte F_j ⁿ und Positionen der Spitzen S_j in den verformten Objekten). Selbstverständlich kann vorgesehen werden, einige der Schritte 12 bis 20 gleichzeitig mit einigen der Schritte 30 bis 36 durchzuführen oder danach. Diese Wahl hängt ab von der Struktur der beiden oben beschriebenen Module, die der vorerwähnte Rechner umfasst, und möglicherweise von der Struktur der Daten.
Vorzugsweise wird ein Test 42 vorgesehen, um zu wissen, ob das Objekt ganz in seine Ruheform zurückgekehrt ist, nämlich: gibt es eine Kante in dem Objekt, deren Länge sich wesentlich von der vorhergehenden Länge unterscheidet?
Wenn der Rechner beim Durchlaufen der Dreiecke keine solche Kante findet, bedeutet das, dass das Objekt in seine ursprüngliche Ruheform zurückgekehrt ist; die Verarbeitung endet dann in 40.
Wenn der Rechner aber beim Durchlaufen der Dreiecke eine solche Kante findet, wird ein Bild des Objekts durch den Monitor VISU angezeigt und eine Zeitschleife (Schritt 44) wird inkrementiert, um zu einer nächsten Zeit überzugehen und eine progressive Rückkehr zu einer Ruheform des Objekts zu zeigen, indem man ggf. die Gleichung E3 anwendet. Diese Zeitschleife kann nach dem ersten Schritt oder Anfangsschritt 10 initialisiert werden.
Selbstverständlich ist die Erfindung nicht auf die oben beispielartig beschriebene Realisierungsform beschränkt, sondern erstreckt sich auf weitere, im Rahmen der beigefügten Patentansprüche definierte Varianten.
So ist verständlich, dass die Erfindung sich nicht auf eine Anwendung zur Simulation der Verformung von körperlichen Weichgeweben beschränkt. Wie oben beschrieben, kann die Verschiebung eines starren Objekts (zum Beispiel bei einem mit einem Muskel verbundenen Knochen), das mit Weichteilgewebe verbunden ist, auch simuliert werden. Gegebenenfalls genügt es, anstatt eines Verformungstensors des oben beschriebenen Typs eher einen Torsor (an sich bekannten Ausdrucks) zu benutzen.
Die Verformungsenergieformel W sowohl einer dreieckigen als auch einer tetraederförmigen Masche wurde oben beispielartig angegeben und hängt von dem Material des Objekts ab, dessen Verformung simuliert wird. Jedoch wird sie allgemein durch eine Kombination der Verläufe des Tensors E und des Tensors E² ausgedrückt. In besonderen Fällen kann es sein, dass nur der Verlauf des Tensors E (oder des Tensors E²) eine Rolle spielt.
ANHANG

Claims

Vorrichtung zur Materialverformungssimulation, insbesondere von weichen Körpergeweben, umfassend: – eine Speicherzone (MEM, MT), fähig von einem Objekt mit dreidimensionaler Form abzuspeichern: – mechanische Parameter des Materials des Objekts (λ, μ), – Positionsdaten des Objekts (Lj), gemessen an den Spitzen (Sj) wenigstens eines gewählten Netzwerks (Tn), und – Kraftdaten (F_j ⁿ), welche Stärke und Position einer auf das Objekt auszuübenden Belastung darstellen, – sowie einen Rechner (μP, MT), fähig mit der Speicherzone (MEM) zu kooperieren, um neue Positionen der Spitzen (Sj) zu ermitteln, in Abhängigkeit von der ausgeübten Belastung und den mechanischen Parametern des Materials, dadurch gekennzeichnet, dass der Rechner umfasst: – einen Modul (12, 14, 16, 18, 20), fähig für jede Masche den Unterschied zwischen der laufenden Länge (lj) einer Kante der Masche und ihrer vorhergehenden Länge und/oder Ruhezustandslänge (Lj) iterativ zu ermitteln sowie für jede Spitze der Masche jeweilige, die potentielle Verfomungsenergie betreffende Kraftdaten (F_j ⁿ) abzuspeichern, wobei man die jeweiligen Kraftdaten aufgrund des genannten Unterschieds erhält, – einen Modul (30, 32, 36), fähig für jede Spitze neue Positionsdaten dieser Spitze (Sj(Qn)) in Abhängigkeit von der Zusammensetzung der dort ausgeübten Kräfte (F_j ⁿ), von Daten bezüglich wenigstens einer früheren bzw. vorhergehenden Position der Spitze (Sj(Tn)) und von mechanischen Parametern des Materials iterativ zu ermitteln.
Vorrichtung nach Anspruch 1, dadurch gekennzeichnet, dass für ein Objekt mit der Form einer hohlen dreidimensionalen Hülle das gewählte Netzwerk dreieckig ist, so dass der Rechner fähig ist, in jeder Spitze des Dreiecks eine Kräftezusammensetzung zu ermiteln, in Abhängigkeit von dem Unterschied zwischen der laufenden Länge jeder Kante des Dreiecks und der Länge dieser Kante im Ruhezustand.
Vorrichtung nach einem der Ansprüche 1 und 2, dadurch gekennzeichnet, dass für ein Objekt mit der Form einer vollen dreidimensionalen Hülle das gewählte Netzwerk tetraederförmig ist, so dass der Rechner fähig ist, in jeder Spitze des Tetraeders eine Kräftezusammensetzung zu ermitteln, in Abhängigkeit von dem Unterschied zwischen der laufenden Länge jeder Kante des Tetraeders und der Länge dieser Kante im Ruhezustand.
Vorrichtung nach einem der vorangehenden Anspreche, dadurch gekennzeichnet, dass der Rechner fähig ist, Differenzen zwischen den Quadraten der laufenden (lj) und der vorhergehenden Länge und/oder Ruhezustandslänge jeder Kante zu ermitteln, um die genannte Kraftezusammensetzung zu bestimmen.
Vorrichtung nach einem der vorangehenden Ansprüche, dadurch gekennzeichnet, dass die Speicherzone fähig ist, für jede Masche, speziell von der gesamten Masche oder von Elementen von dieser, mechanische Parameter des Materials der Masche (λ_n, μ_n) zu speichern, die zumindest teilweise lokal definiert sind.
Vorrichtung nach Anspruch 5, dadurch gekennzeichnet, dass der Rechner fähig ist, den genannten Unterschied zwischen laufender Länge und Ruhezustandslänge zu ermitteln, um eine Ableitung der potentiellen Verformungsenergie (W) jeder Masche in Bezug auf die Position jeder Maschenspitze zu schätzen, wobei die potentielle Verformungsenergie ausgedrückt wird in Abhängigkeit von einem Green-St.Venant-Tensor (E) und von für das Material in dieser Masche typischen mechanischen Koeffizienten, was die genannten Kraftdaten für die Spitze der Masche liefert.
Vorrichtung nach Anspruch 6, dadurch gekennzeichnet, dass die genannten mechanischen Parameter die Lamé-Koeffizienten des Materials in jeder betrachteten Masche (λ_n, μ_n) umfassen.
Vorrichtung nach einem der vorangehenden Ansprüche, dadurch gekennzeichnet, dass sie einen Modul umfasst, der dazu dient, eine Datenstruktur auszuarbeiten und Maschenkanten (Ai) wegfallen zu lassen, die zwei sogenannte "virtuelle" Spitzen (Si) verbinden.
Vorrichtung nach Anspruch 8, dadurch gekennzeichnet, dass der Modul zur Ausarbeitung einer Datenstruktur außerdem fähig ist, zu verifizieren, dass das Netzwerk vorher definierten Konformitätseigenschaften entspricht.
Vorrichtung nach einem der vorangehenden Ansprüche, dadurch gekennzeichnet, dass der genannte Modul (30, 32, 36), der fähig ist, die neuen Positionsdaten der Spitzen (Q'n) in Abhängigkeit von der Kräftezusammensetzung in jeder Spitze zu ermitteln, so konzipiert ist, dass er die genannten neuen Positionsdaten in Abhängigkeit von der Zeit (E3) ermittelt, was ermöglicht, eine zeitliche Entwicklung der jeweiligen Positionen der Spitzen zu verfolgen.
Vorrichtung nach Anspruch 10, dadurch gekennzeichnet, das die Ermittlung der neuen Positionsdaten der Spitzen (Q'n) ein auf jede Spitze angewendetes Differentialgleichungsauflösungsmodell benutzt.
Vorrichtung nach einem der vorangehenden Ansprüche, dadurch gekennzeichnet, dass der Rechner fähig ist, hinsichtlich der Ermittlung einer zeitlichen Entwicklung (t + Δt) der genannten Positionen die Positionsdaten der Spitzen des Netzwerks wiederholt zu ermitteln.
Vorrichtung nach Anspruch 12, dadurch gekennzeichnet, dass sie ein Visualisierungs-Interface (VISU, IG) umfasst, fähig das Objekt entsprechend einer vorher festgelegten Form und maschig gemäß einem gewählten Muster darzustellen, und dadurch, dass das Visualisierungs-Interface fähig ist, mit dem Rechner zu kooperieren, um die Form des Objekts in Bewegung darzustellen.
Vorrichtung nach einem der vorangehenden Ansprüche, dadurch gekennzeichnet, dass sie ein Benutzer-Interface (IU) umfasst, ausgestatte mit einem Erfassungs- bzw. Zugriffsorgan (CLA, MO), um eine oder mehrere global auf das Objekt ausgeübte Kräfte zu simulieren.