DE69513300T2

DE69513300T2 - Bestimmungsverfahren

Info

Publication number: DE69513300T2
Application number: DE69513300T
Authority: DE
Inventors: Knut Conradsen; Knud Fredfeldt; Michael Grunkin; Lars Hyldstrup
Original assignee: TORSANA A S SKODSBORG
Current assignee: TORSANA A S SKODSBORG
Priority date: 1994-08-29
Filing date: 1995-08-21
Publication date: 2000-03-23
Anticipated expiration: 2015-08-22
Also published as: AU3341295A; US5915036A; WO1996007161A1; EP0905638A1; DE69513300D1; EP0777892A1; EP0777892B1

Description

Die vorliegende Erfindung betrifft das Gebiet der nichtinvasiven Bestimmung der Knochengualität von Vertebraten. Diese Bestimmung kann bei der Diagnose von z. B. Osteoporose und anderen Knochenerkrankungen verwendet werden, die eine Knochenbrüchigkeit hervorrufen und somit das Risiko einer Knochenfraktur erhöhen.
Beschleunigter Knochenverlust, der zur Osteoporose führt, ist ein häufig anzutreffendes Phänomen bei Frauen nach der Menopause. Frauen leiden oftmals unwissentlich an einem beschleunigten Knochenverlust, und die verminderte Festigkeit der Knochen wird nicht entdeckt, bis ein Knochen gebrochen ist oder ein Vertebrat aufgrund einer Last zusammenbricht, der ein gesunder Knochen oder Vertebrat widerstehen können sollte. Somit ist sich ein großer Teil der Osteoporose- Patienten der verminderten Festigkeit ihrer Knochen nicht bewußt, bis eine Fraktur die Krankheit und deren Ausmaß zu Tage bringt.
Heutzutage betrifft die Osteoporose mehr als ein Drittel der älteren Frauen in dem industrialisierten Teil der Welt. Die Verbreitung dieser Krankheit ist weiterhin ansteigend, was teilweise durch den proportionalen Anstieg älterer Menschen und teilweise aus unerklärlichen Gründen hervorgerufen wird. In den weniger entwickelten Teilen der Welt, wie in Südostasien und Südamerika, wird vorhergesagt, daß die Osteoporose innerhalb der nächsten 20-30 Jahre zu einer enormen sozio-ökonomischen Last führen wird. Wenn jedoch Individuen identifiziert werden können, die gefährdet sind, Osteoporose zu entwickeln, können Präventivmaßnahmen vorgenommen werden. Dies erfordert ein zuverlässiges, kostengünstiges und sicheres Verfahren für die Identifikation solcher risikobehafteter Individuen.
Da das Risiko einer Osteoporose eng mit der Knochenwiderstandsfähigkeit gegen Frakturen verbunden ist, können Knochenfestigkeits- oder Knochenqualitätsmessungen wahrscheinlich die notwendige Information bereitstellen.
Da die Mittel zum Verhindern von Osteoporose sehr viel effizienter als die zum Behandeln von Osteoporose sind, ist die Identifikation von Individuen, die gefährdet sind, diese Krankheit zu entwickeln, entscheidend. Lediglich durch eine frühe Prävention von Osteoporose können die individuellen sowie die sozio-ökonomischen Folgen der Krankheit minimiert werden.
Das erfindungsgemäße Verfahren schafft ein Verfahren zum Detektieren von Abnahmen in der Knochenqualität, wie der Architektur und der Festigkeit, in einem frühen Stadium. Somit schafft das Verfahren der Erfindung ein Werkzeug, potentielle Patienten auf Knochenerkrankungen zu überprüfen, und schafft somit die Möglichkeit einer frühen Detektion von Knochenerkrankungen in einem Stadium, in dem noch keine Krankheitssymptome von dem Patienten bemerkt werden.
Ein typisches Maß für die Knochenfestigkeit ist die Knochenmineraldichte (KMD) des Knochens gewesen. KMD-Messungen werden typischerweise erhalten, indem ein Knochen zusammen mit einem Standardkeil geröntgt wird. Nachdem bestimmt wurde, welcher Teil des Keiles die Röntgenstrahlen im gleichen Ausmaß wie der bestrahlte Knochen dämpft, kann ein Maß der KMD erhalten werden.
Ein differenzierteres Verfahren zum Bestimmen der KMD ist die Röntgenstrahl-Absorptionsmeßtechnik bei zwei unterschiedlichen Wellenlängen. Durch Verwendung von zwei Wellenlängen kann das Verfahren die Effekte von weichem Gewebe etc. um den Knochen herum kompensieren. Somit kann mit dem Verfahren eine exaktere Bestimmung der Röntgenstrahldämpfung des Knochens erhalten werden.
Wie jedoch aus dem folgenden deutlich wird, ist die Knochenmineraldichte eines Knochens nicht notwendigerweise mit der tatsächlichen Festigkeit des Knochens verbunden. Der Grund hierfür liegt in dem Aufbau eines Knochens eines Vertebraten.
Ein Knochen eines Vertebraten besteht aus einer kortikalen äußeren Schicht und einer abbaubaren ("cancellous") inneren Struktur. Ein Weglassen der abbaubaren inneren Struktur eines Knochens würde zu einem ziemlich gebrechlichen Knochen führen. Die abbaubare innere Struktur des Knochens besteht aus sog. Trabekulae. Dickere vertikale Trabekulae sind in dem Knochen in Richtung der Hauptlast (Hauptkompression oder -zug in dem Knochen) angeordnet und dünnere horizontale Trabekulae verbinden die vertikalen Trabekulae miteinander. Somit setzt sich die Hauptdichte eines Knochens aus der kortikalen Schicht und den vertikalen Trabekulae zusammen. Dies ist ziemlich natürlich, da der Knochen den größten Teil des Materials dafür bereithält, die am häufigsten auftretenden Lasten auszuhalten. Eine KMD-Messung, welche die Dichte eines Knochens bestimmt, wird daher primär den Knochengehalt in der kortikalen Schicht und in den vertikalen Trabekulae schätzen und daher nur zu einem geringen Maß die Fähigkeit des Knochens, die Lasten auszuhalten, denen der Knochen ausgesetzt ist.
Wenn jedoch Lasten nicht in Richtung der vertikalen Trabekulae angelegt werden, würde die Struktur des Knochens ohne die dünneren horizontalen Trabekulae, welche die vertikalen Trabekulae miteinander verbinden, gebrechlich sein. Die horizontalen Trabekulae definieren die Fähigkeit des Knochens, Lasten nicht in Richtung der vertikalen Trabekulae auszuhalten. Somit wird die Festigkeit des Knochens in Richtungen, die nicht in Richtung der vertikalen Trabekulae verlaufen, durch einen sehr geringen Anteil der Knochenmasse definiert.
Verschiedene Knochen variieren in der Relation zwischen der kortikalen und der trabekularen Knochenstruktur. Ferner variiert die Architektur der Trabekulae entsprechen der potentiellen Belastung dieses Teils des Knochens. Die Anzahl und der Durchmesser longitudinaler Trabekulae und das Vorliegen transversaler Trabekulae sind für die Brüchigkeit des Knochens ebenfalls von fundamentaler Wichtigkeit. Diese Unterschiede sind für die Festigkeit des Knochens äußerst wichtig, sie können jedoch lediglich zu einem geringen Ausmaß durch KMD-Messungen detektiert werden.
Da der trabekulare Knochen ein größeres Oberflächen/Gewichts-Verhältnis hat, ist er in höherem Maße einem beschleunigten Knochenverlust als der kortikale Knochen ausgesetzt. Daher werden osteoporotische Frakturen hauptsächlich in Knochen mit überwiegend trabekularer Struktur gesehen, wie dem distalen Unterarm, dem proximalen Oberarm, dem Brust- und dem Lendenwirbel sowie dem Oberschenkelhals. Änderungen in der KMD, die auf einen Verlust des trabekularen Knochens bezogen sind, unterschätzen stark den Verlust der Widerstandsfähigkeit gegen Frakturen. Als Folge zeigen KMD-Messungen von Patienten mit osteoporotischen Frakturen eine starke Überlappung mit KMD-Messungen normaler Subjekte der selben Kohorte.
Die obige Information erklärt, warum Personen mit niederenergetischen Frakturen normale KMD-Werte haben können, verglichen mit altersangepaßten normalen Kontrollen. Somit kann der Bedarf bestehen, zwischen der Knochenquantität und der Knochenqualität zu unterscheiden, wobei die Knochenqualität enger auf die mechanische oder biomechanische Festigkeit des Knochens bezogen ist.
Ein Maß der gesamten Knochenfestigkeit kann natürlich aus einer Knochenprobe erhalten werden, die dem potentiellen Patienten entnommen und einem mechanischen Test unterworfen wurde. Dies erfordert jedoch eine Knochenbiopsie, die schmerzhaft ist und ein leichtes Komplikationsrisiko mit sich führt. Um eine bequeme, kostengünstige, schnelle und sichere Überprüfung der sehr großen Gruppe potentieller Patienten (meistens Frauen nach der Menopause) zu erhalten, sollte die Schätzung der Knochenqualität nicht-invasiv, auf sichere und schnelle Weise durchgeführt werden.
Die vorliegende Erfindung bietet ein nicht-invasives Verfahren, das Maße der Gesamtfestigkeit von Knochen liefert und das sicher, bequem für den potentiellen Patienten, schnell und kostengünstig ist.
In "A New Method for Automatic Recognition of the Radiographic Trabekular Pattern", Wil G.M. Geraets et al., Journal of Bone and Mineral Research, Vol. 5, Nr. 3, 1990, Seiten 227-233, ist ein Verfahren zum Erkennen des trabekularen Musters aus Röntgenbildern offenbart. Entsprechend diesem Verfahren werden zwei Rauscharten aus einem Röntgenbild entfernt, das in einen Bildspeicher abgetastet worden ist: ein hochfrequentes Rauschen wird unter Anwendung eines Median-Filters entfernt, und ein niederfrequentes Rauschen wird unter Anwendung einer Operation einer lokalen Mittelung an dem Bild entfernt.
Nach der Rauschverminderung wird eine Segmentierung der Bilddaten durchgeführt, wobei das Bild binärisiert und anschließend ausgewaschen ("eroded") wird (indem eine Median- Achsentransformation angewandt wird), um lediglich zentrale Linien mit einer Dicke von einem Pixel zurückzubehalten. Insgesamt werden sieben Merkmale aus den segmentierten und ausgewaschenen Bildern abgeleitet. Diese Merkmale werden anschließend mit den Messungen der Knochenmineraldichte der Knochen korreliert.
Aufgrund der nahezu unwesentlichen Korrelation zwischen den erhaltenen Bildmerkmalen und den KMD-Messungen in dem oben erwähnten Bezugsdokument, ist es jedoch offensichtlich, daß eine höhere Korrelation und somit eine präzisere Schätzung der Knochenqualität erforderlich ist.
In "computerized Radiographic Analysis of Osteoporosis: Preliminary Evaluation", Philip Caligiuri et al., Radiology 1993, 186; 417-474, ist ein Verfahren analysiert, bei dem Röntgenbilder des Lendenwirbels abgetastet werden. In dieser Referenz wurde das Leistungsspektrum ermittelt und der quadratische Mittelwert (QMW) und das erste Moment (EM) davon mit typischen, an den Knochen durchgeführten KMD- Messungen verglichen. Es wurde sowohl geschlossen, daß die KMD nicht sehr gut mit dem EM und dem QMW korrelierten, als auch, daß dies an dem EM und dem QMW liegen könne, die Zusatzinformation enthalten. In dieser Referenz wurden lediglich das Leistungsspektrum und daraus ableitbare Merkmale untersucht.
Die Medical Physics, v. 21, Nr. 8, Aug. 94, Woodbury, NY (US); Seiten 1293-1300 (S.N.C. Cheng et al.), beschreibt ein Verfahren zum Schätzen der Knochenqualität auf der Basis zweidimensionaler Handröntgenbilder, welche Information betreffend die trabekulare Struktur eines Teils eines Knochens umfassen. Das Verfahren unterwirft die Bilddaten einer statistischen Analyse, die umfaßt:
- eine Prozedur, ein automatischer adaptiver Schwellwertvergleich, bei der niederfrequente Intensitätsvariationen bezüglich Bilddaten vermindert werden, die auf die trabekulare Struktur und die Knochenmasse bezogen sind;
- eine Bildmanipulations- und Merkmalsextraktionsprozedur, bei welcher die lokale Bildinformation und deren lokale Intensitätsvariation verwendet werden, Information zu extrahieren, die auf die trabekulare Struktur des Teils des Knochens bezogen ist; wobei diese umfaßt:
- Ermitteln einer Schätzung der Pixelintensitätsverteilung der Bilddaten;
- Ermitteln einer Schätzung des projizierten Musters der Bilddaten, um einige Information in den Bilddaten hervorzuheben, und Unterziehen der manipulierten Bilddaten einer morphologischen Operation und Extrahieren von Merkmalen;
- Ermitteln, auf der Basis einer Frequenzanalyse der Bilddaten, von Merkmalen, die auf die Periodizität des Unterschieds zwischen normalen und abnormalen Fällen bezogen sind;
- eine Schätzprozedur, bei der die Knochenqualität auf der Basis der abgeleiteten Merkmale und anderer Merkmale geschätzt wird.
Die vorliegende Erfindung schafft ein Verfahren und ein System zum Bereitstellen einer leicht durchführbaren, jedoch hoch zuverlässigen Schätzung der Knochenqualität.
Somit bezieht sich die vorliegende Erfindung auf ein Verfahren gemäß Anspruch 1.
Wie oben angegeben, ist in dem vorliegenden Kontext die "Knochenqualität" nicht gleichzusetzen mit der "Knochenquantität", wie eine Knochenmineraldichte, da selbst ein geringer Verlust der Knochenquantität zu einem großen Verlust der Knochenqualität führen kann, falls der Knochenverlust an kritischen Stellen in der trabekularen Struktur stattgefunden hat. Somit ist in dem vorliegenden Kontext die "Knochenqualität" ein Maß, das eng mit dem Risiko einer Fraktur des Knochens verbunden ist, wie es in Belastungs/Zug-Auswertungen von Knochenbiopsien demonstriert worden ist, bei denen ebenfalls die mikroskopische Struktur ausgewertet wurde (siehe z. B. Lis Mosekilde bezüglich dieses Themas).
Selbst wenn die Knochen und somit die trabekularen Strukturen von Natur aus dreidimensional sind, vermittelt eine Projektion dieser Struktur auf zwei Dimensionen, wie ein Röntgenbild, das ein sog. projiziertes trabekulares Muster ergibt, ausreichend Information über die trabekulare Struktur und somit die Knochenqualität, um eine zufriedenstellende Schätzung der Knochenqualität zu liefern.
Information bezüglich der trabekularen Struktur können lokale Variationen in beispielsweise einer Graustufen- Information oder irgendeine andere Information sein, aus der Merkmale bezüglich der trabekularen Struktur abgeleitet werden können.
Niederfrequente Intensitätsvariationen in den Bilddaten hängen natürlich davon ab, wie die Bilddaten erhalten wurden. Falls die Bilddaten auf der Grundlage von Röntgenbildern erhalten wurden, können die niederfrequenten Intensitätsvariationen von beispielsweise Streuungen der elektromagnetischen Strahlung, von anatomischen Strukturen, welche den belichteten Teil der Trabekulae umgeben, wie ein kortikaler Knochen, Fettgewebe und Muskeln variierender Dicke, sowie z. B. von einer inhomogenen Röntgenbestrahlung des Teils des Knochens hervorgerufen werden.
Gemäß dem ersten Aspekt der vorliegenden Erfindung werden eine Information über die lokale Bildintensität und eine Variation in der lokalen Intensität verwendet, Information bezüglich der trabekularen Struktur des Teils des Knochens zu extrahieren. In beispielsweise digitalisierten Röntgenbildern können die einzelnen Pixelwerte die Information über die lokale Bildintensität sein, während eine Variation der lokalen Intensität auf die Texturinformation bezogen ist, die in z. B. Inhomogenitäten in den Bilddaten enthalten ist.
Die extrahierten Merkmale, die aus der Bildmanipulations- und Merkmalsextraktionsprozedur resultieren, quantifizieren Eigenschaften der trabekularen Struktur und somit der Knochenqualität. Die extrahierten Merkmale werden nachfolgend in eine Schätzprozedur eingeführt, bei der eine vorbestimmte Beziehung zwischen Merkmalen und der Knochenqualität die Schätzprozedur in die Lage versetzen, die Knochenqualität des Vertebraten zu schätzen.
In dem vorliegenden Kontext bedeutet Information, wie Betragsinformation, "hervorzuheben", vorherrschende Frequenzinformation herauszustellen. Dies kann entweder dadurch herbeigeführt werden, daß die vorherrschende Information verstärkt wird, oder daß die weniger dominante Information vermindert wird - optional beides.
Gemäß dem ersten Aspekt der Erfindung wird eine Schätzung der Knochenqualität in einer Schätzprozedur auf der Grundlage einer vorbestimmten Beziehung zwischen erhaltenen Merkmalen und Parametern für die Referenzknochenqualität erhalten. Diese vorbestimmte Beziehung wird typischerweise durch statistische Modellbildung aufgestellt, bei welcher erklärende Variablen (Bildmerkmale und optional andere erklärende Merkmale bezüglich der Knochenqualität) verwendet werden, um entsprechende Daten über die Referenzknochenqualität (Antwortvariable) zu modellieren.
Wie oben beschrieben, sind, selbst wenn in dem vorliegenden Kontext die Knochenqualität nicht gleich der Knochenquantität ist, die Parameter für die Referenzknochenqualität vorzugsweise Parameter, die auf die Festigkeit des Knochens bezogen sind, anstatt z. B. auf die KMD oder den KMG, die sich, wie oben erläutert, auf die Dichte des Knochens beziehen, anstatt auf dessen Festigkeit. Um jedoch eine noch höhere diagnostische Empfindlichkeit zu schaffen, wird eine Kombination der vorliegenden Erfindung mit einer Messung des KMG/KMD als geeignet angesehen. KMG/KMD-Messungen können entweder durch die oben erwähnte Zwei-Wellenlängen- Technik erhalten werden oder sie können, falls die Bilddaten auf der Grundlage eines Röntgenbildes erhalten werden, aus den selben Bilddaten verfügbar sein, bei denen z. B. ein Standard-Aluminiumkeil zusammen mit dem zu untersuchenden Knochen bestrahlt worden ist.
Zur Zeit basiert die am häufigsten anzutreffende Art, nicht-invasive Information über die Knochenquantität und Knochenqualität zu erhalten, auf Röntgenbestrahlung. Rönt genbilder können leicht verwendet werden, Bilddaten für die Verwendung in der vorliegenden Erfindung zu erzeugen. Da bevorzugt wird, die statistische Analyse auf einem Computer durchzuführen, kann das analoge Röntgenbild digitalisiert werden und in einen Computer eingeführt werden, indem das Röntgenbild abgetastet wird.
Natürlich wird die Auflösung des zum Erhalt des Röntgenbildes verwendeten Röntgenfilms eine Auswirkung auf die Qualität der Bilddaten haben. Daher wird z. Zt. bevorzugt, daß die Auflösung des Röntgenfilms wenigstens 4 Linienpaare pro Zentimeter hat, wie wenigstens 5 Linienpaare pro Zentimeter. Selbst wenn dies zum Erhalt einer zufriedenstellenden Schätzung der Knochenqualität ausreichend sein kann, wird bevorzugt, daß die Auflösung des Röntgenfilms wenigstens 10 Linienpaare pro Zentimeter, wie wenigstens 25 Linienpaare pro Zentimeter, vorzugsweise wenigstens 50 Linienpaare pro Zentimeter hat. Es können Auflösungen des Röntgenfilms bis zu wenigstens 100 Linienpaaren pro Zentimeter, wie wenigstens 250 Linienpaaren pro Zentimeter und wahrscheinlich so viel wie 500 Linienpaaren pro Zentimeter, wie wenigstens 600 Linienpaaren pro Zentimeter erhalten werden. Natürlich stellt eine Auflösung dieser bevorzugten Größe sicher, daß in dem Röntgenbild ein großer Informationsgehalt vorliegt.
Natürlich wird bevorzugt, daß das Abtasten der Röntgenbilder mit einer ausreichend großen Auflösung durchgeführt wird, um sicherzustellen, daß nur ein Minimum an relevanter Information des Röntgenbildes bei der Übertragung der Bilddaten verloren geht. Somit wird bevorzugt, daß das Abtasten mit einer Auflösung von wenigstens 10 Linien pro cm, wie wenigstens 25 Linien pro cm, vorzugsweise wenigstens 100 Linien pro cm, wie wenigstens 200 Linien pro cm, wie wenigstens 250 Linien pro cm, durchgeführt worden ist.
Ferner wird bevorzugt, daß die Auflösung des Scanners besser als 4 echte Bits pro Pixel, wie besser als 6 echte Bits pro Pixel, vorzugsweise besser als 8 Bits pro Pixel, ist.
Damit der Scanner das Bild wirklich abtasten kann, sollte er natürlich in der Lage sein, die Röntgenbilder zu durchleuchten.
Es wurde herausgefunden, daß der ESKOSCAN 2450 von Eskofot A/S die obigen Kriterien erfüllt und für die Verwendung in dem Verfahren der Erfindung bestens geeignet ist.
Es wird z. Zt. bevorzugt, daß die Hintergrundkorrekturprozedur wenigstens das Vermindern oder optional das Entfernen niederfrequenter Information mit einer Frequenz, die wesentlich niedriger als der Abstand der projizierten Trabekulae ist, umfaßt. Wie oben beschrieben, wird der niederfrequente Spektralgehalt der Bilddaten typischerweise von dem kortikalen Knochen, Fettgewebe und Muskeln variierender Dicke, sowie z. B. inhomogener Röntgenbestrahlung des Teils des Knochens, verursacht. Diese Art unerwünschten Effekts ist natürlich nicht vermeidbar, wenn nicht-invasive Bilddaten auf der Grundlage von Röntgenbildern erhalten werden. Andere nicht-invasive Bilderfassungstechniken, wie beispielsweise eine MR- und CT-Bildwiedergabe, können jedoch Bilddaten erzeugen, die nicht automatisch oder nicht in dem selben Ausmaß an dieser Art unerwünschten Effekts "leiden".
Selbst wenn es schwierig ist, unabhängig von der Bildauflösung exakt die Grenze zwischen unerwünschter niederfrequenter Information und der gewünschten höherfrequenten relevanten Information zu quantifizieren, wird z. Zt. bevorzugt, daß Information mit Frequenzen, die die Hälfte oder weniger der Beabstandung der projizierten Trabekulae sind, wenigstens vermindert oder optional entfernt wird. Noch bevorzugter wird Information mit Frequenzen, die ein Viertel oder weniger, wie ein Zehntel oder weniger, der Beabstandung der projizierten Trabekulae sind, wenigstens vermindert oder optional entfernt.
Eine bevorzugte Art für die Verminderung oder vorzugsweise für die Entfernung der niederfrequenten Information ist die Anwendung einer Hintergrundkorrekturprozedur, die das Erzeugen sekundärer Bilddaten als ein Ergebnis eines Durchführens eines Median-Filterns mit einer vorbestimmten Kerngröße und Subtrahieren dieses Ergebnisses von den ursprünglichen Bilddaten umfaßt. Einer der Vorteile der Anwendung eines Median-Filters liegt darin, daß diese Operation kantenerhaltend ist.
Eine andere Art für die Verminderung oder vorzugsweise für die Entfernung der niederfrequenten Information ist die Anwendung einer Hintergrundkorrekturprozedur, welche die Erzeugung sekundärer Bilddaten als ein Ergebnis des Durchführens eines Mittelwertfilterns mit einer vorbestimmten Kerngröße und Subtrahieren dieses Ergebnisses von den ursprünglichen Bilddaten umfaßt. Ein Mittelwertfiltern ist typischerweise viel schneller als das Median-Filtern. Das Mittelwertfiltern ist jedoch nicht kantenerhaltend. In diesem Teil des Prozesses erzeugte Kanten können einen nachteiligen Effekt haben und später in dem Prozeß falsche oder fehlerhafte Information erzeugen, was zu falschen Schätzungen der Knochenqualität führt.
Es wird typischerweise bevorzugt, daß die Kerngröße höchstens 1/2 der Bilddaten, wie höchstes 1/4 der Bilddaten, vorzugsweise höchstens 1/10 der Bilddaten, noch bevorzugter höchstens 1/20 der Bilddaten, ist.
Eine dritte Art für die Verminderung oder vorzugsweise für die Entfernung der niederfrequenten Information liegt in der Anwendung einer Hintergrundkorrekturprozedur, welche das globale Anpassen eines zweidimensionalen Polynoms an die Bilddaten und das Erzeugen hintergrundkorrigierter Bilddaten auf der Grundlage der Residuen bzw. Reste der Anpaßprozedur umfaßt. Abgesehen von potentiellen Schwierigkeiten beim Bestimmen der optimalen Ordnung des Polynoms, die der Durchschnittsfachmann kennt, bietet dieses Verfahren potentiell eine extrem schnelle Hintergrundkorrektur der Bilddaten.
Es wird erwogen, daß eine geeignete Ordnung des Polynoms höchstens 15, wie höchstens 10, noch bevorzugter höchstens 5, ist.
Nach Durchführung einer Hintergrundkorrektur der Bilddaten wird eine Bildmanipulations- und Merkmalsextraktionsprozedur durchgeführt, um Merkmale zu erhalten, welche Textureigenschaften der trabekularen Struktur quantifizieren.
Eine wichtige, sich auf die Knochenqualität beziehende Textureigenschaft ist der Grad der Anisotropie des projizierten trabekularen Musters. Diese Eigenschaft kann in einer intuitiv greifbaren Weise beschrieben werden, indem das Leistungsspektrum der Bilddaten verwendet wird. Das Leistungsspektrum kann erhalten werden, indem direkte Methoden, Auto-Kovarianzmethoden oder parametrische Methoden angewandt werden.
Eine klassische spektrale Schätzung jedoch, die direkte Methoden und Auto-Kovarianzmethoden umfaßt, liefert typischerweise relativ nicht-konsistente spektrale Schätzungen mit einer geringen Auflösung. Alternativ können Glättungsmethoden angewandt werden. Dies führt jedoch zu den bekannten Varianz/Auflösungsnachteilen.
Daher wird erfindungsgemäß eine parametrische spektrale Schätzung angewandt, da dieser Schätzer ein konsistenter spektraler Schätzer mit ausgezeichneten Auflösungseigenschaften ist.
Eine Art für den Erhalt der klassischen Schätzung des Leistungsspektrums ist, die Bilddaten, optional mit einem Fenster gewichtet, einer Fast-Fourier-Transformation zu unterziehen. Dies ist eine sog. direkte Methode, von der es verschiedene Varianten gibt, wobei die wahrscheinlich populärste die sog. Welch-Methode ist.
Eine andere Art für den Erhalt der klassischen Schätzung des Leistungsspektrums ist, eine Schätzung der Auto- Kovarianzfunktion der Bilddaten, optional gewichtet mit einem Fenster, einer Fast-Fourier-Transformation zu unterziehen.
Wie oben erwähnt, bieten die parametrischen Methoden potentiell eine höhere Auflösung als die klassischen Methoden. Das Grundprinzip parametrischer Methoden liegt darin, ein geeignetes Modell der Bilddaten zu identifizieren, von denen angenommen wird, daß sie eine homogene Textur sind, und nachfolgend die Parameter des Modells zu schätzen. Nachdem diese Parameter geschätzt worden sind, kann die parametrische Schätzung des Leistungsspektrums des angepaßten Models ermittelt werden. In einer Dimension hat sich herausgestellt, daß dieser Typ einer spektralen Schätzung besser ist als die klassischen spektralen Schätzungen.
Die Methoden, auf deren Basis die parametrische Schätzung erhalten wird, werden z. Zt. vorzugsweise aus der Gruppe ausgewählt, die besteht aus: kausalen simultanen autoregressiven Modellen mit gleitendem Mittelwert (SARMA = "Simultaneous Auto-Regressive Moving Average", d. h. simultane auto-regressive mit gleitendem Mittelwert) und nichtkausalen SARMA-Modellen.
Wenn räumliche Daten betrachtet werden, wie Bilddaten, ist eine Kausalität nicht natürlich. Das Aufgeben einer Kausalität führt jedoch zu ernsthaften Problemen bei der Schätzung der Parameter, da der Schätzer gemäß der kleinsten Quadrate nicht länger konsistent ist. Daher wird oftmals bevorzugt, den Bilddaten eine künstliche Richtungsabhängigkeit zu überlagern. Daher kann mit einem kausalen Modell das Prinzip der Schätzung gemäß den kleinsten Quadraten angewandt werden.
Ein Nachteil, den Bilddaten eine künstliche Richtungsabhängigkeit zu überlagern, ist jedoch der sog. "richtungsabhän gige Bias", der z. B. bei der Schätzung der parametrischen Schätzung des Leistungsspektrums ziemlich deutlich sichtbar werden kann. Folglich ist die Anwendung kausaler Modelle in dem vorliegenden Kontext keine selbstverständliche Wahl. Dieses Problem wird teilweise gelöst, indem der "Parallel- Widerstand" ("parallel-resistor") gemittelte spektrale Schätzer angewandt wird, der von Chien ("Two-dimensional Spectral estimation from Autoregressive Models with varying Areas of Support" von 1981) eingeführt worden ist.
Zur Zeit wird bevorzugt, daß das Verfahren, auf dessen Grundlage das parametrische Spektrum erhalten wird, ein nicht-kausales SARMA-Modell ist, das den sog. MORSE- Schätzer anwendet. Dieser Schätzer wird weiter unten beschrieben.
Nach Erhalt der parametrischen Schätzung des Leistungsspektrums der Bilddaten werden ein oder mehrere Merkmale extrahiert, welche Eigenschaften der Knochenqualität quantifizieren.
Die Verteilung der Energie in der parametrischen Schätzung des Leistungsspektrums bezieht sich direkt auf die Anisotropie der Bilddaten. Somit können mehrere Typen an Merkmalen vorgeschlagen werden, welche diese Anisotropie quantifizieren.
Zur Zeit wird bevorzugt, wenigstens ein Merkmal zu ermitteln, das auf Parameter einer Kontur bezogen ist, welche einen vorbestimmten Prozentsatz der Energie der parametrischen Schätzung des Leistungsspektrums umgibt.
Diese Kontur kann bestimmt werden durch
(a) Definieren einer Kontur in Ausdrücken einer Achse in jeder Dimension um das Zentrum der parametrischen Schätzung des Leistungsspektrums, wobei alle Punkte auf der Kontur im wesentlichen den selben Abstand von dem Zentrum der parametrischen Schätzung des Leistungsspektrums haben, und die Kontur weniger als den vorbestimmten Prozentsatz der Energie der parametrischen Schätzung des Leistungsspektrums umgibt,
(b) für jede Dimension der Daten, Berechnen des Prozentsatzes der Energie der parametrischen Schätzung des Leistungsspektrums, das von einer gedehnten Kontur umgeben ist, in der die Achse in der in Frage kommenden Dimension um einen vorbestimmten Abstand verlängert ist,
(c) Verlängern der Achse der Kontur in der Dimension, in der die Vergrößerung der von der gedehnten Kontur umgebenen Energie am größten ist, um den vorbestimmten Abstand, und
(d) Wiederholen der Schritte (b) und (c), bis der von der Kontur umgebene Prozentsatz den vorbestimmten Prozentsatz überschreitet oder angleicht.
Die Energieverteilung in der parametrischen Schätzung des Leistungsspektrums ist typischerweise länglich. Tests haben gezeigt, daß sich die Verteilung der Energie in der parametrischen Schätzung der Leistungsspektren osteoporotischer Patienten wesentlich von denjenigen nicht-osteoporotischer Personen unterscheidet. Dieser Unterschied kann beispielsweise in einem Unterschied in der Fläche, die von der oben erwähnten Kontur bedeckt ist, oder in der Länglichkeit der oben erwähnten Kontur gesehen werden.
Eine Art für den Erhalt geeigneter Merkmale aus der Kontur ist, die Form der Kontur zu definieren und nachfolgend Merkmale aus Parametern der Kontur abzuleiten.
Ein geeigneter Weg für die Definition der Kontur ist, die Achsen der Kontur derart zu definieren, daß diese orthogonal und entlang von Hauptrichtungen der parametrischen Schätzung des Leistungsspektrums ausgerichtet sind.
Eine natürliche Form dieser Kontur kann eine ellipsoidale Kontur sein und ein oder mehrere Merkmale können aus deren Halbachsen abgeleitet werden. Es kann jedoch auch eine rechteckige Kontur verwendet werden.
Das oben erwähnte Verfahren für die Definition der Kontur wird vorzugsweise so durchgeführt, daß die kleinstmögliche Kontur der gewünschten Form erhalten wird. Dies stellt sicher, daß die Kontur eindeutig definiert ist und daß sich die aus Parametern der Kontur abgeleiteten Merkmale so gut wie möglich auf die Anisotropie der Bilddaten beziehen.
Ein alternativer Weg für die Extraktion von Information in dem Frequenzbereich ist, Information aus beispielsweise der Höhe, Breite oder Gesamtfläche unter Peaks eines geglätteten Periodogramms abzuleiten. Aus dieser Information können sich auf die Periodizität der Trabekulae beziehende Merkmale abgeleitet werden. Auf diese Weise wird Information bezüglich der Periodizität der Trabekulae erhalten, indem eine Frequenzanalyse der Bilddaten durchgeführt wird.
Es können mehrere andere Wege für die Extraktion von Merkmalen aus der parametrischen Schätzung des Leistungsspektrums angewandt werden, wie das Anpassen einer Gauß'schen Kurve an die normalisierte parametrische Schätzung des Leistungsspektrums und das Ableiten von Merkmalen aus dieser Gauß'schen Kurve. Andere Verfahren verwenden Ringe und Keile, wie von Weszka et al. (1976) beschrieben.
Ein Bildmanipulations- und Merkmalsextraktionsverfahren, das alternativ oder ergänzend zu der obigen Anwendung der parametrischen Schätzung des Leistungsspektrums ist, liegt darin, das projizierte trabekulare Muster der Bilddaten unter Anwendung eines Fourier-Verfahrens zu schätzen, um die Information bezüglich der trabekularen Struktur in den Bilddaten hervorzuheben, und die manipulierten Bilddaten einer morphologischen Operation zu unterziehen, und Merkma le bezüglich der trabekularen Struktur aus dem geschätzten projizierten trabekularen Muster zu extrahieren.
Die Hervorhebungsprozedur ist wesentlich, um das projizierte trabekulare Muster robust herauszubringen. An dieser Stelle haben sich Fourier-Verfahren als geeignet erwiesen, und es wird bevorzugt, daß das Fourier-Verfahren umfaßt
- Unterziehen der Bilddaten einer Fourier-Transformation,
- Durchführen einer nachfolgenden mathematischen Transformation der Fourier-transformierten Daten,
- Umwandeln der transformierten Daten zurück in den Ortsbereich.
Natürlich kann die nachfolgende mathematische Transformation linear oder nicht-linear sein.
Selbst wenn bevorzugt wird, daß die nachfolgende mathematische Transformation durchgeführt wird, indem der Betrag der Fourier-transformierten Daten auf eine Potenz größer als 1, wie eine Potenz im Bereich von 1,1-10, vorzugsweise 1,5-4, wie ungefähr 2,0, erhöht wird, können auch andere Transformationen, wie andere Filterarten, durchgeführt werden.
Um die Information bezüglich der trabekularen Struktur weiter hervorzuheben, kann eine dritte mathematische Transformation durchgeführt werden, um lediglich einen Teil der Betragsinformation der transformierten Daten beizubehalten. Dies kann beispielsweise durch Entfernen wenigstens eines Teils der Information in dem Betrag des Spektrums geschehen.
Nachdem die Information bezüglich der trabekularen Struktur in den Bilddaten hervorgehoben worden ist, werden die resultierenden Bilddaten für die Vorbereitung auf die Merkmalsextraktion vorzugsweise morphologischen Graustufen- Operationen ausgesetzt.
Zur Zeit liegt der Zweck der morphologischen Graustufen- Operationen darin, eine binäre Darstellung des projizierten trabekularen Muster (PTM) zu extrahieren. Die Vordergrundfarbe stellt die Trabekulae dar, während die Hintergrundfarbe die Hohlräume darstellt.
Das obige PTM wird vorzugsweise erhalten, indem die sog. morphologische Top-Hat-Operation angewandt wird, gefolgt von einem Schwellwertvergleich. Das Ergebnis dieser Operationen wird vorzugsweise weiter "gesäubert", indem kleine isolierte Pixelansammlungen entfernt werden.
Auf der Grundlage des erhaltenen PTM kann eine Anzahl verschiedener Merkmale ermittelt werden. Nachfolgend wird die Erzeugung lediglich einiger einer immensen Anzahl möglicher Merkmale beschrieben.
Ein Verfahren, Merkmale aus dem PTM zu erzeugen, ist das Durchführen einer Operation (der sog. "Abstandstransformation"), die für jeden Hintergrundpixel die Bestimmung des Abstands gemäß einer gegebenen Metrik von dem Hintergrundpixel zu dem nächsten Vordergrundpixel umfaßt. Die aus dieser Operation extrahierten Merkmale sind tatsächlich intuitiv greifbar, da sie sich direkt auf den zwischentrabekularen Abstand beziehen.
Die obige Metrik ist z. Zt. der euklidische Abstand. Andere Metriken können sich jedoch in dem vorliegenden Kontext evtl. als geeigneter erweisen.
Auf der Abstandstransformation basierende Merkmale können beispielsweise aus einem Mittelwert und/oder einer Standardabweichung und/oder dem Variationskoeffizienten und/oder dem Schiefe-Koeffizienten und/oder dem Kurtosis-Koeffizienten der bestimmten Abstände abgeleitet werden.
Ein weiteres Verfahren zum Erzeugen von Merkmalen aus dem PTM ist das Durchführen einer Operation (die nachfolgend als "Sternflächentransformation" bezeichnet wird), die für jeden Hintergrundpixel das Bestimmen eines Maßes umfaßt, basierend auf einer Bestimmung des Abstandes von dem Hintergrundpixel zu dem nächsten Vordergrundpixel, in einer Anzahl vorgegebener Richtungen in den Bilddaten.
Wiederum können sich aus der Sternflächentransformation ableitbare Merkmale beispielsweise auf einen Mittelwert und/oder eine Standardabweichung und/oder den Variationskoeffizienten und/oder den Schiefe-Koeffizienten und/oder den Kurtosis-Koeffizienten, optional den maximalen Abstand, der bestimmten Maße beziehen.
Zusätzlich zu den aus den Bilddaten ableitbaren, oben erwähnten Merkmalen kann bevorzugt werden, zusätzliche Daten in die Schätzprozedur einzugeben. Es ist selbstverständlich, daß eine große Anzahl an Merkmalen, die sich auf den bestrahlten Knochen beziehen, die jedoch nicht aus den Bilddaten ableitbar sind, die Genauigkeit der Schätzung der Knochenqualität verbessern können.
Da das vorliegende Verfahren verwendet werden kann, die Knochenqualität irgendeines Knochens von irgendeinem Vertebraten zu schätzen, kann zusätzliche Information, wie das Alter und/oder Geschlecht und/oder Spezies und/oder Rasse und/oder der spezifische, in dem Vertebraten betrachtete Knochen und/oder eine geschätzte Knochenmineraldichte und/oder ein geschätzter Knochenmineralgehalt in die Schätzprozedur einbezogen werden.
Selbst wenn die KMD in die Schätzprozedur eingeführt werden kann, kann dieses Maß ebenfalls optional aus den Bilddaten bestimmt werden. Die KMD kann bestimmt werden, indem Daten aus einem Referenzobjekt bei der Belichtung des Knochens mit der elektromagnetischen Strahlung und auf der Grundlage der Absorption der elektromagnetischen Strahlung des Knochens und des Referenzobjektes einbezogen werden.
Der Zweck der Schätzprozedur liegt darin, vorgegebene biomechanische Eigenschaften des Knochens auf der Grundlage der aus den Bilddaten extrahierten, eingeführten Merkmale und anderer erklärenden Variablen zu schätzen.
Ein Weg für den Erhalt einer Schätzprozedur vom obigen Typ ist, die Schätzprozedur auf ein statistisches Modell zu basieren, unter Berücksichtigung der Korrelationsstruktur in dem Datensatz, um den wesentlichen Merkmalen entsprechend der vorbestimmten Beziehung geeignete Gewichtungen zuzuordnen.
Das oben erwähnte Modell kann auf verschiedene Weise bestimmt werden. Das Aufstellen eines Modells von diesem Typ erfordert jedoch den Erhalt entsprechender Werte aller relevanten Merkmale und einer Antwortvariabel bezüglich der biomechanischen Eigenschaft, die durch die Schätzprozedur geschätzt werden soll. Die in Frage kommende biomechanische Eigenschaft kann ein absolutes oder ein relatives Maß der Knochenqualität sein.
Beispiele von in Frage kommenden biomechanischen Eigenschaften des Knochens sind die mechanische Knochenfestigkeit und/oder eine Messung der Knochenmineraldichte und/oder eine Messung des Knochenmineralgehalts und/oder eine Bewertung durch einen erfahrenen Radiologen. Natürlich wird das Verfahren der vorliegenden Erfindung einen Wert ausgeben, der der Kalibrierung des Verfahrens entspricht: falls das Verfahren in Richtung eines Festigkeitsparameters kalibriert ist, wird sich die Ausgabe des Verfahrens auf die Festigkeit des in Frage kommenden Knochens beziehen.
Ein wichtiger Zweck der vorliegenden Erfindung liegt darin, ein Verfahren zum Bewerten der Knochenfestigkeit und des Frakturrisikos mit signifikantem klinischem Wert zu schaffen, also eine bessere Schätzung des Frakturrisikos zu schaffen als die besten heutzutage verfügbaren Techniken, die Messungen des Knochenmineralgehalts über die Doppel- Röntgenabsorptionsmessung (DXA-BMC-Messungen = Double X-ray Absorptiometry Bone Mineral Content Measurements).
Da die Kalibrierung der Schätzprozedur die geschätzten Parameter des Knochens bestimmt, sollten diese Parameter derart ausgewählt werden, daß sie mit der Knochenqualität korrelieren. Die Beurteilung des Ausmaßes, mit dem die Ausführungsbeispiele der vorliegenden Erfindung die Knochenqualität schätzen können, kann durchgeführt werden durch:
1. Kohorten-Studien, bei denen das Auftreten neuer Frakturen aufgezeichnet und auf die anfängliche Schätzung bezogen wird,
2. Querschnittsstudien, bei denen Frakturfälle und alters- und geschlechtsangepaßte Kontrollen mit den von dem erfindungsgemäßen Verfahren erhaltenen Schätzungen und optional ferner mit DXA-BMC-Messungen verglichen werden, oder
3. Anwendung biomechanischer Tests von Knochenproben, denen eine Abschätzung der Knochenqualität gemäß der Erfindung und optional zusätzlich gemäß DXA-BMX-Messungen des in Frage kommenden Knochens vorangegangen sind.
Durch Verwendung entweder der distalen Speiche, des Lendenwirbels, des Oberschenkelhalses oder irgendeines Knochens, in dem zweidimensionale Bilddaten erhalten werden können, die Information bezüglich der trabekularen Struktur des Knochens umfassen, können kohärente Werte für die bestimmte Knochenfestigkeit, DXA-BMC und die erfindungsgemäße Schätzung auf der Grundlage menschlicher post-mortem Knochenproben oder Proben von anderen Vertebraten durchgeführt werden.
Die Kalibrierung des Verfahrens der vorliegenden Erfindung kann ebenfalls entsprechend der oben dargestellten Verfahren kalibriert werden, da diese Verfahren eine Information und Schätzung über die Knochenqualität auf der Grundlage der entsprechenden Bilddaten erzeugen.
Die Knochenfestigkeit kann z. B. bewertet werden, indem die Plastizität und die Maximallast eines Knochens in einem Belastungs-Zug-Diagramm direkt gemessen werden (siehe z. B. Lis Mosekilde). Diese Parameter werden vorzugsweise in sowohl der Richtung der vertikalen Trabekulae als auch in der orthogonalen Richtung erhalten, um die vollständigste Schätzung der Knochenqualität zu haben.
In dem Verfahren der Erfindung können in der Schätzprozedur im wesentlichen alle Modelltypen angewandt werden. Zur Zeit wird ein bevorzugtes Modell aus der Gruppe ausgewählt, die besteht aus: einem allgemeinen linearen Modell, einem verallgemeinerten linearen Modell, einem künstlichen neuronalen Netzwerk, einem kausalen Wahrscheinlichkeitsnetz oder Klassifikations- und Regressions-Bäumen.
Wie oben beschrieben, kann das vorliegende Verfahren verwendet werden, die Knochenqualität des Knochens in irgendeinem Vertebrat zu schätzen. Somit kann der Vertebrat ein Mensch, ein Pferd, ein großer Affe, ein Menschenaffe, ein Schwein, eine Kuh etc. sein und der aktuelle Knochen kann natürlich praktisch irgendein Knochen in dem Vertebraten sein, wie der Radius, der Femur, der Corpus vertebrae, der Calcaneus, der Talus, der Os carpi, der Metatarsi, der Metacarpi, der Falanges, der Tibia, der Fibula, der Patella, der Ulna, der Humerus, der Mandible, der Clavicula, der Scapula, der Os coxae, der Os naviculare, der Os cuboideum, der Os cuneiform I, der Os cuneiform II oder der Os cuneiform III.
Die Knochenqualität des Vertebraten wird vorzugsweise mit einem Mehrfach-Korrelationskoeffizienten besser als 0,5 oder 0,55, wie besser als 0,6, geschätzt. Natürlich ist die Korrelation zwischen der tatsächlichen Knochenqualität und der geschätzten Knochenqualität um so höher, je höher der Mehrfach-Korrelationskoeffizient ist. Daher wird ein Mehrfach-Korrelationskoeffizient besser als 0,65, wie besser als 0,7 oder besser als 0,8, wie besser als 0,85, bevorzugt, um eine Korrelation mit dem tatsächlichen Wert zu erhalten, die so hoch wie möglich ist.
Zum Schätzen nicht-kausaler simultaner autoregressiver Modelle mit gleitendem Mittelwert in zwei oder mehreren Dimensionen ist ebenfalls ein Verfahren geeignet, wobei das Verfahren die Optimierung eines gegebenen Direktmaßes der Flachheit des Restspektrums des Modells umfaßt. In dem vorliegenden Kontext bedeutet eine "Optimierung eines gegebenen Direktmaßes der Flachheit", daß die Flachheit maximiert wird.
Es werden typischerweise nicht-lineare Optimierungsprozeduren angewandt, um Maße vom obigen Typ zu maximieren. Eine geeignete nicht-lineare Prozedur kann die folgenden Schritte umfassen
(a) Erzeugen eines Satzes Anfangsparameter für das Modell,
(b) Erzeugen des Restspektrums des Modells auf der Grundlage der Parameter,
(c) Ermitteln des Maßes der Flachheit des Restspektrums,
(d) Ermitteln einer neuen Iteration der Parameter auf der Grundlage des Maßes der Flachheit und einer Suchrichtung in dem Parameterraum,
(e) Wiederholen der Schritte (b)-(d), bis ein vorgegebenes Abbruchkriterium erreicht ist.
Es können verschiedene Maße der Flachheit der Restparameterschätzung des Leistungsspektrums ausgewählt werden. Ein bekanntes Maß der Flachheit von Wahrscheinlichkeitsverteilungen ist die sog. Entropie. Daher wird in dem vorliegenden Kontext die Entropie auf der normalisierten Restparameterschätzung des Leistungsspektrums gemessen, wobei der Wert bei der Frequenz null unbeachtet bleibt.
Bevorzugte Ausführungsbeispiele von Aspekten der Erfindung werden nunmehr mit Bezug auf die Zeichnung beschrieben, in der
Fig. 1 ein typisches Röntgenbild eines Handgelenks zeigt, das von einem ESKOFOT-Scanner erhalten wurde (siehe Beispiel 1),
Fig. 2 ein extrahiertes Unterbild des Bildes aus Fig. 1 zeigt,
Fig. 3 das Bild der Fig. 2 zeigt, nachdem es mit einer 25 · 25-Kerngröße median-gefiltert wurde,
Fig. 4 das letzte hintergrundkorrigierte Bild zeigt,
Fig. 5 das Leistungsspektrum des hintergrundkorrigierten Bildes der Fig. 4 zeigt,
Fig. 6 das wiederhergestellte Ergebnis des Bildes aus Fig. 4 zeigt, nachdem es einer Hervorhebungsprozedur unterworfen worden ist,
Fig. 7 das Ergebnis einer morphologischen Top-Hat-Operation auf dem wiederhergestellten Bild der Fig. 6 zeigt,
Fig. 8 das Bild der Fig. 7 zeigt, bei dem kleine isolierte Gruppen an Pixeln (Rauschen) entfernt worden sind (dieses Bild wird als das projizierte trabekulare Muster (PTM) bezeichnet),
Fig. 9 das Ergebnis einer Abstandstransformation zeigt, die für das PTM der Fig. 8 berechnet worden ist,
Fig. 10 die Sternflächentransformation für einen einzelnen Pixel darstellt,
Fig. 11 das Ergebnis der Sternflächentransformation aller Hintergrundpixel in dem PTM der Fig. 8 zeigt,
Fig. 12 das Ergebnis der Maximalabstandstransformation des PTM der Fig. 8 zeigt,
Fig. 13 die kleinstmögliche Ellipse zeigt, die 70% der spektralen Energie für einen 'normalen' Patienten enthält,
Fig. 14 die kleinstmögliche Ellipse zeigt, die 70% der spektralen Energie für einen osteoporotischen Patienten enthält,
Fig. 15 den systematischen Bias der LS-Schätzungen für simulierte isotropische Texturen darstellt (die vertikalen Linien entlang der Kurve skizzieren die empirischen Konfidenzintervalle für den Mittelwert der geschätzten Parameter),
Fig. 16 den systematischen Bias der approximierten ML- Schätzungen darstellt (nahe der Nicht-Stationarität) für simulierte isotropische Texturen (die vertikalen Linien entlang der Kurve skizzieren die empirischen Konfidenzintervalle für den Mittelwert der geschätzten Parameter),
Fig. 17 einen Teil der Likelihood-Funktion für ein bilaterales SAR-Modell erster Ordnung darstellt (unter Vorgabe einer Realisierung des Prozesses mit den Parametern α = β = -0,22735. Die Likelihood-Funktion ist in dem Intervall (α,β) [-0,275, -0,20] · [-0,275, -0,20] gezeichnet),
Fig. 18 darstellt, daß gute Schätzungen erhalten werden können, indem der Torus-ML-Schätzer mit den genäherten ML- Schätzungen als Anfangswerte, kombiniert mit einer Stationaritätsprüfung verwendet wird (falls ein nicht-stationäres Lösungsmodell gefunden wird, wird der Optimierer mit neuen Anfangswerten neu gestartet; die vertikalen Linien entlang der Kurve skizzieren das empirische Konfidenzintervall für den Mittelwert der geschätzten Parameter),
Fig. 19 darstellt, daß gute Schätzungen erhalten werden können, indem der MORSE-Schätzer angewandt wird (die vertikalen Linien entlang der Kurve sind die empirischen 95%- Konfidenzintervalle für den Mittelwert der geschätzten Parameter),
Fig. 20 das Zielfunktional des MORSE-Schätzers nahe der Nicht-Stationarität zeigt (im Gegensatz zu beispielsweise dem ML-Schätzer ist dieses Zielfunktional nicht wellig, was die Durchführung der nicht-linearen Optimierung sehr viel leichter macht),
Fig. 21 den Periodogramm-Schätzer des Leistungsspektrums für ein normales Individuum (links) und für ein osteoporotisches Individuum (rechts) darstellt. Das Spektrum ist invers dargestellt, so daß die Werte höherer Intensität dunkler als die Werte niedrigerer Intensität sind,
Fig. 22 die parametrische NSHP-Schätzung des Leistungsspektrums für ein normales Individuum (links) und für ein osteoporotisches Individuum (rechts) darstellt. Das Spektrum ist invers dargestellt, so daß die Werte höherer Intensität dunkler sind als die Werte niedrigerer Intensität,
Fig. 23 die parametrische WQP-Schätzung des Leistungsspektrums für ein normales Individuum (links) und für ein osteoporotisches Individuum (rechts) darstellt. Das Spektrum ist invers dargestellt, so daß die Werte hoher Intensität dunkler sind als die Werte niedrigerer Intensität,
Fig. 24 ein Experiment zusammenfaßt, das durchgeführt wurde, um die Anwendbarkeit des MORSE-Schätzers zu beschreiben,
Fig. 25 Histogramme darstellt, basierend auf Parametern, die für jede von 1000 künstlich erzeugten 128 · 128 Texturen geschätzt wurden,
Fig. 26 optimale Festigkeiten darstellt, die gegen die vorhergesagten Festigkeiten gezeichnet sind, wobei in dem Modell das Alter nicht einbezogen ist,
Fig. 27 optimale Festigkeiten darstellt, die gegen die vorhergesagten Festigkeiten gezeichnet sind, wobei in dem Modell das Alter einbezogen ist, und
Fig. 28 die Korrelation der logarithmischen optimalen Bruchbelastung und einem einzelnen geeignet gewählten Merkmal darstellt.

BEISPIEL 1: Texturanalyse der Knochenstruktur

In dem nachfolgenden Beispiel wird ein bevorzugtes Ausführungsbeispiel des Verfahrens der Erfindung beschrieben, bei dem Röntgenbilder der distalen Speiche und eines Wirbelkörpers von dem Lendenwirbel (L3) unter Anwendung einer Bildverarbeitung wiederhergestellt werden. Die wiederhergestellten Bilder werden als eine Grundlage zum Extrahieren von Texturmerkmalen verwendet, für die sich gezeigt hat, daß sie mit der Struktur, der Dichte und der Bruchbelastung des trabekularen Knochens gut korrelieren. Es werden mehrere Wege zum Extrahieren vorherrschender Texturmerkmale beschrieben. Ein auf Texturmerkmalen (und ggf. Alter und Geschlecht) basierendes Modell wird beschrieben und es wird gezeigt, daß es gut zwischen osteoporotischen und nichtosteoporotischen Fällen unterscheiden kann.

Ermitteln geeigneter Bilddaten

Die in diesem Beispiel verwendeten Röntgenbilder wurden unter Verwendung eines standardgemäßen Röntgengerätes mit einem Fokusbereich von 0,6 mm², einer Röhrenspannung von 45 kV, einem Fokus-zu-Folien-Abstand von 100 cm, einer einzelnen dünnen Röntgenfolie (es sei bemerkt, daß keine Folie erforderlich ist), unter Verwendung eines Doppelemulsionsfilms (es sei bemerkt, daß ein Einzelemulsionsfilm gleichfalls anwendbar ist) und unter Einstellung des MAS- Produktes erhalten, um einen geeignet belichteten Film zu erhalten.
Die Röntgenbilder wurden in einen Computer unter Verwendung eines ESKOFOT ESKOSCAN 2450 Scanners abgetastet, der bis zu 5000 Linien/cm abtasten kann. Vorliegend wurden jedoch 600 Linien/cm verwendet. Dieser Scanner wurde für die vorliegende Verwendung als höchst geeignet befunden. Ein typisches Bild ist in Fig. 1 gezeigt.
Die Anforderungen an einen für diese Zwecke geeigneten Scanner sind, daß er natürlich den Röntgenfilm durchleuchten kann, daß er vorzugsweise eine Auflösung besser als die Körner in der Film/Folien-Kombination hat (damit keine Information bei der Abtast-Prozedur verloren geht) und daß er vorzugsweise Grauwerte mit wenigstens 8 echten Bits erzeugen sollte, um einen geeigneten dynamischen Bereich der Bilddaten zu erhalten.

Ortstexturmodelle

Vorverarbeitung

Bei dieser Studie wird lediglich die Texturinformation berücksichtigt. Daher wird ein Unterbild in einer Region des Knochens extrahiert, das lediglich einen trabekularen Knochen enthält, wie beispielsweise in Fig. 2 gezeigt. Niederfrequente Intensitätsvariationen aufgrund von beispielsweise Streuung der Strahlung, anatomischen Strukturen, wie dem kortikalen Knochen, Muskeln, Fettgewebe variierender Dicke etc., werden durch Subtrahieren einer median-gefilterten Version des Bildes (in Fig. 3 gezeigt) von dem Bild selbst entfernt (eine Technik, die als "unscharfes Filtern" bekannt ist). In diesem Fall wird ein 31 · 31 Median-Filter verwendet. Das hintergrundkorrigierte Bild ist in Fig. 4 gezeigt. Dieses Bild ist sehr verrauscht und sollte vor weiteren Transformationen und nachfolgender Merkmalsextraktion weiterverarbeitet werden.
Das Rauschen in dem Bild ist deutlich sehr hochfrequent, während die Information, die vorzugsweise betrachtet wird (die trabekulare Struktur), eher niederfrequenter Natur ist.
Das Leistungsspektrum des extrahierten, hintergrundkorrigierten Unterbildes ist in Fig. 5 gezeigt. Nachfolgend wird eine Hervorhebungsprozedur beschrieben, die nicht allgemein für das Entfernen eines hochfrequenten Rauschens anwendbar ist, jedoch gewöhnlich die dominierenden (in diesem Fall niederfrequenten) Merkmale in den Bildern hervorhebt, während sie weniger dominante (in diesem Fall hochfrequentes Rauschen) Merkmale unterdrückt. Diese Eigenschaft macht sie für diesen speziellen Zweck attraktiv.
Algorithmus 1 Die folgenden Schritte skizzieren den Restaurierungsalgorithmus, der hochfrequentes Rauschen entfernt und die relevanten Strukturen in dem Bild hervorhebt.
1. Führe eine 2-D FFT durch und ermittle die Phase Xp(fr, fc) und den Betrag XM(fr, fc).
2. Der Betrag wird zur Potenz p erhöht:
S(fr, fc) = XM(fr, fc) p (1)
3. Das Histogramm von S(fr, fc) wird zur Anpassung an eine Gauß'sche Funktion gedehnt: SG(fr, fc)
4. Ein Kontextbild C(fr, fc) wird durch einen Schwellwertvergleich von SG(fr, fc) gebildet, indem T% der Pixel beibehalten werden. Objekte mit weniger als n Pixeln werden aus dem Schwellwertbild entfernt. (Es sei bemerkt, daß mehr als T% der Energie beibehalten wird.)
5. Ein 'neues' Betragsbild XN(fr, fc) wird unter Verwendung des Kontextbildes und des Potenzspektrums (nicht des Betrags) gebildet:
6. Das neue Betragsbild XN(fr, fc) wird mit dem Phasenbild Xp(fr, fc) kombiniert und das wiederhergestellte Bild χR(r, c) wird als der Realteil der inversen Fourier-Transformation erhalten. Es wird bevorzugt, den Schwellwertvergleich im Bereich T = 1-6 und p = 1-2 anzuwenden.
Es sei bemerkt, daß der zur Potenz p erhöhte Betrag verwendet wird, um das neue Betragsbild zu bilden. Dies muß nicht naheliegend sein, aber die Erfahrung hat gezeigt, daß dies im allgemeinen sehr nützliche Ergebnisse liefert. Der Grund hierfür ist, daß beim Erhöhen des Betrags auf eine gegebene Potenz größer als eins die Teile des Spektrums, die die meiste Energie enthalten, gewöhnlich hervorgehoben werden, während Spektralkomponenten mit weniger Energie gewöhnlich weiter abgeschwächt werden. Das Ergebnis der Hervorhebungsprozedur unter Verwendung von T = 2 und p = 2 ist in Fig. 6 gezeigt.

Extraktion des PTM

Nachfolgend werden verschiedene Bildtransformationen und Texturmerkmale diskutiert, die aus transformierten Bildern extrahiert werden. Mehrere Merkmale basieren auf einem Bild, das nachfolgend als das projizierte trabekulare Mu ster (PTM) bezeichnet wird, das wie nachstehend skizziert erhalten wird:
Algorithmus 2 Die nachstehend skizzierten Schritte beschreiben die gewünschten Transformationen des wiederhergestellten Bildes, um einen relevanten Merkmalssatz auf dem resultierenden PTM zu messen.
1. Führe die morphologische Graustufen-Top-Hat- Operation auf dem wiederhergestellten Bild durch. Die Top-Hat-Operation ist das ursprüngliche Bild minus einer geöffneten Version des Bildes. Dieses Bild wird nachfolgend mit einem Schwellwert verglichen. Das Ergebnis dieser Operationen ist in Fig. 7 gezeigt.
2. Dieses Bild enthält deutlich eine Menge Rausch- Pixel. Kleine isolierte Objekte (in dem vorliegenden Beispiel alle Objekte mit weniger als 10 Pixeln) werden entfernt. Das resultierende PTM ist in Fig. 8 gezeigt.

Extraktion statistischer Merkmale

Natürlicherweise wird bevorzugt, Merkmale zu extrahieren, die einigermaßen die Eigenschaften (wie Struktur und Dichte) der trabekularen Struktur quantifizieren. Nicht überraschend sind die Maschen des PTM-Netzes viel größer und oftmals gebrochen für osteoporotische Individuen als für nicht-osteoporotische Individuen. Dieser Unterschied kann auf mehrere Weisen quantifiziert werden. Solche Eigenschaften sind beispielsweise in dem mittleren Abstand von jedem Hintergrundpixel zu dem nächsten Vordergrundpixel (Knochen- Pixel) in dem projizierten trabekularen Muster wiedergegeben.
Die Größe der dunklen Maschen in dem Muster ist ein relevanter Parameter, der geeigneterweise quantifiziert wird, indem die sog. Abstandstransformation für das projizierte trabekulare Muster berechnet wird. Das Ergebnis der Abstandstransformation auf dem PTM ist in Fig. 9 gezeigt und mit χD(r, c) bezeichnet.
Zusätzlich wird für jeden Hintergrundpixel in dem projizierten trabekularen Muster eine Linie von dem Pixel in eine gegebene Richtung gezeichnet, bis ein Vordergrundpixel ('Knochen-Pixel') gefunden wird. Der Abstand zwischen diesen Pixeln wird gemessen. Dies wird für eine festgelegte Anzahl an Richtungen durchgeführt, wie in Fig. 10 dargestellt. Danach kann der mittlere Abstand, der als die Sternfläche bezeichnet wird, als χA(r, c) angegeben, oder die maximale Sternlänge, als χX(r, c) angegeben, berechnet werden. Diese Bilder sind jeweils in Fig. 11 und Fig. 12 gezeigt.
Für alle Hintergrundpixel kann eine Anzahl statistischer Momente aus den resultierenden Bildern aus dem gemessenen Mittelwert, der Standardabweichung, dem Variationskoeffizienten und dem Schiefe-Koeffizienten berechnet werden:
wobei das tiefgestellte Zeichen verschiedene betrachtete Transformationen bezeichnet:
- D: Abstandstransformation
- M: Mittelwert der Sternlänge
- S: Standardabweichung der Sternlänge
- X: Maximale Sternlänge
- A: Sternfläche
Wir skizzieren kurz die Interpretation der betrachteten Parameter:
Die Merkmale und . sind entweder auf die allgemeine Maschengröße (für D, M, A und X) oder die Verschlechterung bzw. den Abbau der Maschen (5) bezogen. Beispielsweise sind größere Maschengrößen in einem größeren Wert von . wiedergegeben. Somit sollte für steigende Werte von . und . der Texturindex absinken.
Das Merkmal cv. ist etwas unterschiedlich. Es mißt die Breite der empirischen Verteilung in den transformierten Bildern. Falls all die Abstände groß sind, würde cv. klein sein. Das selbe würde der Fall sein, falls all die Abstände klein wären. Typischerweise finden wir jedoch, daß die empirische Verteilung breit ist, lediglich wenn das PTM-Netz gebrochen ist. Somit sollte für steigende Werte von c . der Texturindex typischerweise sinken.
Das Merkmal . bezieht sich ebenfalls auf die Maschengröße. Falls das Histogramm der Abstände auf der rechten Seite schwerer ist, ist .< 0. Falls die empirische Verteilung symmetrisch ist, ist . = 0. Für große Maschen in dem PTM würde die Verteilung der Abstände typischerweise für große Werte schief sein. Somit sollte für steigende Werte von . der Texturindex sinken.

In dem Ortsfrequenzbereich extrahierte Texturmerkmale

Man würde intuitiv denken, daß in dem Frequenzbereich gemessene Merkmale relevante Information über die trabekulare Struktur tragen sollten. Die trabekulare Struktur für stark osteoporotische Patienten ist deutlich niederfrequenter als die für 'normale' Patienten. In dem vorliegenden Abschnitt werden die betrachteten Fourier-Merkmale nicht an dem projizierten trabekularen Muster gemessen. Statt dessen werden die Merkmale an dem Spektrum des hintergrundkorrigierten Bildes gemessen.

Vorverarbeitung

Bevor das Spektrum geschätzt wird, werden Nicht- Stationaritäten des extrahierten Bildes durch Subtrahieren einer 31 · 31 median-gefilterten Version des Bildes von sich selbst entfernt, wie oben beschrieben.

Spektralschätzung

Da das Hauptinteresse in Richtung der relativ niederfrequenten Eigenschaften des Leistungsspektrums gerichtet ist, kann eine gewichtete Schätzung der Dispersion verwendet werden, wobei p*(u, v) das normalisierte Leistungsspektrum ist, das eine Gewichtungsfunktion berücksichtigt:
und
w(u,v) = 1 / u² + v²
Gegeben sei eine positiv definite symmetrische Matrix
Die Eigenwerte sind
und die entsprechenden (normalisierten) Eigenvektoren
Eine alternative Schätzung des Leistungsspektrums kann aus dem folgenden Prozeß erhalten werden:
Der wahrscheinlich bestbekannteste Schätzer des Leistungsspektrums ist der Periodogramm-Schätzer:
wobei (u, v) Ortsfrequenzen sind.
Es ist bekannt, daß diese Schätzung eine nicht-konsistente Schätzung des Leistungsspektrums ist. Wenn beispielsweise die Proben aus einem Gauß'schen Prozeß weißen Rauschens verwendet werden, hat die Periodogramm-Schätzung eine Varianz, die proportional zu der Varianz des Prozesses weißen Rauschens ist.
Ein Ansatz für die Verminderung der Varianz ist, die Daten zu segmentieren, z. B. in nicht-überlappende Blöcke, wobei das Spektrum in jedem Block unter Anwendung dieses Periodogramm-Schätzers geschätzt wird, und die endgültige Schätzung erhalten wird, indem die Periodogramme der Blöcke gemittelt werden. Auf diese Weise wird die Varianz um einen Faktor gesenkt, der gleich der Anzahl der Blöcke ist. Selbstverständlich geht dies auf Kosten der Frequenzauflösung.
Nachfolgend wird eine iterative Prozedur für die Schätzung der Dispersionsmatrix vorgeschlagen, deren Schätzung empfindlicher auf die Anisotropien ist, die in den niederfrequenten Teilen der Orientierung der Textur wiedergegeben sind.
Algorithmus 3 Die Dispersionsmatrix der genäherten Gauß'schen Funktion wird iterativ auf folgende Weise geschätzt:
1. Σ&sub0; = I
2. i = 1
3. wiederhole
3a. Schätze die Dispersion:
wobei
und
3b. i = i-1
bis Σi - Σi-1 max < δ
Aus diesen beiden Techniken kann die Hauptrichtung des Leistungsspektrums ermittelt werden unter Verwendung von:
Diese Hauptrichtung kann verwendet werden, die Ellipse mit der kleinstmöglichen Fläche zu finden, die einen festgelegten Bruchteil der Energie des ursprünglichen (nichtgewichteten) Leistungsspektrums enthält. Die mit dem Winkel α um den Ursprung rotierte Ellipse wird in Polarkoordinaten ausgedrückt:
R(θ)² = a²b² / b²cos²(θ-α)+a²sin²(θ-α)
wobei a und b jeweils die Haupthalbachse und die Nebenhalbachse sind und wobei die oben erwähnten λ's die Eigenvektoren sind.
In dem vorliegenden Beispiel können die Haupthalbachse und die Nebenhalbachse gefunden werden, indem das folgende Schema angewandt wird:
Algorithmus 4 Die nachstehend skizzierten Schritte werden zum Schätzen der Haupthalbachse und der Nebenhalbachse in der Ellipse mit der kleinstmöglichen Fläche verwendet, die einen gegebenen Bruchteil der spektralen Energie enthält. Der Winkel, mit dem diese Ellipse um den Ur sprung rotiert werden sollte, wird wie oben skizziert ermittelt.
1. Starte mit einer ausreichend kleinen Nebenhalbachse und Haupthalbachse, z. B. a = b = 10.
2. Erhöhe a und b um eins und berechne den Bruchteil der Energie in jeder Ellipse (mit dem Winkel α um den Ursprung rotiert).
3. Die Ellipse mit dem größten Bruchteil der Energie wird ausgewählt.
4. Wiederhole die Schritte 2 und 3, bis der gesuchte Bruchteil der Energie erreicht ist.
Die hier betrachteten Merkmale sind die geschätzten Werte:
â, , â- , /â sowie die Merkmale λ&sub1;, λ&sub2;, λ&sub2;/λ&sub1;, ρ = σuv/σuσv, σv²/σu² und
In Fig. 13 ist die geschätzte Ellipse, die 70% der Energie enthält, (und die Richtung) für einen 'normalen' Patienten gezeigt. In Fig. 14 ist die 70% Ellipse für einen osteoporotischen Patienten gezeigt. Der Frequenzgehalt für den osteoporotischen Patienten ist deutlich weit mehr niederfrequent. Dies bedeutet, daß die Energie nahe dem Zentrum des Bildes konzentriert ist und folglich die Flächen der Ellipsen sehr unterschiedlich sind, obwohl sie den selben Bruchteil der Energie enthalten.

Ortsinteraktionsmodelle

Oben ist beschrieben worden, daß ein starker Satz an Merkmalen aus dem Ortsfrequenzbereich extrahiert werden kann. Die Leistungsfähigkeit der in dem Frequenzbereich extrahierten Merkmale hängt jedoch stark von dem verwendeten Spektralschätzer ab. Das Periodogramm ist ein nichtkonsistenter Spektralschätzer. Dieses Problem kann behoben werden, indem die Daten in nicht-überlappende Segmente aufgeteilt werden, das Periodogramm in jedem Segment geschätzt wird und schließlich über diese Periodogramme gemittelt wird.
Die Varianzreduktion ist direkt proportional zu der Anzahl an Segmenten. Dies geht jedoch auf Kosten der Frequenzauflösung.
Nachfolgend werden wir zwei Typen kausaler simultaner autoregressiver (SAR = "Simultaneous AutoRegressive", d. h. simultane auto-regressive) Modelle für den Erhalt von parametrischen Spektralschätzungen betrachtet; das Modell mit nicht-symmetrischen Halbebenen (NSHP = "Non-Symmetrical Halfplane", d. h. mit nicht-symmetrischen Halbebenen) und das Viertelebenen-Modell (QP = "Quarter-Plane", d. h. Viertelebenen).
Es wird ein über weite Bereiche stationärer stochastischer Prozeß mit dem Mittelwert null angenommen, der auf einem rechtwinkligen Gitter an Pixelplätzen Ω = {s = (sr, sc): 0 ≤ sr ≤ M, 0 ≤ sc ≤ N} definiert ist, und einem SARMA-Model gehorcht, das definiert ist als:
wobei r = (rr, rc) NAR oder NMA ist. NAR und NMA sind jeweils die Trägerregionen für die Parameter-Arrays {φ} und {θ}.
Wir beziehen auf diese Sätze ebenfalls einen Nachbarsatz, der einigermaßen die Ordnung des Modells definiert.
Es wird angenommen, daß der Innovationsprozeß s ein IID N(0, Σ²) oder eine Gauß'sche Folge weißen Rauschens ist. Falls NMA = ist, haben wir ein reines SAR-Modell und falls NAR = ist, haben wir ein reines SMA-Modell.
Nachfolgend werden wir aus folgenden Gründen auf reine Auto-Regressionen fokussieren: Die Verwendung gemischter SARMA-Modelle würde wahrscheinlich geizigere Modelle geben, die zum Einfangen der Natur der betrachteten Textur erforderlichen Modelle sind jedoch von einer solch niedrigen Ordnung, daß Geiz nicht wirklich ein wichtiger Gesichtspunkt ist. Es stellt sich ebenfalls heraus, daß das z. Zt. wichtige Merkmal die Form des Spektrums in Ausdrücken von beispielsweise Peak-Stellen ist. Der LS-Schätzer für das reine SAR-Modell ist ebenfalls ein Satz linearer Gleichungen, während eine nicht-lineare Optimierung erforderlich ist, um Schätzungen gemischter SARMA-Modelle zu erhalten.
Die (nicht-normalisierte) spektrale Dichtefunktion ist gegeben als das quadrierte Modul der Frequenzantwortfunktion multipliziert mit der Rauschvarianz:
wobei die Ortsfrequenzen (fr, fc) = f Ωf sind, und
Ωf = (f = fr, fc) = (-1/2 -r/M, -1/2 +c/M) / (r, c) Ω (27)
Die Vektornotation ist definiert als:
φ = col[φr, r NAR]
θ = col[θr, r NMA]
(col = Spalte)
C(N.)f = col[cos(ω(f, r)),r N.]
S(N.)f = col[sin(ω(f, r)), r N.]
ω(f, r) = 2π(rr · ff + rc · fc)
Bis jetzt ist die Frage der Kausalität noch nicht angesprochen worden. In zwei Dimensionen ist die Vorstellung der Kausalität nicht natürlich. Ein Aufgeben der Kausalität führt jedoch zu einer Reihe Schwierigkeiten, von denen eine ist, daß der Schätzer nach den kleinsten Quadraten der Parameter nicht länger konsistent ist. Folglich ist es populär geworden, dem Bild eine künstliche Richtungsabhängigkeit aufzuerlegen, um somit ein kausales Modell zu erhalten (tatsächlich ist rekursiv berechenbar der korrektere Ausdruck).
Die beiden Typen in der Literatur typischerweise betrachteten kausalen Modelle sind die Viertelebenen-Trägermodelle (QP) und die Trägermodelle mit nicht-symmetrischen Halbebenen (NSHP). In Fig. 22 sind parametrische Schätzungen unter Anwendung des NSHP-Modells dargestellt. Es ist ersichtlich, daß das Spektrum in einer bestimmten Richtung leicht in die Länge gezogen ist, im Vergleich zu der entsprechenden aus Fig. 21 ersichtlichen Periodogramm-Schätzung. Dieses Phänomen hängt von dem richtungsabhängigen Bias ab, was ein Effekt ist, der für den QP-Modus sogar noch stärker hervortritt. Um dieses Problem zu lösen, ist ein gewichteter QP- Schätzet (WQP = "Weighted Quarter Plane", d. h. gewichteter Viertelebenen) vorgeschlagen worden:
wobei die tiefgestellten Zeichen ++ und -+ einfach den betrachteten Quadranten bezeichnen, d. h. den ersten und den zweiten Quadranten. Es herrscht allgemein die Meinung, daß dieser parametrische Schätzer einer unter den besten parametrischen spektralen Schätzern zum Schätzen von Sinusoiden ist, die in weißem Rauschen eingebettet sind (KAY, S. M., 1988, "Modern Spectral Estimation", Prentice-Hall, Signal Proc. Series). Er hat jedoch den Nachteil, daß er keine Modellinterpretation hat.
In Fig. 23 ist eine WPQ-Schätzung dargestellt, bei der kein richtungsabhängiger Bias ersichtlich ist.
Unter Verwendung dieser Schätzungen können Merkmale des oben beschriebenen Typs, d. h. λ&sub1;, λ&sub2;, λ&sub2;/λ&sub1;, ρ = σrc/σrσc, σc²/σr² und.
verwendet werden.
Da α die Hauptorientierung der genäherten Gauß'schen Funktion beschreibt, kann diese Richtung verwendet werden, die oben erwähnten Energieellipsen zu finden. Somit können ebenfalls die Haupthalbachse und die Nebenhalbachse dieser Ellipsen â und als auch die Merkmale â· , /â verwendet werden.

Extraktion von Merkmalen bezüglich der Periodizität

Ein alternatives Verfahren zum Ableiten von Merkmalen aus dem vorverarbeiteten Bild der Fig. 6 ist wiederum, dieses Bild einer zweidimensionalen Fourier-Transformation zu unterziehen und das Fourier-Leistungsspektrum (das Periodogramm-Spektrum) zu extrahieren.
Es wird z. Zt. bevorzugt, daß die Auflösung des Bildes unter Anwendung einer Gauß'schen Pyramide reduziert wird, wodurch so viel wie möglich an Information in der hohen Auflösung vor der Fourier-Transformation bewahrt wird.

Merkmalsextraktion

Es ist nicht überraschend, daß das Leistungsspektrum Peaks umfaßt, welche der Periodizität der vertikalen Trabekulae entsprechen. Es ist jedoch erstaunlich, daß ebenfalls Peaks - jedoch klein -, die einer Periodizität der horizontalen Trabekulae entsprechen, in dem Leistungsspektrum ausgemacht werden können.
Information bezüglich der Dichte der Trabekulae, deren Dicke und des Abstands zwischen diesen kann aus der Höhe dieser sich auf die Periodizitäten beziehenden Peaks, der Gesamtfläche unter diesen Peaks und/oder deren "Breite" extrahiert werden.
Zur Zeit wird bevorzugt, die Peak-Positionen nicht aus dem Fourier-Leistungsspektrum zu extrahieren, da das Periodogramm oder das Leistungsspektrum sehr rauschbelastet ist, aufgrund der Tatsache, daß der Periodogramm-Schätzer des Leistungsspektrums nicht-konsistent ist. Statt dessen wird ein kombinierter Algorithmus bevorzugt, bei dem das parametrische Spektrum zum Identifizieren dieser Peak-Positionen verwendet wird und das geglättete Periodogramm zum Schätzen der zugehörigen Energie etc. verwendet wird.
In dem Leistungsspektrum liegen typischerweise insgesamt vier Peaks vor - oder genauer gesagt, lediglich zwei Peaks, da das parametrische Spektrum symmetrisch ist. Eine Linie wird durch die zwei größten Peaks (die sich auf die vertikale Trabekulae beziehen) durch den Wert bei der Frequenz null gezogen. Von der Mitte dieser Linie wird eine Linie senkrecht dazu ebenfalls durch den Wert bei der Frequenz null gezeichnet. Es liegt in der Natur des Modells, daß die Peaks der horizontalen Trabekulae rechtwinklig zu denjenigen der vertikalen Trabekulae angeordnet sind. Die Maxima auf den neuen Linien werden den Peaks der horizontalen Trabekulae entsprechen.
Dieses Verfahren wird bevorzugt, damit es unabhängig und unempfindlich gegenüber irgendeiner Rotation des Bildes ist.
Ein Glätten des Spektrums kann unter Anwendung eines k·k- Kerns, wie einem 5 · 5-Kern, unter Verwendung beispielsweise eines Mittelwertfilters, eines Median-Filters, eines Gauß'schen Filters etc. durchgeführt werden.
Aus diesem Spektrum können die Peakhöhen, das Volumen unter den Peaks usw. bestimmt und bei der Vorhersage der Knochenqualität verwendet werden.
Zusätzlich wird in Erwägung gezogen, daß ebenfalls die Krümmung der Peaks in dem parametrischen Leistungsspektrum als Vorhersage der Knochenqualität verwendet werden kann, da sich diese Krümmung ebenfalls auf Abweichungen der Periodizität bezieht.

BEISPIEL 2: Studie der Merkmalsleistungsfähigkeit

Es wurde ein Datensatz verwendet, um die verschiedenen Bildmerkmale zu untersuchen und ihre Leistungsfähigkeit miteinander zu vergleichen. Es wurden Röntgenbilder von 97 Individuen erhalten und jedes Bild wurde von erfahrenen Radiologen bewertet. Die Analyse dieses Datensatzes ist nachfolgend beschrieben.
Jedem Patienten wurde das Ergebnis +1 zugeordnet, falls der Radiologe den Patienten als 'normal' betrachtete, und das Ergebnis -1 zugeordnet im Falle eines beschleunigten Knochenverlustes (osteoporotisch). Die vier Radiologen stimmten bei der Diagnose in 52 Fällen überein, und lediglich diese wurden in dieser vorläufigen Studie verwendet. In dieser Studie der Merkmalsleistungsfähigkeit wird ein allgemeines lineares Modell mit dem Ergebnis als Antwortvariable und dem Alter, dem Geschlecht und den Texturmerkmalen als erklärende Variablen verwendet. Ferner sind rigorose statistische Tests der Signifikanz der individuellen Merkmale unter Anwendung einer Rückwärtseliminations-Prozedur durchgeführt worden.
Die verwendeten Merkmale sind oben beschrieben worden und basieren alle auf Statistiken der ersten, der zweiten und der dritten Ordnung, die an den Transformationen des projizierten trabekularen Musters (PTM) oder angepaßten Ellipsen auf dem Fourier-Leistungsspektrum der wiederhergestellten Bilder gemessen wurden.
Die angepaßten Modelle sind vom Typ
si = β&sub0; + β&sub1;Alteri + β&sub2;Geschlechti + β&sub3;f&sub1;,i + ··· + i
wobei si das Ergebnis des i-ten Patienten ist, Alteri das Alter des i-ten Patienten ist, Geschlechti das Geschlecht des i-ten Patienten ist (+1: männlich, -1: weiblich, um beide Geschlechter gleichmäßig am Mittelwert beitragen zu lassen) und fj,i das j-te Bildmerkmal des i-ten Patienten ist. Es wird angenommen, daß ein Gauß'sches weißes Rauschen ist (??). Unter dieser Annahme sind die LS- Schätzungen gleich den ML-Schätzungen. Es werden Teststatistiken für jeden Parameter in dem Modell berechnet. Unter Anwendung eines Rückwärtseliminationsansatzes wird der am wenigsten signifikante Parameter des Modells eliminiert und das Modell wird nachfolgend mit den verbleibenden Parametern wieder neu geschätzt. Dies wird so lange fortgeführt, bis kein Parameter mehr aus dem Modell eliminiert werden kann.

Merkmale basierend auf dem projizierten trabekularen Muster

Unter Anwendung der Rückwärtseliminationsprozedur wurden von 12 Bildmerkmalen 9 eliminiert. Die Reihenfolge, in der diese eliminiert wurden, ist in der nachstehenden Tabelle zusammengefaßt:
Die fettgedruckten, mit einem Stern markierten Merkmale sind signifikante Parameter, die nicht eliminiert werden konnten.
Das resultierende Modell ist gegeben durch:
si = β&sub0; + β&sub1;Alteri + β&sub2;Geschlechti + β&sub3;uD,i + β&sub4;c D,i + β&sub5;σs,i + i
Eine Zusammenfassung des resultierenden Modells ist nachstehend gegeben:
call: Im(formula = Ergebnis[idx4] ~ Alter + Geschlecht [idx4] + f[idx4, 1] + f[idx4, 3] + f[idx4, 6])
Reste bzw. Residuen:
Reststandardfehler: 1,746 auf 45 Freiheitsgrade
Mehrfach-R-quadriert: 0,7516
Korrelation der Koeffizienten:
(Achsenabsch.)Alter Geschlecht[idx4] f[idx4,1] f[idx4,3]
Alter 0,0063
Von dem erklärenden Effekt ist das meiste im Alter und im Geschlecht enthalten. Der Rest der Variabilität (ungefähr 16% Punkte) wird durch die Texturmerkmale erklärt. Es kann erwünscht sein, einige der Merkmale zu transformieren oder sogar neue abzuleiten.
Aus dem obenstehenden kann geschlossen werden, daß die betrachteten Statistiken, die wenigstens an der Abstandstransformation und der Sternfläche des projizierten trabekularen Musters gemessen werden, ziemlich starke Merkmale sind. Diese Untersuchung schlägt jedoch nicht die Maximalabstands-Transformation als ein optimales Merkmal vor.

Spektralmerkmale

Nachfolgend werden vier aus dem Leistungsspektrum extrahierte Merkmale betrachtet. Das einzige aus diesem Merkmalssatz eliminierte Merkmal ist die Länge der Nebenhalbachse (â). Dies steht in guter Übereinstimmung mit der Erwartung, da die Auswaschung der trabekularen Struktur bei osteoporotischen Patienten hauptsächlich in der vertikalen Richtung stattfindet. Die Leistungsfähigkeit des resultierenden Modells ist mit den Ergebnissen vergleichbar, die unter Verwendung von auf der Abstandstransformation basierenden Merkmalen erhalten wurden.
Wiederum wurden die Daten in dem statistischen Software- Paket Splus analysiert und es wurde die folgende Ausgabe für das resultierende Modell erhalten:
Reste bzw. Residuen:
Reststandardfehler: 1,843 auf 45 Freiheitsgraden Mehrfach-R-quadriert: 0,7232
Korrelation von Koeffizienten:
Aus den obigen Ergebnissen ist ersichtlich, daß ein Mehrfach-Korrelationskoeffizient besser als 0,7 erhalten werden kann, indem entweder die aus dem Leistungsspektrum extrahierten Merkmale oder die aus dem PTM extrahierten Merkmale verwendet werden; lediglich ein Teil der von dem Verfahren der Erfindung angebotenen Merkmale. Es wird somit erwogen, daß selbst höhere Mehrfach-Korrelationskoeffizienten erwartet werden können, wenn diese Merkmalsgruppen in den Schätzprozeduren kombiniert werden. In diesem Zusammenhang sei bemerkt, daß Mehrfach-Korrelationskoeffizienten von dieser Größe ziemlich zufriedenstellend sind.

BEISPIEL 3: Der MORSE-Schätzer für gemischte nicht-kausale SARMA-Modelle

Es herrscht allgemein die Ansicht, daß das Texturkonzept in vielen Anwendungsgebieten äußerst wichtig ist, wie z. B. medizinischer Bildwiedergabe, industrieller Qualitätsuntersuchung und Fernerfassung. Visuelle Eigenschaften einer gegebenen Textur sind jedoch nicht sehr greifbar und viele Versuche sind unternommen worden, gute und bedeutungsvolle Merkmale zu extrahieren, um solche Eigenschaften zu beschreiben.
Ein solcher Ansatz ist, die Ortskorrelationsstruktur unter Verwendung von Texturmodellen zu modellieren. Die beiden Hauptmodelltypen sind die Gauß-Markov-Zufallsfeld-Modelle (GMRF = "Gaussian Markov Random Field", d. h. Gauß-Markov- Zufallsfeld) und die simultanen auto-regressiven Modelle mit gleitendem Mittelwert (SARMA) ("Simultaneous Auto- Regressive Moving Average", d. h. simultane auto-regressive mit gleitendem Mittelwert). Der letztere Modelltyp ist ein Versuch, die sehr mächtigen ARMA-Modelle (aus der Zeitreihenanalyse und der eindimensionalen Signalverarbeitung bekannt) auf zwei Dimensionen zu erweitern.
Während eindimensionale Signale von Natur aus kausal sind, ist dies keine natürliche Beschränkung in zwei Dimensionen. Das Aufgeben von Kausalität führt jedoch zu ernsthaften Problemen bei wichtigen Modellierungsschritten, wie der Identifikation und der Schätzung. Daher ist es sehr geläufig, den Daten eine künstliche Richtungsabhängigkeit aufzuerlegen und dadurch konsistente Schätzungen zu erhalten, indem die Summe quadrierter Schätzungsfehler - d. h. Schätzungen nach den kleinsten Quadraten - minimiert wird.
Die Analyse nicht-kausaler SARMA-Modelle ist stark eingeschränkt worden, hauptsächlich aufgrund der Tatsache, daß die Nicht-Konsistenz des LS-Schätzers einen zwingt, eine sehr komplizierte Likelihood-Funktion zu berücksichtigen. Als Konsequenz hat sich die Analyse allein auf reine SAR- Modelle konzentriert. Für reine SAR-Modelle sind der toroidale ML-Schätzer und ein genäherter ML-Schätzer formuliert worden. Es kann gezeigt werden, daß die Leistungsfähigkeit der genäherten Schätzer nicht geeignet sein kann, wenn die betrachteten Texturen nahe der Nicht-Stationarität liegen.
Eine Berücksichtigung der exakten ML-Schätzungen für toroidale SAR-Modelle gibt im allgemeinen sehr gute Ergebnisse. Es ist jedoch gezeigt worden, daß die Likelihood-Funktion in der Nähe der Nicht-Stationarität extrem wellig ist. Dies führt in vielen Fällen dazu, daß die Optimierungsroutine bei einer sub-optimalen Lösung "stecken" bleibt. Dieses Problem kann tatsächlich ernsthaft werden, da visuelle Texturen durch eine starke positive Korrelationsstruktur charakterisiert sind, die Schätzungen nahe der Nicht- Stationarität ergeben.
Ein weiteres Problem, das sich auf den Maximum-Likelihood- Schätzer und die traditionellen Maximum-Entropie-Schätzer bezieht, ist, daß diese den hochfrequenten Gehalt des Betrachteten gewöhnlich überbetonen.
Daher wird ein neuer Entropie-basierender Schätzer eingeführt. Anstatt die Entropie des zugehörigen Modellspektrums zu maximieren (wie in der MEM-Spektralschätzung, eingeführt in "A new algorithm for two-dimensional maximum entropy power spectrum estimation", Lim, J. S. & Malik, N. A. (1981) IEEE Transactions on Acoustics, Speech and Signal Processing, 29, 401-413), wird die Entropie des Restspektrums maximiert. Dieser Schätzer wird MORSE-Schätzer genannt ("Maximum Of Residual Spectral Entropy", d. h. Maximum der Restspektralentropie). Dieser Schätzer verhält sich sogar nahe der Nicht-Stationarität gut, was es viel leichter macht, die nicht-lineare Optimierung des einbezogenen Zielfunktionals (bezüglich der Parameter) durchzuführen. Es ist möglich, das Zielfunktional für gemischte SARMA-Modelle (d. h. einschließlich der Teile mit gleitendem Mittelwert) zu schreiben und sehr gute Schätzungen zu erhalten. Nach unserem besten Wissen gibt es keine Ergebnisse in der Literatur, welche die Folgerung auf nicht-kausale SARMA-Modelle berücksichtigen.
Die Maximierung der Restspektralentropie scheint intuitiv ein vernünftiges Kriterium zu sein, da die zugrundeliegende Annahme der simultanen Modelle ist, daß der Innovationsprozeß weiß ist. Es ist offensichtlich, daß das Optimum eines gegebenen Modells eng auf das ausgewählte Kriterium des Optimums bezogen ist und die Wahl der Entropie keineswegs die einzige vorstellbare ist. Es sei bemerkt, daß keine Annahmen bezüglich der Folge des antreibenden Rauschens (Innovationsprozeß ????) gemacht werden. Wir verlangen lediglich, daß es eine symmetrische Verteilung hat.

Das nicht-kausale SARMA-Modell

Betrachte einen weitläufig stationären, bilateralen SARMA- Prozeß mit Mittelwert null, der auf einem quadratischen MxM-Pixelgitter s = (r, c) Ω abgetastet wird
wobei r = (r&sub1;, r&sub2;) N, ist und N. die Ordnung des Modells definiert. Die einzige Annahme, die über den Innovationsprozeß {y.} gemacht wird, ist, daß dieser unabhängige identisch verteilte Zufallsvariablen mit einer symmetrischen Dichte hat.
Im Gegensatz zu den GMRF-Modellen muß das bilaterale SARMA- Modell nicht symmetrisch sein, sondern gegenüberliegende Nachbarn müssen identischen Parametern zugehören (d. h. φr = φ-r). Dies ist erforderlich, um die Identifizierbarkeit der Parameter sicherzustellen. Aus Gründen der Einfachheit wird sich nachfolgend darauf beschränkt, lediglich symmetrische Modelle zu betrachten.
Zum Vereinfachen der Notation, falls gewünscht wird, die Ordnung des Modells genauer anzugeben, wird die in den folgenden Definitionen skizzierte Notation eingeführt:
Definition 1 Ein Rechteck mit den Seitenlängen 2a+1 und 2b+1, das um einen Winkel α um den Ursprung rotiert ist, wird mit R(a, b, α) bezeichnet. Nunmehr wird folgendes definiert
· R(a, b) = R(a, b, 0)
· R(a) = R(a, a, 0)
Definition 2 Eine Ellipse mit der Haupthalbachse a und der Nebenhalbachse b, die um einen Winkel α um den Ursprung rotiert ist, wird mit E(a, b, α) bezeichnet. Es wird folgendes definiert
· E(a, b) = E(a, b, 0)
· E(a) = E (a, a, 0)
Somit kann die Ordnung eines gegebenen Modells in Ausdrücken dieser geometrischen Formen definiert werden, wie beispielsweise SARMA(E(1),E(1)), was bedeutet, daß der autoregressive Teil sowie der Teil mit gleitendem Mittelwert die vier nächsten Nachbarn sind.
Falls das bilaterale Modell berücksichtigt wird, ist die Stationaritäts-Bedingung stark abgeschwächt; alles, was ge fordert wird, ist, daß keine Pole der Übertragungsfunktion auf den Einheitskreis fallen:
für alle Zr = 1 und Zc = 1.
Die Übertragungsfunktion ist als die zweidimensionale z- Transformation der Impulsantwortfunktion gegeben.
und die Frequenzantwortfunktion wird als die zweidimensionale Fourier-Transformation der Impulsantwortfunktion erhalten
H(f(s)) = (1+ΘTMAC(NMR)f(s)) / (1+ΘTARC(NAR)f(s))
wobei
f(s) = (r/M - 1/2, c/M -1/2)T
die Ortsfrequenz ist und
ΘMA = col[θr, r NMA]
ΘAR = col[φr, r NAR]
C(N.)f(s) = col[cos(2πf(s)Tr), r N.]
Bei der Zeitreihenanalyse und der Signalverarbeitung hat die erfolgreiche Verwendung von ARMA-Modellen für eine spektrale Schätzung ein Interesse für die Erweiterung die ser Verfahren auf zwei Dimensionen gefördert. Wie für ARMA- Modelle ist das zweidimensionale parametrische Spektrum proportional dem Betragsquadrat der Frequenzantwort:
S(f(s)) = σ² (f(s)) ²
wobei σ² die Varianz der Innovations-Folge ist.
Bei Anwendung dieses Ansatzes wird das Problem der Spektrumsschätzung, ein geeignetes Modell an die beobachteten Daten anzupassen. Meistens hat die parametrische spektrale Schätzung in zwei Dimensionen versucht, unilaterale Modelle zu verwenden. Dies führt typischerweise zu einem richtungsabhängigen Bias. Unter Verwendung eines bilateralen Modells sollte dieses spezielle Problem unter der Voraussetzung gelöst werden, daß ein praktisch anwendbares Schätzprinzip bereitgestellt wird. Unsere vorläufigen Ergebnisse zeigen, daß das MORSE-Schätzprinzip ein sehr nützliches Werkzeug für den Erhalt parametrischer spektraler Schätzungen ohne richtungsabhängigen Bias ist.
Verschiedene Ansätze sind verfolgt worden, um Schätzungen für das nicht-kausale SAR-Modell zu erhalten. Es gibt nach unserem besten Wissen keine Ergebnisse bezüglich der Schätzung gemischter nicht-kausaler SARMA-Modelle. Selbst für das reine SAR-Modell ist die Schätzung der Parameter kein einfaches Problem. Nachstehend werden verschiedene Schätzer für das nicht-kausale SAR-Modell dargestellt und ihre Nachteile diskutiert.
Die nicht-kausalen (oder bilateralen) simultanen Auto-Regressions-Modelle (SAR = "Simultaneous Auto-Regressions", d. h. simultane Auto-Regressions), die manchmal als NCSAR bezeichnet werden, wurden ursprünglich von P. Whittle in "On stationary processes in the plane", Biometrika, 41, 343-449 eingeführt. In dieser Druckschrift wird das Problem der Nicht-Konsistenz des LS-Schätzers behandelt, und es wird ein Ansatz für große Proben, basierend auf spektralen Methoden, vorgestellt. Erstens ist dieses Verfahren lediglich für sehr einfache Modelle verwendbar und zweitens umfassen viele praktische Situationen kleine Proben, für welche die Nichtberücksichtigung von Kanteneffekten nicht vernachlässigt werden kann. Für die Darstellung, wie die Nicht-Konsistenz des LS-Schätzers die Schätzungen beeinflußt, wird ein einfaches Experiment durchgeführt. Für ein isotropes Modell des nächsten Nachbarn werden 100 Realisierungen für verschiedene Parameter erzeugt. Die Parameter werden in jedem Fall geschätzt (als ob das Modell anisotrop wäre). In Fig. 15 ist der Mittelwert des geschätzten Parameters gegen den 'wahren' Parameter aufgetragen (d. h. der Parameter, mit dem die Textur simuliert wurde). Es ist ersichtlich, daß der LS-Schätzer keine zufriedenstellenden Schätzungen liefert. Tatsächlich liegt in vielen Fällen die Größenordnung um einen Faktor zwei daneben.
Kashyap, R. & Chellappa, R. schlugen eine Näherung für den exakten (toroidalen) ML-Schätzer in "Estimation and choice of neighbours in spatial interaction models of image", IEEE Transactions on Information Theory, 29(1), 60-72 vor. Die Näherung basiert auf einer Taylor-Reihenentwicklung der Determinante der Jacobi-Matrix (die blockumlaufend ("block circulant") mit umlaufenden Blöcken unter der Torusannahme wird). Es sind lediglich Teile erster und zweiter Ordnung (in den Parametern) einbezogen. Die Verbesserung der Ergebnisse ist definitiv umfassend, wie in Fig. 16 gezeigt ist, wo das obige Experiment wiederholt ist. Nichtsdestotrotz ist ersichtlich, da der Bias sehr deutlich nahe der Nicht- Stationarität ist, daß dieser Schätzer in vielen praktischen Situationen nicht geeignet sein kann, bei denen die Texturen gewöhnlich nahe der Nicht-Stationarität liegen, und es wird stark bevorzugt, einen geeigneten Schätzer ohne diese Nachteile zu ermitteln.
Grunkin untersuchte den exakten (toroidalen) ML-Schätzer in "On The Analysis of Image Data Using Simultaneous Interaction Models" in IMSOR, Doktorarbeit Nr. 67, Institute of Mathematical Statistics and Operations Research, Technical University of Denmark, Lyngby, Seiten 223ff.. Die einbezogene Likelihood-Funktion ist kompliziert und rechnerisch sehr aufwendig. Die Likelihood-Funktion ist ebenfalls in der Nähe der Nicht-Stationarität äußerst wellig, wie in Fig. 17 gezeigt. Ob die nicht-lineare Optimierungsroutine bei einer sub-optimalen Lösung "stecken" bleibt oder nicht, hängt sehr stark von der Anfangsschätzung und den Parametern des Optimierungsalgorithmus (wie z. B. die Schrittlänge) ab. Daher kann es mehrere erneute Schätzungen erfordern, um überzeugt zu sein, daß die korrekte Lösung gefunden ist. Dies ist selbstverständlich bloß ein technisches Problem und in dem einfachen, oben beschriebenen Experiment können diese Probleme behoben werden, und die exakten (toroidalen) ML-Schätzungen liefern sehr gute Ergebnisse, wie in Fig. 18 gezeigt.

Der MORSE-Schätzer

Wie oben erwähnt, führt die Nicht-Konsistenz des LS- Schätzers zu nicht-sinnvollen Ergebnissen. Die praktische Anwendung des exakten toroidalen ML-Schätzers wird in vielen Fällen eingeschränkt sein, da die Likelihood-Funktion in der Nähe der Nicht-Stationarität gewöhnlich extrem wellig ist. Daher kann es mehrere erneute Schätzungen erfordern, um sicherzustellen, daß das 'wahre' Optimum gefunden ist. Wie ebenfalls zuvor festgestellt, kann der ML-Schätzer den hochfrequenten Gehalt der Textur überbetonen. Außerdem ist die Rechenbelastung, die für den Erhalt von Schätzungen für gemischte SARMA-Modelle unter Verwendung der ML- Schätzung aufgebracht wird, für alle praktischen Zwecke untragbar.
Ein hier aufgenommener Ansatz liegt darin, den Spektralgehalt der Reste bzw. Residuen zu betrachten. Schließlich wird es wahrscheinlich extrem schwierig sein, ML- Schätzungen gemischter nicht-kausaler SARMA-Modelle zu erhalten.
In diesem Beispiel ist ein neuer Schätzer beschrieben, der sehr gute Schätzergebnisse liefert, eine Antwortfläche hat, die sich (nach unserer Erfahrung) sehr gut verhält und in gemischten SARMA-Modellen Schlußfolgerungen zuläßt.
Da die zugrundeliegende Annahme ist, daß die Reste bzw. Residuen weiß sein sollten, d. h. unkorreliert, sollten alle Frequenzen in dem Spektrum des Restprozesses mit dem gleichem Gewicht dargestellt werden. Falls das Restspektrum normalisiert ist, nimmt die Entropie des Restspektrums ihr Maximum an, wenn das Spektrum so flach wie möglich ist. Diese Beobachtung ist die Grundlage des Schätzers gemäß dem Maximum der Restspektralentropie. (MORSE "Maximum of Residual Spectral Entropy", d. h. Maximum der Restspektralentropie).
Betrachte wiederum die Differenzgleichung, die das System definiert:
Die Ausgabe kann als die Konvolution des Innovationsprozesses mit der Impulsantwort des Systems angesehen werden. Somit wird im Frequenzbereich die folgende Darstellung gefunden:
Y(f(s)) = (f(s)) · E(f(s)) (50)
wobei Y(f(s)) die Fourier-Transformierte der Ausgabe ist, E(f(s)) die Fourier-Transformierte des (nicht-observablen) Innovationsprozesses ist und (f(s)) die Fourier-Transformierte der Impulsantwortfunktion ist, die ebenfalls als die Frequenzantwortfunktion bekannt ist.
Unter Anwendung von Gleichung (50) kann das Leistungsspektrum des Innovationsprozesses in Ausdrücken des Leistungsspektrums der Beobachtungen und des parametrischen Spek trums ausgedrückt werden (unter der Annahme, daß die Parameter bekannt sind). Damit wird gefunden, daß:
Es kann eine Normalisierungsfunktion, die von den Beobachtungen und den Parametern abhängt, definiert werden:
Die Parameter werden nachfolgend ermittelt, indem die Entropie des Restspektrums maximiert wird.
Definition 3 Die Schätzungen gemäß dem Maximum der Restspektralentropie (MORSE-Schätzungen) werden ermittelt, indem die Entropiefunktion L(y, θAR, θMA) bezüglich der Parameter θAR und θMA, maximiert wird.
( AR, MA) = arg(ΘAR, ΘMA)max L(y, ΘAR, ΘMA) (53)
wobei
Eine Schätzung der Restvarianz wird erhalten als
² = N(y, ΘAR, ΘMA)/ M² (55)
Diese Funktion ist in den Parametern eindeutig stark nichtlinear. Die Optimierung kann durchgeführt werden, indem z. B. eine Quasi-Newton-Methode angewandt wird. In dieser Arbeit wird die IMSL-Implementation der BFGS-Methode angewandt.

Experimentelle Ergebnisse

Wie oben dargestellt, liefern die exakten (Torus-)ML- Schätzungen sehr gute Ergebnisse. Das Problem ist hier mehr von technischer Natur. Erstens ist die einbezogene Likelihood-Funktion (oder das Zielfunktional) extrem kompliziert und sehr aufwendig zu berechnen. Zweitens ist die Likelihood-Funktion in der Nähe der Nicht-Stationarität stark wellig. Folglich kann eine Optimierungs-Routine leicht bei einer sub-optimalen Lösung 'steckenbleiben'.
Oben sind wichtige Probleme, wie die Zentralität, Konsistenz und Wirksamkeit des vorgeschlagenen Schätzers nicht behandelt worden. Offensichtlich sind diese Fragen sehr schwierig analytisch zu beantworten, aufgrund der unbearbeitbaren Form der einbezogenen Ausdrücke. Unten werden diese Fragen empirisch durch Simulation und Schätzexperimente behandelt.

Schätzung reiner SAR-Modelle

In Fig. 19 ist das Zielfunktional des MORSE-Schätzers in der selben Region wie für den ML-Schätzer gezeigt. Es ist deutlich, daß das Zielfunktional des MORSE-Schätzers sehr viel glatter ist - d. h., daß dort keine Wellen sind. Unserer Erfahrung nach zeigt das Zielfunktional des MORSE- Schätzers niemals Wellen im Parameterraum - nicht einmal außerhalb der zulässigen Region. Dies macht die nichtlineare Optimierung zu einer einfachen Aufgabe, die von den meisten Standard-Softwarepaketen automatisch vorgenommen werden kann.
Das für die verschiedenen Schätzer oben in Verbindung mit Fig. 17 durchgeführte Experiment ist ebenfalls für den MORSE-Schätzer durchgeführt worden. Das Ergebnis ist in Fig. 20 gezeigt und es ist ersichtlich, daß die Ergebnisse sehr akzeptabel sind. Es sei bemerkt, daß die Varianz bei Annähern an die Nicht-Stationarität abnimmt. Dies hängt von der Tatsache ab, daß eine Annäherung an die Pole der Übertragungsfunktion stattfindet. Solche Pole (oder Nah-Pole) erzeugen scharfe Peaks im Zielfunktional und dadurch eine sehr viel kleinere Varianz.

Schätzung gemischter SARMA-Modelle

Ein Hauptvorteil dieses Ansatzes liegt darin, daß ein Schätzer für gemischte SARMA-Modelle 'umsonst' erhalten wird. Um zu zeigen, daß für gemischte SARMA-Modelle ebenfalls gute Ergebnisse erhalten werden, wird das folgende Experiment durchgeführt. Es wird ein Satz Parameter für ein SARMA(E(1.0),E(1,0)-Modell ausgewählt. Diese Parameter sind ferner nicht zu weit von der Nicht-Stationarität entfernt gewählt worden, so daß starke positive räumliche Auto- Kovarianzen erhalten werden. Nunmehr sind 100 Simulationen durchgeführt worden und das Histogramm jedes Parameters ist mit dem 'wahren' Parameter aufgetragen worden, markiert wie in Fig. 20 gezeigt.
Bei einem weiteren Test sind 1000 Realisierungen mit einer Größe von 64 · 64 Pixeln und 1000 Realisierungen mit einer Größe von 128 · 128 Pixeln erzeugt worden - ebenfalls nicht zu weit von der Nicht-Stationarität entfernt. Für jede dieser Realisierungen ist der MORSE-Schätzer angewandt worden, um eine Schätzung der Modellparameter zu erhalten. Es sind die empirischen 95%-Konfidenzintervalle für den Mittelwert jedes Parameters berechnet worden. In Fig. 24 ist ersichtlich, daß die bekannten Werte der Parameter sehr schön innerhalb diese Konfidenzintervalle fallen, was ferner anzeigt, daß der Schätzer zentral ist. Es sei ebenfalls bemerkt, daß die Breite der Konfidenzintervalle mit einem An stieg der Größe der künstlich erzeugten Texturproben abnimmt. Dies zeigt an, daß der Schätzer konsistent ist.
Daher ist ersichtlich, daß sehr zufriedenstellende Ergebnisse auch für das gemischte nicht-kausale SARMA-Modell erhalten werden.

Verteilungseigenschaften des Schätzers

Es könnte intuitiv erwartet werden, daß die Verteilung der geschätzten Parameter Gauß-verteilt ist, falls das antreibende Rauschen Gauß-verteilt war. Da dies jedoch formal schwer zu beweisen ist, ist es über Experimente wahrscheinlich gemacht worden. In Fig. 25 sind Histogramme der Parameter gezeigt, die aus den oben beschriebenen 1000 128 · 128 simulierten Texturen geschätzt worden sind. Die Hypothese, daß die gezeigten Histogramme Proben von Gauß'schen Schätzungen zeigen, wird nunmehr formal unter Anwendung des χ²- Tests geprüft. Für jeden Parameter werden die Schätzungen in K Klassen aufgeteilt, wobei K zu K - 1+3.3log&sub1;&sub0;(N) ausgewählt wird (wobei N die Anzahl an Beobachtungen ist). Damit der Test so stark wie möglich gemacht wird, werden die Klassen derart ausgewählt, daß die erwartete Anzahl an Beobachtungen in jeder Klasse gleich ist. Dies wird durchgeführt, indem die Trennungen als das (i·100/k)%-Quantil erzeugt werden (i = 1, ..., k-1). Die Teststatistik wird nunmehr berechnet zu:
wobei Obst die Anzahl an Beobachtungen in der i-ten Klasse ist. Diese Teststatistik ist χ²(K-3) verteilt (da sowohl der Mittelwert als auch die Standardabweichung geschätzt werden sollten). Die Tests werden in der unten stehenden Tabelle zusammengefaßt, aus der ersichtlich ist, daß die Verteilungshypothese für alle Parameter einfach angenommen wird.

Schlußfolgerung

Das dieser Schätzung zugrundeliegende Prinzip ist intuitiv sehr greifbar. Die für diesen Schätzer erhaltenen empirischen Ergebnisse sind tatsächlich sehr vielversprechend. Unter Anwendung dieses Schätzers können ebenfalls nichtkausale SARMA-Modelle betrachtet werden.

BEISPIEL 4: Test des z. Zt. bevorzugten Verfahrens

In dem vorliegenden Beispiel wurde ein Test durchgeführt, um die Anwendbarkeit des vorliegenden Verfahrens bei der Vorhersage der Knochenfestigkeit auf der Grundlage von Röntgenbildern des Knochens zu bewerten.

Datensatz

Der Datensatz bestand aus 20 Proben des dritten Lendenwirbels, der aus Leichenproben entnommen wurde, von denen jede zum Erhalt eines Röntgenbildes und einer KMD-Messung jedes Knochens einer Röntgenbestrahlung ausgesetzt wurde.
Nachfolgend wurde aus jedem Knochen eine vordere und eine hintere Kernprobe des trabekularen Knochens entnommen, und es wurde für jede dieser Proben die optimale Bruchbelastung bestimmt, indem ein Drei-Punkt-Testaufbau verwendet wurde, wie er per se bekannt ist.
Die Röntgenbilder wurden mit einer räumlichen Auflösung von 600 Linien pro cm abgetastet. Bei der endgültigen Auflösung wurde eine Fläche ausgewählt, die zentral zwischen den Positionen der ausgebohrten Proben angeordnet ist. Die Positionen der Kernproben können teilweise oder vollständig von der ausgewählten Fläche umfaßt sein.
Es wurde entschieden, als Festigkeit des Knochens die mittlere optimale Bruchbelastung der beiden Knochenproben eines Knochens auszuwählen. Dies entspricht weitestgehend der Festigkeit des Knochens in dem zentralen Abschnitt, wie es durch lineare Interpolation der Festigkeiten der beiden Proben bestimmt wird. Die Ergebnisse des vorliegenden Tests scheinen zu bestätigen, daß diese Entscheidung vernünftig ist.

Vorverarbeitung des Bildmaterials

Damit Bilder in der endgültigen, sehr viel geringeren Auflösung erhalten werden, die so gut wie möglich sind, hatten die abgetasteten Röntgenbilder eine große räumliche Auflösung (600 Linien pro Zentimeter).
Die vorliegenden Bilder waren so groß, daß eine Skalenverminderung unter Anwendung einer Burt-Adelson-Gauß-Pyramide unmöglich war, da der Computer zu wenig Speicher hatte. Somit wurde eine vorläufige Reduktion durchgeführt, indem das Bild unter Anwendung eines 3 · 3 Mittelwertfilters und durch nachfolgendes Verwerfen jedes zweiten Pixels in dem Bild geglättet wurde.
Nach dieser vorläufigen Reduktion hatte das Bild eine Größe, die für den Erhalt der endgültigen Auflösung von der Burt-Adelson-Gauß-Pyramide reduzierbar war (150 Linien pro Zentimeter).
Nach der Reduktion des Bildes wurde das Bild hintergrundkorrigiert, indem eine 31 · 31 median-gefilterte Version von diesem subtrahiert wurde.
Aus dem korrigierten Bild wurde eine in Frage kommende Region (FR) mit einer Größe von 256 · 256 Pixeln ausgewählt. Vorzugsweise wurde diese FR nahe dem Zentrum des Knochens ausgewählt, aufgrund der Tatsache, daß ein großer Gradient der Dicke der Trabekulae nahe dem vorderen und dem hinteren Ende des Knochens vorlag. Dieser Effekt und die Tatsache, daß der Knochen in Richtung seiner Enden abgerundet ist, erzeugte eine Anzahl nicht-linearer Effekte in den Bilddynamiken, die nicht auf triviale Weise behandelt werden konnten.
Nach Extrahieren der in Frage kommenden Region wurde das Bild gefiltert, um hochfrequentes Rauschen zu entfernen. Da das Signal-zu-Rausch-Verhältnis des vorliegenden Bildes sehr hoch war, konnte das hochfrequente Rauschen auf sehr 'behutsame' Weise entfernt werden:
fourtr > fft.$$
htof -CR M < fft.$$ powerpix -p 1,8 > mag.-$$
htof -CR P < fft.$$ > pha.$$
combine -pm pha.$$ < mag.$$ inv.fourtr htof -CR R skaliere
rm * .$$
Dieser Weg ist eine Alternative zu dem Verfahren des Algorithmus 1 auf den Seiten 30 und 31, da in diesem Verfahren T = 2 und p = 1,8 ist.
Aus diesem vorverarbeiteten Bild wird das projizierte trabekulare Muster (PTM) auf die selbe Weise extrahiert, wie oben in Verbindung mit Algorithmus 2 beschrieben wurde. Im wesentlichen wurden alle im Beispiel 1 beschriebenen Merkmale berechnet.

Statistische Verfahren

Es wurde versucht, ein allgemeines lineares Modell für die Vorhersage der endgültigen Festigkeit zu konstruieren. In diesem Zusammenhang ist es unbedingt erforderlich, die wesentlichen Parameter zu identifizieren und die unwesentlichen zu eliminieren. Da es eine große Anzahl an Variablen gibt und lediglich eine geringe Anzahl an Beobachtungen bewertet wurde, war es nicht möglich, ein Rückwärts- Eliminationsverfahren direkt anzuwenden, da (XTX) bei dieser Situation typischerweise singulär sein wird. Zusätzlich gibt es ebenfalls das Risiko, die falschen Variablen zu eliminieren. Daher wird nachfolgend ein Verfahren angewandt, das einem Vorwärts-Eliminationsverfahren ähnelt, das wie folgt durchgeführt wird:
1. Durchsuche die Liste möglicher Merkmale, bis das erste signifikante Merkmal identifiziert ist.
2. Fahre mit der Suche so lange fort, bis ein weiteres signifikantes Merkmal identifiziert und in das Modell einbezogen ist.
3. Prüfe, ob vorhergehende Merkmale eliminiert werden sollten.
Die Schritte 2. und 3. werden wiederholt, bis alle Merkmale untersucht worden sind. Natürlich sollte ein Signifikanzpegel für die beschriebenen Tests bestimmt werden. In der vorliegenden Untersuchung wird der Pegel zu p = 0,10 gewählt, da die Datenmenge einen sehr konservativen Test nicht möglich macht.
Es wurde deutlich, daß die optimale Bruchbelastung exponentiell auf die Parameter des Modells bezogen war, wobei das Modell als eine Regression mit dem Logarithmus der optima len Bruchbelastung als eine abhängige Variable durchgeführt wurde.
Unter Auswahl der Parameter gemäß dem obigen Prinzip wurden lediglich drei signifikante Parameter gefunden: â (siehe Algorithmus 4), Formel (23) quadriert und u²s (siehe Formel (3)).
Nachfolgend werden diese Parameter in zwei Modellen bewertet, von denen eines das Alter als eine Variable verwendet und das andere das Alter nicht verwendet.

Modellbildung ohne Alter

Es wurde ein allgemeines lineares Modell aufgestellt, das das Alter als eine beschreibende Variable nicht verwendet. Dieses Modell versucht daher, die optimale Bruchbelastung auf der Grundlage der extrahierten Strukturindizes vorherzusagen. Das Ergebnis ist in dem nachfolgenden Programmausdruck zusammengefaßt:
Reststandardfehler: 0,2901 auf 15 Freiheitsgraden
Mehrfach-R-quadriert: 0,8158
F-Statistik: 22,14 auf 3 und 15 Freiheitsgraden, der p-Wert ist 9,17e-06
Korrelation der Koeffizienten:
Erklärendes Modell:
Favmax = exp(-70,8473·(Formel23)²-0,0287·â +143625.7515·u²s+11.2930)
In Fig. 26 ist ein Plot der logarithmischen optimalen Bruchbelastung gegen die vorhergesagte Festigkeit dargestellt. In dieser Figur sind die Meßpunkte durch ihre Anzahl in der Datenbank dargestellt.
An dieser Stelle ist es natürlich interessant zu wissen, ob die an allen Proben durchgeführten KMD-Messungen als ein statistisch signifikantes Merkmal teilnehmen, das zu dem Anstieg des Mehrfach-R-quadriert beiträgt. Das Ergebnis ist auch in dem untenstehenden Ausdruck ersichtlich.
Reste bzw. Residuen:
Reststandardfehler: 0,3002 auf 14 Freiheitsgraden
Mehrfach-R-quadriert: 0,8158
F-Statistik: 15,5 auf 4 und 14 Freiheitsgraden, der p-Wert ist 8,43e-05
Korrelation der Koeffizienten:
Es ist offensichtlich, daß die KMD kein signifikantes Merkmal ist und daß sie aus dem Modell eliminiert werden kann. Unter Verwendung der KMD als einziger beschreibender Parameter wird ein Mehrfach-R-quadriert von R² = 0,2142 erhalten, der anzeigt, daß die KMD ein schlechter Indikator für die Qualität des trabekularen Knochens ist. Somit werden die Aspekte des sich auf den trabekularen Knochen beziehenden Frakturrisikos schlecht von der KMD beschrieben.

Modellbildung mit Alter

Unter Verwendung des Alters als beschreibender Parameter wird das folgende Modell erhalten:
Reste bzw. Residuen:
Reststandardfehler: 0,2533 auf 14 Freiheitsgraden
Mehrfach-R-quadriert: 0,8689
F-Statistik: 23,2 auf 4 und 14 Freiheitsgraden, der p-Wert ist 4,707e-06
Korrelation der Koeffizienten:
Erklärendes Modell:
Favmax = exp(-68,3085·(Formel23)²-0,0236·â +152951.1278·u²s-0.016·Alter+11.8882)
In Fig. 27 ist ein Plot der logarithmischen optimalen Bruchbelastung gegen die vorhergesagte Festigkeit dargestellt. In dieser Figur sind die Meßpunkte durch ihre Anzahl in der Datenbank dargestellt.

Diskussion

Aus dem Ergebnis des vorliegenden Beispiels wird deutlich, daß die optimale Bruchbelastung des trabekularen Knochens lediglich unter Anwendung von Texturparametern vorhergesagt werden kann - und dies sogar mit einer hohen Genauigkeit.
Drei Bildmerkmale werden verwendet: die Formel 23 quadriert, u²s und â, welche jeweils die Dichte der horizontalen Trabekulae (die dicksten), den Zusammenbruch der trabekularen Struktur und die globale trabekulare Dichte beschreiben. Außerdem wird ersichtlich, daß weder das Geschlecht noch die KMD Information zu der Beschreibung hinzufügen können, die unter Verwendung der obigen drei Bildmerkmale erhalten wird.
Es ist jedoch ersichtlich, daß das Alter ein signifikanter Parameter ist. Das Einbeziehen des Alters in das Modell erhöht den Mehrfach-R-quadriert auf R² = 0,8689. Es sei bemerkt, daß die beiden Extrembeobachtungen (258 und 267) zu dem Wert von R² stark beitragen. Ein Weglassen dieser Beobachtungen ergibt ein R² = 0,7538, das Alter einbezogen, das weiterhin als ein hoher Mehrfach-R-quadriert angesehen werden sollte.
Wie jedoch im Zusammenhang mit dem Algorithmus 1 beschrieben wurde, haben die beim Wiederherstellen des Bildes verwendeten Parameter eine starke Auswirkung auf die Anwendbarkeit der Merkmale, die aus dem resultierenden Bild extrahiert werden. In der Tat gibt das Merkmal σM allein einen R² von 0,92, wenn im Algorithmus 1 p = 1,4 und T = 4 anstatt p = T = 2 verwendet wird.
In Fig. 28 ist die Korrelation zwischen σM und dem Logarithmus der optimalen Bruchbelastung für einen Test dargestellt, basierend auf dem Bild und der optimalen Bruchbelastung - oben beschriebenes Material. Selbst wenn der äußerste Punkt in der unteren rechten Ecke der Figur entfernt wird, wird eine Korrelation von 0,88 erhalten.
Somit kann ein großer Unterschied gesehen werden, abhängig von der tatsächlichen Weise, wie das PTM des Bildes erhalten wird.

Claims

1. Verfahren zum Schätzen der Knochenqualität eines Vertebraten auf der Basis zweidimensionaler Bilddaten, die auf die trabekulare Struktur wenigstens eines Knochenteils des Vertebraten bezogene Information umfaßt, wobei die Bilddaten Daten sind, die dadurch erhalten werden, daß wenigstens der Knochenteil elektromagnetischer Strahlung ausgesetzt wird, wobei das Verfahren das Unterziehen der Bilddaten einer statistischen Analyse umfaßt, die folgendes umfaßt:

- einen Hintergrundskorrekturschritt, der wenigstens das Vermindern oder optional das Entfernen niederfrequenter Information mit einer Frequenz umfaßt, die wesentlich niedriger als der Abstand der projizierten Trabekulae ist,

- einen Bildmanipulations- und Merkmalsextraktionsschritt, bei dem wenigstens die Information über die lokale Bildintensität sowie die Variation der lokalen Intensität verwendet wird, um auf die trabekulare Struktur des Knochenteils bezogene Information zu extrahieren, wobei der Bildmanipulations- und Merkmalextraktionsprozeß folgendes umfaßt: Ermitteln einer parametrischen Schätzung des zweidimensionalen Leistungsspektrums der Bilddaten unter Verwendung eines Modells, das aus der Gruppe ausgewählt ist, die besteht aus: kausalen simultanen autoregressiven Modellen mit gleitendem Mittelwert (SARMA), in denen das Spektrum entweder unter Verwendung von gewichteten Viertelebenen-Modellen oder von Modellen mit nicht-symmetrischen Halbebenen bestimmt wird; und nicht-kausalen SARMA-Modellen, in denen das Verfahren zum Schätzen des nicht kausalen simultanen auto-regressiven Modells mit gleitendem Mittelwert in zwei oder mehreren Dimensionen die Optimierung eines gegebenen Direktmaßes der Flachheit des Restspektrums des Modells und das Extrahieren von auf die Anisotropie des projizierten trabekularen Musters der Bilddaten bezogenen Merkmale aus der parametrischen Schätzung umfaßt,

- und einen Knochenqualitäts-Schätzschritt, in dem die Knochenqualität des Vertebraten auf der Grundlage der abgeleiteten Merkmale und optional anderer auf den Knochen oder den Vertebraten bezogener Merkmale und einer vorbestimmten Beziehung zwischen den Merkmalen und Parametern der Referenzknochenqualität geschätzt wird.

2. Verfahren nach Anspruch 1, bei welchem die elektromagnetische Strahlung Röntgenstrahlung ist.

3. Verfahren nach einem der vorhergehenden Ansprüche, bei welchem die Bilddaten von einem Röntgenfilm abgetastet werden, wie einem Röntgenfilm mit einer Auflösung von wenigstens 4 Linienpaaren pro Zentimeter, wie wenigstens 5 Linienpaaren pro Zentimeter, vorzugsweise wenigstens 10 Linienpaaren pro Zentimeter, wie wenigstens 25 Linienpaaren pro Zentimeter, vorzugsweise wenigstens 50 Linienpaaren pro Zentimeter, wie wenigstens 100 Linienpaaren pro Zentimeter, wie wenigstens 250 Linienpaaren pro Zentimeter, vorzugsweise wenigstens 500 Linienpaaren pro Zentimeter, wie wenigstens 600 Linienpaaren pro Zentimeter.

4. Verfahren nach Anspruch 3, bei welchem die Abtastung mit einer Auflösung von wenigstens 10 Linien pro Zentimeter durchgeführt worden ist, wie wenigstens 25 Linien pro Zentimeter, vorzugsweise wenigstens 100 Linien pro Zentimeter, wie wenigstens 200 Linien pro Zentimeter, wie wenigstens 250 Linien pro Zentimeter.

5. Verfahren nach Anspruch 1, bei welchem Information mit Frequenzen, die verglichen mit dem Abstand der projizierten Trabekulae halb so groß oder weniger, wie ein Viertel so groß oder weniger, vorzugsweise ein Zehntel so groß oder weniger, wenigstens vermindert oder optional entfernt wird.

6. Verfahren nach Anspruch 5, bei welchem der Hintergrund-Korrekturprozeß das Erzeugen sekundärer Bilddaten als Ergebnis einer Durchführung eines Median- Filterns mit einer vorbestimmten Kerngröße und einer Subtraktion dieses Ergebnisses von den ursprünglichen Bilddaten umfaßt.

7. Verfahren nach Anspruch 5, bei welchem der Hintergrund-Korrekturprozeß das Erzeugen sekundärer Bilddaten als Ergebnis einer Durchführung eines Mittelwert-Filterns mit einer vorbestimmten Kerngröße und einer Subtraktion dieses Ergebnisses von den ursprünglichen Bilddaten umfaßt.

8. Verfahren nach Anspruch 5, bei welchem der Hintergrund-Korrekturprozeß das globale Anpassen eines zweidimensionalen Polynoms an die Bilddaten und das Erzeugen hintergrundkorrigierter Bilddaten auf der Basis der Reste des Anpaßprozesses umfaßt.

9. Verfahren nach Anspruch 8, bei welchem die Ordnung des Polynoms höchstens 15, wie höchstens 10, vorzugsweise höchstens 5 ist.

10. Verfahren nach Anspruch 6 oder 7, bei welchem die Kerngröße höchstens 1/2 der Bilddaten, wie höchstens 1/4 der Bilddaten, vorzugsweise höchstens 1/10 der Bilddaten, noch bevorzugter höchstens 1/20 der Bilddaten ist.

11. Verfahren nach Anspruch 1, bei welchem der Bildmanipulations- und Merkmalsextraktionsschritt zusätzlich umfaßt:

Ermitteln einer Schätzung des projizierten trabekularen Musters der Bilddaten, indem

ein Fourier-Verfahren angewandt wird, welches das Unterziehen der Bilddaten einer Fourier-Transformation, das Durchführen einer nachfolgenden mathematischen Transformation der Fourier-transformierten Daten und das Rückkonvertieren der transformierten Daten in den Ortsbereich umfaßt, um die sich auf die trabekulare Struktur beziehende, vorherrschende Information in den Bilddaten hervorzuheben,

die Information, die in den auf die trabekulare Struktur bezogenen Bilddaten hervorgehoben worden ist, einer morphologischen Graustufenoperation unterzogen wird, und

eine Operation durchgeführt wird, welche das Bestimmen eines Maßes für jeden Hintergrundpixel, basierend auf einer Bestimmung des Abstandes von dem Hintergrundpixel in einer Anzahl vorgegebener Richtungen in den Bilddaten zu dem nächsten Vordergrundpixel, und das Extrahieren wenigstens eines Merkmals, das sich auf die trabekulare Struktur des geschätzten projizierten trabekularen Musters bezieht, aus dem Ergebnis der Operation umfaßt,

wobei der Knochenqualitäts-Schätzprozeß ausgelegt ist, die Knochenqualität auch auf der Grundlage des wenigstens einen so abgeleiteten Merkmals zu schätzen.

12. Verfahren nach Anspruch 1, bei welchem die Ordnung des Modells identifiziert wird, die Parameter des Modells geschätzt werden und das Spektrum des geschätzten Modells erzeugt wird.

13. Verfahren nach Anspruch 12, bei welchem das Verfahren, auf dessen Grundlage die parametrische Schätzung ermittelt wird, ein nicht-kausales SARMA-Modell ist, und bei welchem das Verfahren zum Schätzen des nicht-kausalen simultanen auto-regressiven Modells mit gleitendem Mittelwert in zwei oder mehreren Dimensionen das Optimieren eines gegebenen Direktmaßes der Flachheit des Restspektrums des Modells umfaßt.

14. Verfahren nach einem der vorhergehenden Ansprüche, bei welchem wenigstens ein Merkmal auf Parameter einer Kontur bezogen ist, welche wenigstens einen vorbestimmten Prozentsatz der Energie des Leistungsspektrums umgibt.

15. Verfahren nach Anspruch 14, bei welchem die Kontur bestimmt wird, indem

(a) eine Kontur in Ausdrücken einer Achse in jeder Dimension um das Zentrum des Leistungsspektrums definiert wird, wobei alle Punkte auf der Kontur im wesentlichen den gleichen Abstand zum Zentrum des Leistungsspektrums haben, und die Kontur weniger als den vorbestimmten Prozentsatz der Energie des Leistungsspektrums umgibt,

(b) für jede Dimension der Daten der Prozentsatz der Energie des Leistungsspektrums berechnet wird, das von einer gedehnten Kontur umgeben ist, in der die Achse in der in Frage kommenden Dimension um einen vorbestimmten Abstand verlängert ist,

(c) die Achse der Kontur in der Dimension, in der die Vergrößerung der von der gedehnten Kontur umgebenen Energie am größten ist, um den vorbestimmten Abstand verlängert wird, und

(d) die Schritte (b) und (c) wiederholt werden, bis der von der Kontur umgebene Prozentsatz den vor bestimmten Prozentsatz überschreitet oder angleicht.

16. Verfahren nach Anspruch 15, bei welchem die Achsen der Kontur orthogonal und entlang von Hauptrichtungen des Leistungsspektrums definiert sind.

17. Verfahren nach Anspruch 11, bei welchem die nachfolgende mathematische Transformation nicht-linear ist und durch Erhöhen der Daten mit einer Potenz größer als 1, wie einer Potenz in dem Bereich von 1,1-10, vorzugsweise 1,5-4, wie ungefähr 2.0 durchgeführt wird.

18. Verfahren nach Anspruch 1 oder 17, bei welchem die Prozedur (c) angewandt wird und bei welchem eine dritte mathematische Transformation durchgeführt wird, um lediglich einen Teil der Betragsinformation der transformierten Daten beizubehalten.

19. Verfahren nach einem der Ansprüche 11, 17 oder 18, bei welchem die Information, die in den auf die trabekulare Struktur bezogenen Bilddaten hervorgehoben worden sind, einer morphologischen Graustufen- Operation unterzogen wird, wie einer top-hat- Operation, bei der allen in den top-hat passenden Pixeln ein Wert zugeordnet wird, der von einem angegebenen Wert für einen Hintergrundpixel unterschiedlich ist.

20. Verfahren nach Anspruch 11, bei welchem die Bilddaten, die der top-hat-Operation unterzogen worden sind, mit einem Schwellwert verglichen werden, um denjenigen Pixeln einen Wert für einen Vordergrundpixel zuzuordnen, die einen Pixelwert haben, der von dem Wert für einen Hintergrundpixel unterschiedlich ist.

21. Verfahren nach Anspruch 11 oder 20, bei welchem die Bilddaten, die der morphologischen Graustufen- Operation unterzogen worden sind, einer weiteren Rauschverminderung unterzogen werden, um isolierte Einzelpixel und Pixelgruppen, die kleiner als eine vorbestimmte Schwellwertanzahl an Pixeln sind, im wesentlichen zu entfernen, die einen Pixelwert haben, der von dem Wert für einen Hintergrundpixel unterschiedlich ist.

22. Verfahren nach einem der Ansprüche 19 bis 21, bei welchem eine Operation durchgeführt wird, welche das Bestimmen, für jeden Hintergrundpixel, des Abstands gemäß einer vorgegebenen Metrik von dem Hintergrundpixel zu dem nächsten Vordergrundpixel umfaßt, und bei welchem wenigstens ein Merkmal aus dem Ergebnis der Operation extrahiert wird.

23. Verfahren nach Anspruch 22, bei welchem sich wenigstens ein Merkmal auf einen Mittelwert und/oder eine Standardabweichung und/oder den Variationskoeffizienten und/oder den Schiefe-Koeffizienten und/oder den Kurtosis-Koeffizienten der bestimmten Abstände von allen Hintergrundpixeln zu dem nächsten Vordergrundpixel bezieht.

24. Verfahren nach Anspruch 11, bei welchem sich wenigstens ein Merkmal auf einen Mittelwert und/oder eine Standardabweichung und/oder den Variationskoeffizienten und/oder den Schiefe-Koeffizienten und/oder den Kurtosis-Koeffizienten der bestimmten Maße für alle Hintergrundpixel bezieht und/oder bei welchem das Maß für jeden Hintergrundpixel ein Mittelwert und/oder eine Standardabweichung und/oder der Variationskoeffizient und/oder der Schiefe-Koeffizienten und/oder der Kurtosis-Koeffizienten und/oder der maximale Abstand der bestimmten Abstände in die gegebenen Richtungen ist.

25. Verfahren nach einem der vorhergehenden Ansprüche, bei welchem wenigstens ein Merkmal, das sich auf den Knochen oder den Vertebraten bezieht, wie das Alter und/oder Geschlecht und/oder die Spezies und/oder die Rasse und/oder der angegebene in dem Vertebraten beobachtete Knochen und/oder eine geschätzte Knochenmineraldichte des Knochens und/oder ein geschätzter Knochenmineralgehalt des Knochens, in den Schätzungsprozeß einbezogen wird.

26. Verfahren nach einem der vorhergehenden Ansprüche, bei welchem der Schätzungsprozeß auf einem statistischen Modell beruht, das die Korrelationsstruktur in dem Datensatz berücksichtigt, um den signifikanten Merkmalen geeignete Gewichtungen gemäß der vorbestimmten Beziehung zuzuordnen.

27. Verfahren nach einem der vorhergehenden Ansprüche, bei welchem die Referenzknochenqualität eine absolute oder relative Knochenqualität von Knochen ist, für die Bilddaten gemäß einem der vorhergehenden Ansprüche verarbeitet worden sind.

28. Verfahren nach Anspruch 27, bei welchem die Referenzknochenqualität auf der Grundlage einer Messung der mechanischen Knochenfestigkeit und/oder einer Messung der Knochenmineraldichte und/oder auf der Grundlage einer Messung des Knochenmineralgehalts und/oder auf der Grundlage einer Bewertung der Knochenfestigkeit durch einen erfahrenen Radiologen bestimmt wird.

29. Verfahren nach einem der vorhergehenden Ansprüche, bei welchem die vorbestimmte Beziehung zwischen den Merkmalen und den Parametern der Referenzknochenqualität auf der Grundlage von Parametern der Referenzknochenqualität und Merkmalen, die gemäß einem der vorhergehenden Ansprüche extrahiert werden, aufgestellt wird.

30. Verfahren nach Anspruch 29, bei welchem die vorbestimmte Beziehung in Ausdrücken eines Modells definiert ist, das aus der Gruppe ausgewählt ist, die besteht aus: einem allgemeinen linearen Modell, einem verallgemeinerten linearen Modell, einem künstlichen neuronalen Netzwerk, einem kausalen probabilistischen Netz und Klassifikations- und Regressions-Bäumen.

31. Verfahren nach einem der vorhergehenden Ansprüche, bei welchem der Vertebrat ein menschliches Wesen, ein Pferd, ein großer Affe, ein Menschenaffe, ein Schwein oder eine Kuh ist.

32. Verfahren nach einem der vorhergehenden Ansprüche, bei welchem der Knochen aus der Gruppe ausgewählt ist, die besteht aus dem Radius, dem Femur, dem Corpus Vertebrae (L1, L2, L3, L4, L5, T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11, T12, C1, C2, C3, C4, C5, C6, C7), dem Calcaneus, dem Talus, dem Os Carpi, dem Metatars, dem Metacarpi, dem Falanges, dem Tibia, dem Fibula, dem Patella, dem Ulna, dem Humerus, dem Mandible, dem Clavicula, dem Scapula, dem Os Coxae, dem Os Naviculare, dem Os Cuboideum, dem Os Cuneiform I, dem Os Cuneiform II, dem Os Cuneiform III, dem Os Sacrum, dem Os Coccygis.

33. Verfahren nach Anspruch 1, bei welchem der Bildmanipulations- und Merkmalsextraktionsprozeß ferner das Unterziehen der Bilddaten dem Prozeß zum Ermitteln, auf der Grundlage von Bilddaten, an denen ein Fourier-Verfahren zum Hervorheben der sich auf die trabekulare Struktur beziehenden, vorherrschenden Information in den Bilddaten angewandt worden ist, einer Schätzung einer Graustufen-Miteintrittswahrscheinlichkeit-Matrix und das Extrahieren wenigstens eines Merkmals auf der Grundlage der bestimmten Miteintrittswahrscheinlichkeit-Matrix umfaßt.

34. Verfahren nach Anspruch 1 oder 2, bei welchem der Bildmanipulations- und Merkmalsextraktionsprozeß ferner das Unterziehen der Bilddaten dem Prozeß zum Ermitteln, auf der Grundlage einer Frequenzanalyse der parametrischen Schätzung des Leistungsspektrums der Bilddaten, von sich auf die Periodizität der trabekularen Struktur des Knochenteils beziehenden Merkmalen umfaßt.

35. Verfahren nach Anspruch 34, bei welchem die sich auf die Periodizität der trabekularen Struktur des Knochenteils beziehenden Merkmale aus einer Information abgeleitet werden, die sich auf die Höhe, die Breite und/oder die Krümmung von Peaks und/oder die Gesamtfläche unter Peaks des Leistungsspektrums der Bilddaten, das geglättete Periodogramm der Bilddaten oder die Schätzung der parametrischen Schätzung des Leistungsspektrums der Bilddaten bezieht.

36. Verfahren nach Anspruch 35, bei welchem die Information abgeleitet wird, die sich auf die Höhe, die Breite und/oder die Krümmung von Peaks des geglätteten Periodogramms der Bilddaten und/oder die Gesamtfläche unter deren Peaks bezieht, und bei welchem die Positionen der Peaks in dem Leistungsspektrum der Bilddaten identifiziert werden.

37. Verfahren nach Anspruch 1, bei welchem die Knochenqualität des Vertebraten mit einem Mehrfachkorrelationskoeffizienten besser als 0,5, wie besser als 0,55, wie besser als 0,6, vorzugsweise besser als 0,65, wie besser als 0,7, vorzugsweise besser als 0,8, wie besser als 0,85 geschätzt wird.

38. Verfahren nach Anspruch 13, bei welchem die Optimierung des vorgegebenen Flachheitsmaßes erhalten wird, indem

(a) ein Satz Anfangsparameter für das Modell erzeugt wird,

(b) das Restspektrum des Modells auf der Grundlage der Parameter erzeugt wird,

(c) das Maß der Flachheit des Restspektrums ermittelt wird,

(d) eine neue Iteration der Parameter auf der Grundlage des Flachheitsmaßes und einer Suchrichtung in dem Parameterraum ermittelt wird,

(e) die Schritte (b)-(d) wiederholt werden, bis ein vorgegebenes Abbruchkriterium erreicht ist.

39. Verfahren nach Anspruch 38, bei welchem das Flachheitsmaß des Restspektrums als die Entropie des normalisierten Rest-Leistungsspektrums erhalten wird, unter Nichtbeachtung des Wertes bei der Frequenz Null.