DE102009022834A1

DE102009022834A1 - Verfahren zur automatischen Analyse von Bilddaten einer Struktur

Info

Publication number: DE102009022834A1
Application number: DE102009022834A
Authority: DE
Inventors: Mausumi Acharyya
Original assignee: Siemens AG; Siemens Corp
Current assignee: Siemens AG; Siemens Corp
Priority date: 2009-05-27
Filing date: 2009-05-27
Publication date: 2010-12-09
Also published as: US20100303358A1

Abstract

Es wird ein Verfahren zur automatischen Analyse von Bilddaten (D) einer Struktur beschrieben, welches Verfahren folgende Verfahrensschritte umfasst: Bereitstellung von Bilddaten (D) in Form eines dreidimensionalen Voxel-Arrays, Durchführung einer Segmentierung des Voxel-Arrays zur Ermittlung einer Voxel-Teilmenge, Durchführung einer Merkmalsextraktion zumindest für bestimmte Voxel der Voxel-Teilmenge zur Erzeugung einer Merkmalskarte (FM), Erzeugung einer skalaren Differenzkarte auf Basis der Merkmalskarte (FM), Durchführung einer Klassifizierung anhand der skalaren Differenzkarte und Identifizieren einer Strukturanomalie ind den Bilddaten (D) auf Basis eines Klassifikationsergebnisses. Weiterhin werden ein Verfahren zur Ansteuerung einer Bildanzeigeeinrichtung, ein Bildverarbeitungssystem (1) und ein bildgebendes System (1) beschrieben.

Description

Die Erfindung betrifft ein Verfahren zur Durchführung einer automatischen Analyse von Bilddaten einer Struktur. Die Erfindung betrifft weiterhin ein Bildverarbeitungssystem und ein Computerprogrammprodukt, um die Schritte nach diesem Verfahren auszuführen.
In heutigen Diagnoseverfahren werden Computertomographie (CT) oder Magnetresonanzbildgebung (MRI) häufig benutzt, um die Ursachen von Krankheiten zu identifizieren, die mit anderen diagnostischen Methoden nicht einfach erkannt werden können. MRI wird dabei zur Bildgebung von Weichteilgewebe wie inneren Organen bevorzugt, wogegen CT besonders zur Bildgebung von hartem Gewebe wie z. B. Knochen geeignet ist. In den bekannten CT-Scan-Verfahren liefert der Computertomograph beispielsweise einen Datenstrom bzw. ein Signal, welches unterschiedliche radiographische Intensitäten (auch ”Massenbelegungsdichten” oder ”Abschwächungen” genannt) in einer Region des untersuchten Patienten repräsentiert. Die Signale werden meist mit Hilfe einer um eine Rotationsachse des Patienten rotierenden Röntgenquelle und einer der Röntgenquelle gegenüberliegenden Detektorfläche erzeugt. Während der Rotation werden Röntgenquelle und Detektorfläche sequenziell oder kontinuierlich entlang der Rotationsachse verfahren. Aus den Signalen werden dann eine Serie von Querschnitt-Bildschichten oder in sonstiger Weise Volumenbilddaten erzeugt. Ein Punkt oder Pixel des Bildes wird dann entsprechend der an dem betreffenden Punkt gemessenen Abschwächung (üblicherweise in Form von Hounsfield-Werten ausgedrückt) ausgegeben. Die Abschwächungen der verschiedenen Gewebe können verwendet werden, um beliebige Schnittbilder graphisch z. B. in Form von Graustufenbildern auszugeben.
Unter Nutzung eines CT- oder MRI-Verfahrens kann also ein dreidimensionales Array von Voxeln erstellt werden, das automatisch analysiert und zur Ausgabe realistischer dreidimensionaler Darstellungen (3D-Darstellungen) für diagnostische Zwecke verwendet werden kann, z. B. 3D-Darstellung des Herzens, des Darms, eines Knochens etc. Dies macht diese Techniken zum einen wertvoll nur zur Analyse von Weichteilgewebe, aber auch zur Erkennung von Brüchen, die mit herkömmlichen Röntgenstrahlentechniken nicht detektiert werden können, da diese auf die Erkennung von Brüchen beschränkt sind, die eine klar sichtbare Fehlausrichtung des Knochens zeigen, wie einen kompletten oder einen zusammengestauchten Bruch. Röntgenbilder sind im Allgemeinen ungeeignet zur Erkennung von schwächeren Haarrissen oder Spannungsbrüchen.
Haarrissbrüche in Knochen des Körpers können aus verschiedenen Gründen auftreten. Typischerweise werden Haarrissbrüche in gewichttragenden Knochen, wie z. B. Fußknochen oder Beinknochen, im Allgemeinen durch wiederholte Spannungsbelastung, z. B. beim Sport, verursacht, und werden daher allgemein als Spannungsbrüche bezeichnet. Im Gegensatz dazu treten Haarrissbrüche in Schädel- und Gesichtsknochen meist aufgrund von traumatischen Verletzungen auf. In einem Spannungs- oder Haarrissbruch sind die gebrochenen Knochenstücke nicht sichtbar fehlausgerichtet.
Eine Knochenbruchidentifizierung mit Hilfe von dreidimensionalen Bilddaten umfasst in erster Linie einen ersten ”Segmentierungsschritt”, in dem mögliche Knochenregionen identifiziert werden, sowie nachfolgende Bearbeitungsschritte auf den identifizierten Regionen. Nach dem Stand der Technik wird die Segmentierung üblicherweise durchgeführt, indem zunächst ein Schwellenwertverfahren (z. B. durch Vergleich der vorhandenen Voxelwerte mit einem bekannten Schwellenabschwächungswert für Knochen und Verwerfen aller Voxel mit einem Wert unter diesem Schwellenwert) und dann ein ”Nachbesserungsverfahren” durchgeführt wird, wie z. B. ein Konnektivitätsverfahren und/oder eine manuelle Nachbearbeitung. Eine Segmentierung von größe ren Knochen, wie z. B. gewichttragenden Knochen, mit einem Schwellenwertverfahrens liefert üblicherweise akzeptable Ergebnisse, da die Hounsfield-Abschwächungswerte für Knochen im Allgemeinen höher sind als die des umgebenden Weichteilgewebes, zudem gut bekannt sind und im Wesentlichen zwischen verschiedenen Aufnahmen auch konstant sind.
Ungünstigerweise ist eine solche Diskriminierung, basierend auf einem direkten Schwellenwertverfahren, in der Praxis nicht für alle Knochenbereiche gut geeignet, insbesondere nicht für dünne, nicht gewichttragende Knochen, wie die Knochen innerhalb des Kopfes. Daher bleibt die Erkennung von Haarrissbrüchen in Schädel- oder Gesichtsschädelknochen auch bei Anwendung von CT- oder MRI-Technologien schwierig. Der Hauptgrund liegt darin, dass das akquirierte Signal ein digital gesampeltes Signal mit einem inhärenten Verlust an Informationen (im Vergleich zu einem analogen Signal) ist. Weiterhin wird durch das Rauschen in den gewonnenen Bildern die Erkennung und nachfolgende Visualisierung von Haarrissbrüchen extrem erschwert. Herkömmliche Techniken wie Kantendetektion oder Kennzeichnung von Punkten von Oberflächendiskontinuitäten scheitern häufig. Insbesondere im Falle von dünnen Knochen wie den Gesichtsknochen funktioniert die herkömmliche Technik zur Detektion von Punkten an Oberflächendiskontinuitäten (typischerweise basierend auf Kantendetektoren wie dem Harris-Detektor) aufgrund der ausgeprägten Intensitätsinhomogenitäten und des Rauschens in den CT-Bildern nicht gut. Zudem kann nach dem Schwellenwertverfahren ein sog. ”Partialvolumeneffekt” auftreten. Dieser Effekt tritt auf, wenn benachbarte unterschiedliche Gewebetypen – z. B. dichte Knorpel und dünne, kompakte Knochen – zu einem gleichen Abschwächungswert in einem Voxel führen. Da herkömmliche Segmentierungsverfahren nur auf einer Graustufenwert-Diskriminierung basieren, führt der Partialvolumeneffekt im Ergebnis zu Schwierigkeiten bei der Entdeckung von dünnen Knochenstrukturen. Ein anderes Problem bei Schwellenwertverfahren wird durch Artefakte verursacht, welche durch abrupte Übergänge zwischen Geweben von hoher und niedriger Dichte auch dann auftreten können, wenn das Signal eng genug gesampelt wurde. Weiterhin müssen unbekannte lokale Abschwächungen im Bild auch berücksichtigt werden. Schließlich macht der hohe Grad an Unterschieden von Knochen unterschiedlicher Subjekte eine automatische Analyse schwierig.
Es ist daher eine Aufgabe der Erfindung, verbesserte Mittel zur sicheren Erkennung von Haarrissbrüchen in dreidimensionalen Bilddaten, insbesondere CT-Bilddaten, anzugeben.
Diese Aufgabe wird mit einem Verfahren gemäß Patentanspruch 1, durch ein Bildverarbeitungssystem gemäß Patentanspruch 14 und durch ein Computerprogrammprodukt gemäß Patentanspruch 15 gelöst.
Bei dem erfindungsgemäßen Verfahren zur Durchführung einer automatischen Analyse von Bilddaten einer Struktur werden zunächst Bilddaten in Form eines dreidimensionalen Voxel-Arrays bereitgestellt, und es wird eine adaptive Segmentierung des Voxel-Arrays durchgeführt, um eine Voxel-Teilmenge zu ermitteln. Nachfolgend wird zumindest für bestimmte Voxel der Voxel-Teilmenge eine Merkmalsextraktion durchgeführt, um eine Merkmalskarte zu erzeugen. Diese Merkmalskarte wird dann als Basis zur Erzeugung einer skalaren Differenzkarte verwendet. Dann wird eine Klassifizierung anhand der skalaren Differenzkarte durchgeführt, und es wird eine Strukturanomalie in den Bilddaten auf Basis eines Klassifikationsergebnisses identifiziert.
Handelt es sich bei den Bilddaten um CT-Bilddaten, kann das dreidimensionale Voxel-Array beispielsweise eine dreidimensionale Matrix von Hounsfield-Graustufenwerten für die Punkte in einer dreidimensionalen interessierenden Region des untersuchten Patienten sein. Jeder Punkt oder Voxel in dem dreidimensionalen CT-Bild kann dabei entsprechend der an diesem Punkt gemessenen Massenbelegungsdichte oder Abschwächung mit einem bestimmten Graustufenwert verknüpft werden.
Der Segmentierungsschritt kann die Anwendung von Computertechniken zur Identifizierung von Regionen des dreidimensionalen CT-Bilds umfassen, die höchstwahrscheinlich Knochen entsprechen, auch wenn der Knochen dünn, beispielsweise ein Gesichtsknochen, ist.
Das erfindungsgemäße Verfahren ist folglich vorteilhaft, da eine rechenzeitintensive Durchführung der Merkmalsanalyse nur in den Regionen des Voxel-Arrays durchgeführt werden muss, welche die spezifischen relevanten Voxel enthalten, die höchstwahrscheinlich Knochen entsprechen.
Die dabei als Ergebnis erhaltene Merkmalsmatrix charakterisiert in Form von Merkmalen, die typisch für die spezifischen Knochenarten sind, die Knochenvoxel des CT-Bilds besser und bietet daher leistungsfähigere Eingangsdaten für die nachfolgende Klassifizierung, in welcher die Voxel des Voxel-Arrays auf Basis ihrer entsprechenden Merkmale als gesund oder gebrochen markiert werden können. Das erfindungsgemäße Verfahren bietet so einen verbesserten Weg zur Identifizierung von feinen Anomalien in Strukturen, wie Haarrissbrüchen in Knochenregionen, und kann daher insbesondere bei der Diagnose von Brüchen helfen, die bisher als ”schwierig” einzustufen waren, da sie mit den bekannten Verarbeitungstechniken im Allgemeinen nicht identifizierbar sind.
In einem erfindungsgemäßen Verfahren zur Kontrolle einer Bildausgabeeinrichtung zur Ausgabe von Bildern einer Struktur wird aus den dreidimensionalen Bilddaten der Struktur ein Bild der Struktur abgeleitet und die Spur einer Strukturanomalie, welche mit Hilfe der oben beschriebenen Methode zur Durchführung einer automatischen Analyse von Bilddaten einer Struktur identifiziert wurde, wird dem Bild graphisch überlagert.
Ein erfindungsgemäßes Bildverarbeitungssystem zur Durchführung einer automatischen Analyse von Bilddaten einer Struktur umfasst zum einen eine Bilddatenquelle zur Bereitstellung von Bilddaten in Form eines dreidimensionalen Voxel-Arrays. Diese Bilddatenquelle kann ein Gerät wie ein Computertomograph sein, welches direkt solche Daten liefert, oder ein Speicher, in dem die Bilddaten zuvor gespeichert wurden. Das Bildverarbeitungssystem umfasst weiterhin ein Bildanalysesystem, welches so ausgebildet ist, dass eine Segmentierung des Voxel-Arrays zur Ermittlung einer Voxel-Teilmenge und eine Merkmalsextraktion zumindest für bestimmte Voxel der Voxel-Teilmenge zur Erzeugung einer Merkmalskarte durchgeführt werden, eine skalare Differenzkarte auf Basis der Merkmalskarte erzeugt wird, eine Klassifizierung anhand der skalaren Differenzkarte durchgeführt wird und eine Strukturanomalie in den Bilddaten auf Basis eines Klassifikationsergebnisses identifiziert wird.
Die Schritte des Verfahrens können mit Hilfe eines Computerprogrammprodukts ausgeführt werden, welches direkt in einen Speicher eines programmierbaren Bildanalysesystems zur Verwendung in einem derartigen Bildbearbeitungssystem ladbar ist. Das Computerprogrammprodukt kann geeignete Programmcodemittel aufweisen, um die Schritte des Verfahrens zur automatischen Analyse der Bilddaten einer Struktur durchzuführen, wenn das Programmprodukt auf dem Bildanalysesystem ausgeführt wird.
Die abhängigen Ansprüche sowie die weitere Beschreibung enthalten jeweils besonders vorteilhafte Weiterbildungen und Ausgestaltungen der Erfindung, wobei das erfindungsgemäße Bildverarbeitungssystem auch entsprechend den abhängigen Verfahrensansprüchen weitergebildet sein kann.
Wie bereits erwähnt, kann das erfindungsgemäße Verfahren auf beliebige geeignete dreidimensionale Bilddaten von Strukturen angewendet werden, um ganz allgemein eine Anomalie in der Struktur zu finden. Im Folgenden wird, ohne die Erfindung in irgendeiner Weise dadurch zu beschränken, jedoch angenommen, dass die Bilddaten in einem CT-Bildgebungsschritt gewonnen wurden und dass die Struktur eine Knochenstruktur – insbeson dere eine dünne Knochenstruktur – ist. Ein Beispiel für einen dünnen Knochen kann ein Gesichtsknochen sein.
Da eine Segmentierung idealerweise jeden Voxel identifizieren sollte, der zu einem Knochen, vorzugsweise dünnen Knochen, gehört, umfasst der adaptive Segmentierungsschritt in dem erfindungsgemäßen Verfahren vorzugsweise eine Analyse des dreidimensionalen Voxel-Arrays zur Gewinnung von lokalen Strukturorientierungsinformationen und dann die Durchführung eines adaptives Schwellenwertverfahrens auf Basis der lokalen Strukturorientierungsinformationen. Unter dem Begriff ”adaptives Schwellenwertverfahren” ist zu verstehen, dass die Voxel nicht einfach mit einem Abschwächungs-Schwellenwert für Knochen verglichen werden. Stattdessen spielt die Charakteristik der lokalen Struktur in der Region des jeweiligen Voxels eine wichtige Rolle bei der Entscheidung, ob ein Voxel zu einer dünnen Knochenstruktur oder zu anderem Gewebe gehört. Durch die Verwendung der lokalen Orientierung bzw. lokalen Strukturinformationen in der beschriebenen Weise kann die Anzahl der Sampling-Artefakte signifikant reduziert werden. Durch Anwendung der oben beschriebenen Segmentierungstechnik im Rahmen des erfindungsgemäßen Verfahrens können Fehler wie Partialvolumenartefakte vermieden oder zumindest signifikant reduziert werden, da die Tensoren Informationen über die Planarität oder Linearität in der Nachbarschaft eines Voxels und deren Isotropie enthalten. Daher ist es viel unwahrscheinlicher, dass eine ”Nicht-Knochen-Region” als eine ”Knochen-Region” identifiziert wird oder umgekehrt. Durch Anwendung dieser Technik wird ein Voxel-Array mit einer genaueren Signalbasis effektiv resampelt. Unter Berücksichtigung der Charakteristiken, insbesondere von dünnen Knochenstrukturen, wird auf diesem Wege der Partialvolumeneffekt minimiert und Voxel, die zu Knochen gehören, werden besser identifiziert.
Nach Durchführung der Segmentierung liegt als Ergebnis ein Datensatz vor, welcher angibt, welche Teile des CT-Bilds höchstwahrscheinlich Knochen entsprechen und welche Teile höchstwahrscheinlich nicht.
Eine Merkmalsextraktion ist eine rechentechnisch aufwändige Prozedur. Daher wird in einem besonders bevorzugten Ausführungsbeispiel der Erfindung zunächst ein Begrenzungsbereich (auch ”Bonnding Box” oder ”interessierender Bereich” genannt) innerhalb des Voxel-Arrays identifiziert, und die Durchführung der Merkmalsextraktion erfolgt nur für bestimmte Voxel innerhalb des Begrenzungsbereichs, vorzugsweise für alle Voxel innerhalb des Begrenzungsbereichs und nur für diese. Auf diesem Wege wird unnötiger Rechenaufwand vorteilhafterweise reduziert.
In einem besonders bevorzugten Ausführungsbeispiel der Erfindung wird bei der Merkmalsextraktion eine Vielzahl von unterschiedlichen Merkmalen (für jeden Voxel innerhalb des Begrenzungsbereichs) extrahiert, welche eine Struktur vom Typus eines dünnen Knochens optimal charakterisieren. Dabei kann es jedoch sein, dass unterschiedliche Merkmale eine Berücksichtigung unterschiedlicher Anzahlen von Voxeln erfordern. Zum Beispiel kann eine Art von Merkmalen eines Voxels für den jeweiligen Voxel unter Berücksichtigung eines Nachbarvoxels extrahiert werden. Mit anderen Worten, um das Merkmal für diesen Voxel zu extrahieren, ist es nötig, auch Informationen eines benachbarten Voxels zu berücksichtigen. Ebenso kann es für eine andere Art von Merkmalen dieses Voxels nötig sein, mehrere Voxel innerhalb einer Nachbarschaft einer bestimmten Größe um diesen Voxel herum zu berücksichtigen. Daher wird in einem weiter bevorzugten Ausführungsbeispiel der Erfindung das Voxel-Array oder der Begrenzungsbereich für ein bestimmtes Merkmal in eine Vielzahl von dreidimensionalen Array-Blöcke unterteilt und die Merkmalsextraktion wird speziell auf den Voxeln dieser Array-Blöcke durchgeführt, wobei die Abmessungen der Array-Blöcke jeweils in Abhängigkeit von einer Merkmalsart und/oder einer Kontur der Struktur, z. B. der Knochenstruktur, in einem dem Array-Block entsprechenden Bereich des Voxel-Arrays gewählt werden. An den ”Kanten” einer Knochenstruktur können z. B. die Array-Blöcke relativ klein sein (z. B. ein 2 × 2 × 2 Array-Block), da jede Region außerhalb der Knochenstruktur nicht von Interesse ist, wogegen in ”zentralen” Regionen die Array-Blöcke größer sein können, wenn dies für den bestimmten Merkmalstyp erforderlich ist. Um die Merkmalsextraktion weiter zu verbessern, kann ein Array-Block auch einen benachbarten Array-Block überlappen, so dass so viele Informationen wie möglich von den Voxeln extrahiert werden. Der Betrag des Überlapps kann wiederum durch die Strukturkontur, z. B. Knochenkontur, in dieser Region und/oder durch die Art des zu extrahierenden Merkmals bestimmt sein. Durch diese Techniken können Variationen der Form hinsichtlich unterschiedlicher Subjekte (Alter, Geschlecht) besser kompensiert werden.
Wie bereits einleitend erwähnt, kann die Ausdehnung der Fehlausrichtung in einem Haarriss- oder Spannungsbruch so gering sein, dass diese Art von Bruch oftmals nicht in einer visuellen Überprüfung von Röntgen- oder CT-Bildern gesehen werden kann. Jedoch sind die Trabekel (Knochenbälkchen) eines Knochens in spezielle Richtungen ausgerichtet, um die auf den Knochen wirkenden Kräfte bestmöglich abzufangen, und ein Bruch des Knochens führt zu einer Störung in diesem trabekulären Muster in der Region. In einem bevorzugten Ausführungsbeispiel der Erfindung umfasst daher ein Merkmal, das innerhalb der Merkmalsextraktion extrahiert wird, ein Merkmal eines trabekulären Texturmusters, d. h. ein Merkmal, welches ein Muster des trabekulären Gewebes beschreibt. Vorzugsweise wird eine Mehrzahl von Merkmalen extrahiert, wobei jedes von den Merkmalen einige Charakteristika des Knochens hinsichtlich seiner trabekulären Natur erfasst bzw. beschreibt. Ein derartiges Merkmal eines trabekulären Texturmusters kann, wie noch detaillierter erläutert wird, beispielsweise ein fraktales Dimensionsmerkmal sein, ein Intensitätsmerkmal, ein Intensitätsgradientenmerkmal, ein Gabor-Orientierungskartenmerkmal, ein Merkmal eines Markov-Netzwerks (Markov Random Field) etc.
Das Ziel der Merkmalsextraktion ist es, eine Merkmalskarte für die Voxel innerhalb der interessierenden Knochenregion zu erhalten und diese Daten mit zuvor gewonnenen Daten für eine vergleichbare ”gesunde” Knochenregion zu vergleichen, um so jede Diskrepanz zu erkennen, welche anzeigt, ob eine zu untersuchende Knochenregion gebrochen ist oder nicht. Daher umfasst die Merkmalskarte in einem bevorzugten Ausführungsbeispiel der Erfindung einen Merkmalsvektor für jeden Voxel des dreidimensionalen Voxel-Arrays, besonders bevorzugt für jeden Voxel innerhalb des Begrenzungsbereichs. Ein Merkmalsvektor kann ein oder mehrere Merkmale für jeden Voxel umfassen, z. B. eine Sequenz der oben erwähnten Merkmale eines trabekulären Texturmusters. Vorzugsweise wird ein Merkmalsvektor in der gleichen Weise aufgebaut, d. h. unter Verwendung der gleichen Arten von Merkmalen, wie ein zugehöriger Merkmalsvektor in einer Merkmalskarte, welche aus gesunden Trainingsmustern erstellt wurde, so dass ein 1:1-Vergleich möglich ist.
Die zuvor gewonnenen Daten für eine vergleichbare ”gesunde” Knochenregion enthalten vorzugsweise also auch eine Merkmalskarte, in der die Merkmalsvektoren in der gleichen Art aufgebaut sind. Dabei kann eine ”mittlere” Merkmalskarte konstruiert werden, indem eine vergleichbare Knochenregion für eine Anzahl von Subjekten untersucht wird und die Ergebnisse für jedes Merkmal in jedem Merkmalsvektor gemittelt werden. Die Wahl der mittleren Merkmalskarte kann weitgehend von der Natur des zu untersuchenden Subjekts abhängen, z. B. kann eine Auswahl sowohl in Abhängigkeit des Knochentyps als auch des Alters und/oder Geschlechts getroffen werden. Naturgemäß steigt, wie dem Fachmann klar ist, der Informationsgehalt einer solchen mittleren Merkmalskarte mit der Anzahl der untersuchten Subjekte an. Ein Vergleich der Merkmalskarte mit einer zuvor gewonnenen mittleren Merkmalskarte berücksichtigt vorzugsweise jeden Voxel in der interessierenden Region, d. h. dass jeder Voxel mit seinem entsprechenden Voxel innerhalb der mittleren Merkmalskarte verglichen wird, um jeden relevanten Unterschied zu finden. Aus diesem Grund wird bei dem erfindungsgemäßen Verfahren die skalare Differenzkarte auf Basis der Merkmalskarte vorzugsweise so gewonnen, dass durch einen Vergleich eines Merkmalsvektors einer Merkmalskarte mit einem zuvor bestimmten zugehörigen mittleren Merkmalsvektor ein Eintrag für die skalare Differenzkarte erhalten wird.
Bei einem Vergleich des gebrochenen Knochens mit einem gesunden Knochen führt die Störung in der trabekulären Struktur zu einer Differenz zwischen den Merkmalen von entsprechenden Voxeln in den Bilddaten und denen der mittleren Merkmalskarte. Dabei ist es klar, dass auch ”gesunde” Voxel innerhalb der Bilddaten Merkmale aufweisen können, die leicht von den entsprechenden Merkmalen der mittleren Merkmalskarte abweichen. Daher wird bei der Gewinnung der skalaren Differenzkarte auf Basis der Merkmalskarte bei einem besonders bevorzugten Ausführungsbeispiel der Erfindung ein Eintrag für die skalare Differenzkarte verworfen, wenn dieser einen Wert unterhalb eines vorgegebenen Schwellenwerts aufweist, so dass die skalare Differenzkarte nur Einträge aufweist, die größer oder gleich dem vorgegebenen Schwellenwert sind. Im Grunde genommen kann ein Eintrag von Null oder nahe Null in der skalaren Differenzkarte so interpretiert werden, dass die durch den entsprechenden Voxel repräsentierte Region des Knochens gesund ist, d. h. nicht gebrochen, wogegen Einträge, die nicht Null sind bzw. einen Wert größer als oder gleich einem vorgegebenen Schwellenwert aufweisen, anzeigen können, dass die durch diesen Voxel erfasste Knochenregion tatsächlich gebrochen ist. Da die mittlere Merkmalskarte über alle gesunden Trainingsmuster berechnet wurde, sollte ein Bild einer gesunden Knochenstruktur eine Merkmalskarte ergeben, welche sehr ähnlich zur mittleren Merkmalskarte ist. Daher ist zu erwarten, dass die Differenzkarte dieses Bilds überwiegend kleine Werte aufweist. Auf der anderen Seite sollten in einem Bild eines gebrochenen Knochens durch den Bruch einige Störungen des trabekulären Musters verursacht sein. Daher sollte die zugehörige Merkmalskarte an einigen Positionen stark unterschiedlich zu der mittleren Merkmalskarte sein, und es ist daher zu erwarten, dass die zugehörige skalare Differenzkarte einige große Werte zeigt.
Auf Basis der Einträge in der skalaren Differenzkarte, die nicht Null sind, können die Voxel innerhalb des CT-Bilds gekennzeichnet werden, ob sie wahrscheinlich zu einem Pfad eines Haarrissbruchs gehören oder nicht. Hierzu wird bei dem erfindungsgemäßen Verfahren eine Klassifizierung durchgeführt und die Durchführung der Klassifizierung umfasst vorzugsweise die Anwendung zumindest eines Klassifikators, d. h. eines geeigneten Mustererkennungswerkzeugs, auf die Einträge der skalaren Differenzkarte, wobei ein Klassifikator einen Voxel, welcher mit einem Eintrag einer skalaren Differenzkarte verknüpft wird, in eine Klasse einer Gruppe von Klassen einklassifiziert, welche zumindest eine Anomalie-(oder ”Gebrochen”-)Klasse und eine Nicht-Anomalie-(oder ”Nicht-Gebrochen”-)Klasse umfasst. Ein Eintrag in der skalaren Differenzkarte ist mit einem Voxel und daher auch mit einem Merkmalsvektor in der Merkmalskarte verknüpft. Somit können Sätze oder Gruppen von Merkmalsvektoren auf Basis ihrer zugehörigen skalaren Differenzkarteneinträge ausgewählt und für jeden der Merkmalsvektoren jedes Satzes oder Gruppen von Voxeln können eine oder mehrere Regeln angewendet werden, um die Klassifizierung durchzuführen. Die Merkmale können unter Anwendung von Standard-Merkmalsselektionstechniken gewählt werden, die dem Fachmann bekannt sind, z. B. die Mahalanobis-Distanz, das Fisher-Kriterium etc. Mit solchen Techniken wird bestimmt, wie diskriminant ein Merkmalssatz in der Identifizierung von unterschiedlichen Klassen ist, in diesem Fall der zwei Klassen ”Anomalie” und ”Non-Anomalie”. Dabei kann jede geeignete Merkmalsselektionstechnik angewandt werden, um einen brauchbaren Merkmalssatz zu identifizieren.
Vorzugsweise werden bei der Klassifizierung im Rahmen der erfindungsgemäßen Methode eine oder mehrere Regeln einer Gruppe von Regeln, umfassend eine Maximum-Regel, eine Minimum-Regel, eine Produkt-Regel, eine Summen-Regel, eine Mehrheits-Regel (Majority Vote-Rule) angewandt werden, wie dies später noch detaillierter erläutert wird. Es wird darauf hingewiesen, dass diese Gruppe von Regeln nicht erschöpfend ist und auch andere geeignete Regeln angewandt werden können. Es sollte auch beachtet werden, dass die Gruppe von Klassen, die in der Klassifizierung genutzt werden, nicht notwendigerweise auf zwei Klassen beschränkt ist, sondern auch eine oder mehrere Klassen umfassen kann, sofern dies erforderlich ist.
Das Ergebnis der Klassifizierung ist, dass jeder Voxel, der wahrscheinlich mit einem Bruch verknüpft ist, als solcher gekennzeichnet oder markiert ist. Sinnvollerweise werden einzelne Voxel, die als gebrochen klassifiziert sind, ignoriert, aber eine Kette oder Gruppe von benachbarten Voxeln, welche alle als gebrochen klassifiziert sind, können verlässlich als Indikator interpretiert werden, dass dort tatsächlich in der Region des Knochens eine Fraktur vorliegt.
In einem besonders bevorzugten Ausführungsbeispiel des erfindungsgemäßen Verfahrens wird die Spur einer Strukturanomalie graphisch in einem Bild angezeigt, welches aus den dreidimensionalen Bilddaten gewonnen ist. Zum Beispiel kann als Diagnosehilfe eine Gruppe von Knochenvoxeln, welche als möglicherweise gebrochen klassifiziert wurde, graphisch mit einer Kontur oder Linie ergänzt werden, um die als gebrochen angenommene zugehörige Knochenregion hervorzuheben. Ein Radiologe oder Arzt kann dann das Bild begutachten, um zu einer Diagnose zu kommen und weitere erforderliche Behandlungen festzulegen.
Weitere Vorteile und Merkmale der vorliegenden Erfindung werden aus der nachfolgenden detaillierten Beschreibung in Verbindung mit den beigefügten Zeichnungen ersichtlich. Es ist klar, dass die Zeichnungen lediglich zum Zwecke der Illustration und nicht zur Definition der Grenzen der Erfindung dienen. Es zeigen:
1 ein Ablaufdiagramm der Schritte eines Beispiels eines erfindungsgemäßen Verfahrens,
2 eine schematische Darstellung einer bei Anwendung des erfindungsgemäßen Verfahrens gewonnenen Merkmalskarte, und
3 ein Blockdiagramm für ein Ausführungsbeispiel eines erfindungsgemäßen Bildverarbeitungssystems.
Wegen der hohen Anzahl von Gleichungen, die zur Beschreibung der Verfahrensschritte benötigt werden, sind die mathematischen Bezeichnungen, die in der nachfolgenden Figurenbeschreibung verwendet werden, nicht notwendigerweise durchgehend konsistent. Jedoch ist der Fachmann mit den verwendeten Gleichungen und Ausdrücken sowie der Art und Weise, diese zu interpretieren, vertraut.
Im Folgenden werden als Beispiel dreidimensionale CT-Bilddaten einer Knochenregion verwendet und das erfindungsgemäße Verfahren wird genutzt, um einen Bruch in der Knochenregion zu erkennen.
1 zeigt ein Ablaufdiagramm zur Durchführung einer Haarrissbruchanalyse nach einem Ausführungsbeispiel des erfindungsgemäßen Verfahrens. Zuerst wird in Schritt 10 ein CT-Datensatz gewonnen, welcher eine Anzahl von ”Schichten” (slices) des untersuchten Objekts umfasst. Dieser Datensatz kann als dreidimensionales Array von Voxeln gesehen werden, wobei jeder Voxel einem ”Punkt” in einem dreidimensionalen Bild des Objekts entspricht und einen für den entsprechenden Punkt rekonstruierten Abschwächungswert enthält.
Im Schritt 11 wird eine dreidimensionale Segmentierung auf dem gesamten 3D-Voxel-Array durchgeführt, um alle Voxel innerhalb des 3D-Bilds zu identifizieren, die höchstwahrscheinlich Knochen entsprechen. Dieser Segmentierungsschritt 11 umfasst eine Anzahl von Verfahrensschritten:
Im Schritt 110 wird für jeden Voxel unter Verwendung von Quadraturfiltern ein dreidimensionaler lokaler Orientierungstensor Γ berechnet, welcher Tensor Γ die Orientierung einer lokalen Nachbarschaft beschreibt und daher als Merkmalsdeskriptor für Kanten, Linien oder Ebene in dieser Region wirkt. Dabei wird eine Repräsentation eines dreidimensionalen lokalen Orientierungstensors durch eine Kombination der Ausgangsdaten von sechs polar unterschiedlichen Quadraturfiltern berechnet. Für jeden Filter wird ein lognormaler Quadraturfilter Q verwendet, welcher sphärisch separierbar ist und auf einer logarithmischen Skala eine Gaußsche radiale Frequenzfunktion aufweist: Q(u) = R(ρ)Dk(û) (1.1)
Hierbei ist u eine multidimensionale Frequenzvariable und R(ρ) und D_k(û) sind die radiale Funktion bzw. die Richtungsfunktion. Das Ergebnis dieser Filter wird dann zur Gewinnung des dreidimensionalen lokalen Orientierungstensors verwendet:
wobei q_k die Ausgangsgröße des k-ten Quadraturfilters ist und M_k eine duale Tensorbasis entsprechend einer Tensorbasis
mit x ^k als Filterrichtungsvektoren.
Es werden dann im Schritt 111 die drei Eigenwerte λ₁, λ₂, λ₃ für den Tensor Γ berechnet und in absteigender Reihenfolge sortiert, so dass λ₁ ≥ λ₂ ≥ λ₃ ≥ 0 gilt. Auf Basis der Verhältnisse zwischen den Eigenwerten des Orientierungstensors kann bestimmt werden, zu welcher ”Formkategorie” (shape category) die Nachbarschaft gehört. Unter Annahme, dass die Eigenwerte des Tensor Γ in der obigen Gleichung (1.2) die Bedingung λ1 ≥ λ2 ≥ λ3 ≥ 0, (1.4) erfüllen, treten drei Nachbarschaftsfälle auf:

(i) planar: λ₁ >> λ₂ ≈ λ₃
(ii) linear: λ₁ ≈ λ₂ >> λ₃
(iii) sphärisch: λ₁ ≈ λ₂ ≈ λ₃

Im Schritt 112 wird unter Verwendung des Verhältnisses von λ₂ zu λ₁ ein Planaritätsmaß O für den Voxel bestimmt. Für den planaren Fall wird die folgende Gewissheitsschätzung (certainty estimate) verwendet:
Für einen großen Wert von α sind alle Gewissheitsschätzwerte für den nicht-planaren Fall nahe Null.
Nachfolgend wird in Schritt 113 ein Schwellenwert T für eine räumliche Position x und für einen gewählten Wert O berechnet: T(x) = T0 – βO(x) (1.6)
Dabei ist T₀ ein globaler Schwellenwert, β ist eine Konstante und O(x) ist das Planaritätsmaß. Ein Planaritätsmaß wird gewählt, da dünne Knochenstrukturen schwierig zu segmentieren sind. Durch Anpassung des Schwellenwerts unter Nutzung eines Gewissheitsschätzwerts O für den planaren Fall können dünne Knochenregionen besser entdeckt werden. Da dünne Knochen eine solche planare Struktur haben, wird jeder Voxel, der ein durch den Gewissheitsschätzwert O definiertes hohes Planaritätsmaß aufweist, auf Basis des globalen Schwellenwerts heraussegmentiert, der durch die implizit durch den Wert O gegebenen lokalen Strukturinformationen modifiziert ist.
In einem letzten Schritt 114 wird ein Voxel als einen Knochen repräsentierend gekennzeichnet, wenn Γ > T ist, und im anderen Fall bleibt er unbeachtet. Insgesamt dient der Segmentierungsschritt 11 dazu, die Teile des 3D-Bilds oder Voxel V zu identifizieren, bei denen es sich höchstwahrscheinlich um Knochen handelt und die daher für die Zwecke der Bruchidentifizierung von Interesse sind. Zumindest eine relevante Teilmenge dieser segmentierten Bilddaten können durch einen Grenzbereich umfasst werden, und die nachfolgende Merkmalsextraktion muss nur noch auf den Voxeln innerhalb dieses Grenzbereichs durchgeführt werden.
Im nachfolgenden Schritt 12 wird dann in einer Schleife, welche einen Merkmalsextraktionsschritt 13 und einen Merkmalskartenhinzufügungsschritt 14 umfasst, für jede Merkmalsart ein adaptives Sampling durchgeführt
Das adaptive Sampling umfasst im Wesentlichen zunächst eine Unterteilung der Ergebnisse der Segmentierung in kleinere 3D-Blöcke von Voxeln (auch Array-Blöcke genannt). Benachbarte Array-Blöcke von Voxeln können dabei um einen bestimmten Betrag überlappen, wobei der Überlappbetrag, der beispielsweise irgendwo zwischen 30% und 70% liegen kann, während des Verfahrens oder vorab auf Basis von experimentiellen Trainingsdaten bestimmt werden kann. Die Größe jedes Blocks wird auf Basis der segmentierten Knochenregion bestimmt (z. B. kann eine Blockgröße basierend auf der Knochenkontur gewählt werden) und in Abhängigkeit von der im nachfolgenden Schritt zu extrahierenden Art des Merkmals. Der Schritt 12 des adaptiven Samplings wird für jedes Merkmal im nachfolgenden Merkmalsextraktionsschritt 13 durchgeführt, um so 3D-Blöcke mit einer Größe zu erhalten, die am besten für den jeweiligen Merkmalstyp geeignet sind.
Im Merkmalsextraktionsschritt 13 werden bestimmte Merkmale für jeden Voxel des Voxel-Arrays innerhalb des Begrenzungsbereichs extrahiert, um eine Merkmalskarte in Form einer Matrix aufzubauen, in der jede Zeile zu einem Voxel innerhalb des Begrenzungsbereichs (und somit zu einem Punkt innerhalb des dreidimensionalen Bildes) und jede Spalte zu einem Merkmal dieses Voxels gehört, so dass eine komplette Zeile einer fertigen Merkmalskarte letztlich einen Merkmalsvektor für diesen Voxel bildet. Bei diesem Ausführungsbeispiel des erfindungsgemäßen Verfahrens wird für jede Iteration des Merkmalsextraktionsschritts eine der folgenden Extraktionen durchgeführt:

(i) eine Extraktion ein Merkmal einer fraktalen Dimension;
(ii) eine Extraktion eines Lakunaritätsmerkmals;
(iii) eine Extraktion eines Gabor-Orientierungsmerkmals;
(iv) eine Extraktion eines Merkmals eines Markov-Netzwerks;
(v) eine Extraktion eines Intensitätsgradientenmerkmals.

13

130

14

Für jeden 2D-Schnitt des CT-Bilds wird durch abschätzen der fraktalen Dimension eine numerische Abschätzung der Rauhigkeit der trabekulären Textur gewonnen:
Für eine Oberfläche z = f(x, y) eines gegebenen Wertepaares (i, j) ist die Variation V_ε(i, j) als die Differenz zwischen den Extremwerten von f in einer ε-Umgebung von (x, y) definiert: Vε(i, j) = maxf(i, j) – minf(i, j) (2.1)wobei max und min über eine Scheibe oder zirkulare Region berechnet werden, die durch (x – i)2 + (y – j)2 ≤ e2 (2.2)gegeben ist. Die Variation Vf(ε) von f ist die Summe von V_ε(i, j) über die gesamte Oberfläche. Dieses Variationsverfahren schätzt die Minkowski-Bouligand fraktale Dimension durch drei minus der Steigung eines Linienfits [ln ε, ln Vf (ε)] nach der Methode der kleinsten Quadrate über einen geeigneten Bereich von ε ab.
Das so für jeden Voxel gewonnene Merkmal der fraktalen Dimension wird dann in Schritt 14 in die Merkmalskarte eingetragen und der Prozess kehrt zum Schritt 12 zurück, in welchem die Voxel-Arrays für das nächste Merkmal neu berechnet werden.
Im Schritt 131 werden Lakunaritätsmerkmale extrahiert, die ein quantitatives Maß für die Verteilung und die Größe der ”Löcher” in der betrachteten Knochenregion angeben. Lakunaritätsmerkmale für einen Voxel können für verschiedene Größenordnungen berechnet werden, so dass mehrere Lakunaritätsmerkmale für einen Voxel gewonnen werden können.
Um ein Lakunaritätsmerkmal zu gewinnen, wird dem Bild ein gleitendes Fenster oder Block mit einer Größe (r × r × r) überlagert und über das Bild als Abtastfenster bewegt, um Lakunaritätsdaten des Bilds an jedem Ort zu sammeln. Für jeden Abtastblock wird die ”Masse” M durch Aufaddierung der gesamten Anzahl von aktiven Stellen (entsprechend den ”Knochen-Voxeln” die einen Hinweis auf die Knochendichte geben) bestimmt. Es wird dann für jede Block-Masse M die Dichteverteilung [n(M, r)] berechnet, die für die spezielle Blockgröße (r) möglich ist. Im Weiteren wird dann die Verteilungsfunktion [Q(M, r)] der Wahrscheinlichkeit berechnet, indem [n(M, r)] durch die Gesamtanzahl der Blöcke der Größenordnung r geteilt wird.
Die Lakunarität Λ bei der Größenordnung r ist dann durch die mittlere quadratische Abweichung der Fluktuationen der Massenverteilungswahrscheinlichkeiten geteilt durch den zugehörigen quadratischen Mittelwert gegeben und wird ausgedrückt durch:
wobei die mittlere Wahrscheinlichkeitsverteilungsfunktion durch die Summe über M·Q(M, r) gegeben ist und die Varianzwahrscheinlichkeitsverteilungsfunktion durch die Summe über M²·Q(M, r) gegeben ist. Für jede Dimension wird die Merkmalskarte durch Hinzufügen des berechneten Lakunaritätsmerkmals zum Merkmalsvektor für diesen Voxel in Schritt 14 erweitert, bevor das Verfahren zum Schritt 12 zurückkehrt, in dem die Array-Blöcke wieder für das nächste Merkmal berechnet werden.
Im Schritt 132 wird dann das Texturmuster eines Voxels unter Verwendung von Gabor-Orientierungsmerkmalen gewonnen. Die Merkmalskarte kann dann in Schritt 14 wieder durch die Gabor-Orientierungsmerkmale für jeden Voxel erweitert werden.
Hierzu wird ein Satz von Gabor-Filtern auf jeden Block bzw. Array-Block angewandt, um die Orientierung des durch die trabekulären Linien gegebenen Texturmusters zu extrahieren. Da das Texturmuster von ”porösem” oder trabekulärem Knochen mit einem engen Band von Frequenzen verbunden ist, werden zwölf Orientierungsfilter verwendet, beginnend mit 0° unter Anstieg der Orientierung in Intervallen von 15°, so dass diese Filter genau den Frequenzinhalt der Knochenstruktur erfassen. Die dominierende Orientierung der Textur wird dann durch die Ori entierung des Gabor-Filters mit der größten Antwort identifiziert. Die Orientierungsvektoren außerhalb der Grenzkontur des Knochens werden auf den Null-Vektor 0 gesetzt. Eine Orientierungskarte ist durch M = [uijk] (2.4)gegeben, wobei u_ijk ein Einheitsvektor ist, der die Orientierung des Trabekulärgewebes an einem Gitterpunkt (i, j, k) repräsentiert. Die Merkmalskarte wird dann durch Hinzufügung des Gabor-Orientierungsmerkmals zum Merkmalsvektor für den jeweiligen Voxel im Schritt 14 erweitert, bevor ein Rückschritt zum Schritt 12 erfolgt, in dem die Array-Blöcke erneut für das nächste Merkmal berechnet werden.
Im Schritt 133 wird die Intensität eines Voxels durch eine Linearkombination der Intensitäten der benachbarten Voxel beschrieben. Die Merkmalskarte wird dann in Schritt 14 für jeden Voxel durch ein Texturmerkmal eines sog. Markov-Netzwerks erweitert.
Ein Markov-Netzwerk- oder Markov Random Field-(MRF)Texturmodell beschreibt die Intensität eines Pixels p als Linearkombination der Intensitäten der Nachbarpixel q:
wobei θ(p, q) die Modellparameter sind und ε(q) einen Null-Mittelwert und konstantes Varianzrauschen repräsentiert. Vorzugsweise wird vor der Extraktion der MRF-Texturmodellmerkmale eine Intensitätsnormalisierung durchgeführt.
Der Modellparameter θ(p, q) am Punkt p wird dann durch Minimierung des Fehlers E berechnet:
Die Modellparameter θ(p, q), mit p = (i, j, k), werden dann auf einen Einheitsvektor u_ijk normalisiert, um eine MRF-Texturkarte gemäß MMRF = [uijk] (2.7)zu bilden. Einträge innerhalb der Knochenkontur werden als Null-Vektoren angesehen. Die Merkmalskarte wird dann durch Hinzufügen der MRF-Intensitätsmerkmale zu den Merkmalsvektoren für die jeweiligen Voxel im Schritt 14 erweitert, bevor zum Verfahrensschritt 12 zurückgekehrt wird, in dem die Array-Blöcke erneut für das nächste Merkmal berechnet werden.
Schließlich werden für jeden Voxel im Schritt 134 die Intensitätsgradientenmerkmale berechnet, welche ein Maß für die Knochendichte in diesem Voxel ergeben. Die Merkmalskarte wird dann für jeden Voxel im Schritt 14 durch ein Intensitätsgradientenmerkmal erweitert.
Hierfür werden die Bilder so normiert, dass ihre mittleren Intensitäten und Kontraste ähnlich sind. Für eine am Punkt p zentrierte, gegebene Region R(p) wird eine Region für einen Punkt q gesucht, dessen Intensitätsdifferenz d_m am größten ist:
Die Intensitätsgradientenrichtung g(p) wird dann als Vektordifferenz
berechnet, wobei sgn() die Vorzeichenfunktion ist. Die Richtung von g wird dann so definiert, dass sie von einem Ort höherer Intensität zu einem Ort niedriger Intensität zeigt. Die Intensitätsgradientenrichtung wird dann an jedem Ort (i, j, k) innerhalb der Knochenkontur berechnet, um eine Intensitätsgradientenrichtungskarte zu erhalten: MIG = [uijk] (2.10)
Die Gradientenrichtung wird dabei für jeden Voxel, der innerhalb des festgelegten Begrenzungsbereichs aber außerhalb der Knochenkontur liegt, als Null definiert. Die Merkmalskarte wird dann für den jeweiligen Voxel durch Hinzufügen der Intensitätsgradientenmerkmale zum Merkmalsvektor im Schritt 14 erweitert.
Nachdem alle gewünschten Merkmalsarten der aus den Schritten 12, 13 und 14 bestehenden Schleife extrahiert wurden, kann die weitere Verarbeitung unter Verwendung der vervollständigten Merkmalskarte fortgesetzt werden.
Die Merkmalskarte kann wie in 2 als eine Matrix dargestellt werden, die so viele Zeilen n hat, wie die Gesamtanzahl von Voxel in dem segmentierten CT-Bild, und die Anzahl der Spalten m entspricht dabei der Anzahl der Merkmale, die für jeden Voxel extrahiert wurden. Ein Merkmalsvektor fv entspricht dabei einem Voxel innerhalb des CT-Bilds und wird durch m Merkmale beschrieben. In dem oben beschriebenen Ausführungsbeispiel entspricht das erste Merkmal f₂₁ des Merkmalsvektors fv beispielsweise der fraktalen Dimension des Voxels und die nachfolgenden Einträge beschreiben die Lakunarität, gemessen für unterschiedliche Größenordnungen, gefolgt von einer Anzahl von Gabor-Orientierungsmerkmalen (z. B. zwölf solcher Merkmale, wenn zwölf Orientierungs-Gabor-Filter benutzt werden). Das vorletzte Merkmal ist ein Markov Random Field-Merkmal und das letzte Merkmal f_2m entspricht dem Intensitätsgradienten dieses Voxels.
Nachdem die Merkmalskarte vervollständigt ist, kann dann im Schritt 15 die Analyse der CT-Daten mit dieser Karte fortgesetzt werden. Zunächst wird im Schritt 150 die Merkmalskarte mit einer mittleren Merkmalskarte verglichen, welche aus ver schiedenen gesunden Trainingsmustern M = [m_ijk] für jede Merkmalsart gebildet wird, die in einem vorherigen Trainingsverfahren ermittelt wurden:
Dabei ist m_ijk der mittlere Merkmalsvektor an der Position (i, j, k) innerhalb der Merkmalskarte M, n ist die Anzahl der gesunden Trainingsmuster, u_sijk ist der Einheitsmerkmalsvektor eines Trainingsbeispiels s am Ort (i, j, k) und c_ijk ist die Anzahl von Beispielen mit Merkmalsvektoren, die am Ort (i, j, k) nicht Null sind.
Für eine bestimmte Position (i, j, k) wird der zugehörige Wert in der mittleren Merkmalskarte als Null-Vektor gesetzt und als insignifikant betrachtet, wenn mehr als die Hälfte aller Trainingsbeispiel-Merkmalskarten an dieser Position einen Wert von Null haben. Diese Situation trifft üblicherweise nahe der Grenzkontur eines Knochens auf.
Die oben beschriebene Merkmalskarte ist eine Vektorkarte, deren Format nicht für eine Klassifikation günstig ist. Daher wird jede Merkmalskarte zunächst in eine skalare Differenzkarte konvertiert, indem jeder Eintrag in der Merkmalskarte mit dem entsprechenden Eintrag der mittleren Merkmalskarte verglichen wird, um die skalare Differenzkarte V = [v_ijk] zu erhalten. Jeder Eintrag v_ijk in der skalaren Differenzkarte ist ein Indikator für die Differenz zwischen einem Voxel an einer Position (i, j, k) in der Merkmalskarte, welche von den Bilddaten gewonnen wurde und einem entsprechenden Voxel für die mittlere Merkmalskarte und ist gegeben durch:
Die Einträge in dieser skalaren Differenzkarte liegen zwischen 0 und 1. Ein großer Wert v_ijk zeigt eine große Differenz an, wogegen ein Wert v_ijk, der nahe Null liegt, eine geringe Differenz anzeigt. Werte nahe bei Null bzw. Werte unterhalb einer vordefinierten Schwelle können folglich ignoriert werden.
Ein signifikanter Wert, der nicht Null ist, zeigt an, dass für diesen Voxel eine Differenz detektiert wurde und dass dieser Voxel deshalb einem Bruch in der Knochenregion zugeordnet werden könnte. Da die Merkmalskarte über alle gesunden Trainingsmuster generiert wurde, sollte ein gesunder Knochen zu einer Merkmalskarte führen, die relativ ähnlich der mittleren Merkmalskarte ist. Daher ist davon auszugehen, dass für ein solches Bild die Differenzkarte meistens kleine Werte aufweist. Auf der anderen Seite treten in einem Bild eines gebrochenen Knochens einige durch die Fraktur verursachte Störungen des trabekulären Musters auf. Daher wird die zugehörige Merkmalskarte sehr unterschiedlich zu der mittleren Merkmalskarte an einigen Positionen sein und es ist zu erwarten, dass die skalare Differenzkarte einige große Werte zeigt.
Diese skalare Differenzkarte wird dann als Basis für den nachfolgenden Klassifizierungsschritt 16 verwendet, indem Regionen des Knochens als ”gebrochen” oder ”gesund” klassifiziert werden. Hierzu werden in den Schritten 160, 161, 162, 163, 164 ein oder mehrere unterschiedliche Klassifikatoren auf die Merkmalsvektoren mit den skalaren Einträgen der Differenzkarte angewandt, welche nicht Null sind, um die a-posteriori-Wahrscheinlichkeit zu bestimmen, dass ein Voxel zu einem gesunden Knochenstück oder einem gebrochenen Knochenstück gehört.
Die Ergebnisse der Klassifikatoren werden dann im Schritt 165 kombiniert und jeder Voxel wird in geeigneter Weise gekennzeichnet oder entsprechend klassifiziert, dass er zur Klasse ”gebrochen” oder zur Klasse ”nicht-gebrochen” (gesund) gehört.
Jeder Klassifikator misst dabei die a-posteriori-Wahrscheinlichkeit P(ω_j|x_i), dass eine ”Probe” bzw. ”Muster” Z unter Verwendung eines gewählten Satzes von Merkmalsvektoren x_i zu einer bestimmten Klasse ω_j gehört. Mit anderen Worten, jeder Klassifikator repräsentiert das gegebene Muster durch ein Vektormaß, in diesem Fall die Wahrscheinlichkeit, dass ein Muster zur Klasse ”gesund” oder zur Klasse ”gebrochen” gehört. Die Probe bzw. das Muster Z enthält im Wesentlichen den Voxel-Satz, der zur Berechnung der Merkmalsvektoren verwendet wurde.
Bei dem hier beschriebenen Ausführungsbeispiel des Verfahrens werden im Klassifizierungsschritt folgende Regeln verwendet, wobei N die Anzahl der Klassifikatoren bezeichnet, P(ω_j) das zugehörige Vektormaß und P(ω_j|x_i) die a-posteriori-Wahrscheinlichkeit ist, dass das betrachtete Muster zu einer bestimmten Klasse gehört:

Maximum-Regel: Ordne Z → ω_k zu, wenn
Minimum-Regel: Ordne Z → ω_k zu, wenn
Produkt-Regel: Ordne Z → ω_k zu, wenn
Summen-Regel: Ordne Z → ω_k zu, wenn
Mehrheits-Regel: Ordne Z → ω_k zu, wenn

Ein Knochen wird als gebrochen klassifiziert, wenn eine vorbestimmte Anzahl der Klassifikatoren ein Muster als gebrochen klassifiziert. In der Maximum-Regel z. B. wird ein Muster (oder eine Gruppe von Voxeln) einer Klasse zugeordnet, die unter den zum Muster gehörigen Wahrscheinlichkeiten die höchste Wahrscheinlichkeit hat.
In einem abschließenden Anzeigeschritt 17 wird ein Pfad eines wahrscheinlichen Haarrissbruches graphisch mit den Voxeln eines oder mehrerer der CT-Schichtbilder dargestellt, so dass der Bruch für diagnostische Zwecke visuell beurteilt werden kann. Diese Ausgabe kann in einer zweidimensionalen Weise erfolgen, z. B. als eine farbige Linie, welche dem kalkulierten Pfad der Fraktur über einem 2D-Bild folgt, welches aus den 3D-Bilddaten rekonstruiert wurde, oder auch in einer dreidimensionalen Weise, wobei der Pfad der Fraktur mit größerer Genauigkeit betrachtet werden kann, z. B. auf einem Computerbildschirm.
3 zeigt ein Blockdiagramm für ein Bildverarbeitungssystem 1. In diesem Ausführungsbeispiel werden dreidimensionale Bilddaten D einer Knochenregion, die zuvor unter Verwendung z. B. eines Computertomographen erzeugt und in einem Speicher 313 hinterlegt wurden, an ein Bildanalysesystem 310 übergeben. Die Bilddaten D werden dann unter Verwendung der oben mit Hilfe der 1 und 2 beschriebenen Verfahrensschritte verarbeitet, um von Voxeln der Bilddaten D eine Merkmalskarte zu erhalten, welche dann bezüglich bekannter gesunder Daten untersucht wird, um mögliche Brüche im Knochen zu lokalisieren. Eine Mittelwertsmerkmalskarte wird aus Trainingsdaten TD gebildet, die wiederum aus einer geeigneten Datenquelle wie beispielsweise einem Speicher 311 mit einer Datenbank für bekannte gesunde Subjekte entnommen wurde. Wenn z. B. die dreidimensionalen Bilddaten von einem Schläfenbein entnommen wurden, kann die gemittelte Merkmalskarte unter Verwendung von Trainingsdaten TD berechnet werden, die von Bildern von Schläfenbeinen von gesunden Personen gemacht wurden. Jede Diskrepanz zwischen den Datensätzen D, TD wird für einen nachfolgenden Ausgabeschritt markiert oder gekennzeichnet. Es wird dann eine entsprechende Teilmenge von Voxeln der dreidimensionalen Bilddaten, welche die vermutete Fraktur umfassen, ebenso wie die Informationen, welche Voxels innerhalb der Teilmenge hervorgehoben werden müssen, an ein Bildausgabemittel 312 übermittelt, so dass die vermutete Fraktur von einem Diagnostiker, beispielsweise einem Arzt oder Radiologen, gesichtet werden kann, um die Informationen zu prüfen und eine entsprechende Diagnose anzufertigen.
Es wird abschließend noch einmal darauf hingewiesen, dass es sich bei dem vorhergehend detailliert beschriebenen Verfahren und der Systemarchitektur lediglich um bevorzugte Ausführungsbeispiele handelt, welche vom Fachmann in verschiedenster Weise modifiziert werden können, ohne den Bereich der Erfindung zu verlassen, soweit er durch die Ansprüche vorgegeben ist. Zum Beispiel ist die Erfindung nicht auf die Erkennung von Brüchen in Knochen beschränkt, sondern kann für anderes Material angewendet werden, z. B. poröses Material, welches zur Gewinnung von dreidimensionalen Bilddaten gescannt werden kann, und entsprechend der beschriebenen Technik analysiert werden kann um eine strukturelle Anomalie im Material zu detektieren.
Es wird der Vollständigkeit halber auch darauf hingewiesen, dass der Gebrauch der unbestimmten Artikel ”ein” oder ”eine” innerhalb dieser Anmeldung nicht ausschließt, dass die betreffenden Merkmale auch mehrfach vorhanden sein können und dass der Begriff ”umfassend” nicht weitere Schritte oder Elemente ausschließt.

Claims

Verfahren zur automatischen Analyse von Bilddaten (D) einer Struktur, welches Verfahren folgende Verfahrensschritte umfasst: – Bereitstellung von Bilddaten (D) in Form eines dreidimensionalen Voxel-Arrays, – Durchführung einer Segmentierung des Voxel-Arrays zur Ermittlung einer Voxel-Teilmenge, – Durchführung einer Merkmalsextraktion zumindest für bestimmte Voxel der Voxel-Teilmenge zur Erzeugung einer Merkmalskarte (FM), – Erzeugung einer skalaren Differenzkarte auf Basis der Merkmalskarte (FM), – Durchführung einer Klassifizierung anhand der skalaren Differenzkarte und – Identifizieren einer Strukturanomalie in den Bilddaten (D) auf Basis eines Klassifikationsergebnisses.
Verfahren nach Anspruch 1, wobei die Durchführung der Segmentierung folgende Schritte umfasst: – Analyse des dreidimensionalen Voxel-Arrays zur Gewinnung einer lokalen Strukturorientierungsinformation, – und Durchführung eines adaptiven Schwellenwertverfahrens auf Basis der lokalen Strukturorientierungsinformation.
Verfahren nach Anspruch 1 oder 2, wobei ein Begrenzungsbereich innerhalb des Voxel-Arrays festgelegt wird und die Durchführung der Merkmalsextraktion für bestimmte Voxel des Begrenzungsbereichs erfolgt.
Verfahren nach einem der vorstehenden Ansprüche, wobei das Voxel-Array oder der Begrenzungsbereich in eine Vielzahl von dreidimensionale Array-Blöcke unterteilt wird und wobei die Merkmalsextraktion für die Voxel in einem Array-Block durchgeführt werden, wobei die Abmessungen der Array-Blöcke in Abhängigkeit von einer Merkmalsart und/oder einer Kontur der Struktur in einem dem Array-Block entsprechenden Bereich des Voxel-Arrays gewählt werden.
Verfahren nach einem der vorstehenden Ansprüche, wobei ein bei der Merkmalsextraktion extrahiertes Merkmal ein Merkmal eines trabekulären Texturmusters umfasst
Verfahren nach Anspruch 5, wobei das Merkmal eines trabekulären Texturmusters eines der folgenden Merkmale umfasst: – eine fraktale Dimension, – ein Lakunaritätsmaß, – eine Gabor-Orientierung, – ein Markow-Netzwerk, oder – einen Intensitätsgradienten.
Verfahren nach einem der vorstehenden Ansprüche, wobei die Merkmalskarte (FM) für jeden Voxel der Voxel-Teilmenge des dreidimensionalen Voxel-Arrays einen Merkmalsvektor (fv) umfasst.
Verfahren nach Anspruch 7, wobei die Erzeugung einer skalaren Differenzkarte auf Basis der Bearbeitung der Merkmalskarte (FM) einen Vergleich eines Merkmalsvektors (fv) der Merkmalskarte (FM) mit einem vorab bestimmten entsprechenden gemittelten Merkmalsvektor zur Gewinnung eines Eintrags für die skalare Differenzkarte umfasst.
Verfahren nach Anspruch 8, wobei die Erzeugung einer skalaren Differenzkarte auf Basis der Bearbeitung der Merkmalskarte (FM) das Verwerfen eines Eintrags für die skalare Differenzkarte aufweist, wenn dieser einen Wert unterhalb eines vorgegebenen Schwellenwerts aufweist.
Verfahren nach Anspruch 8 oder 9, wobei die Durchführung der Klassifizierung zumindest die Anwendung eines Klassifikators auf Einträge der skalaren Differenzkarte umfasst, wobei ein Klassifikator ein Voxel, welches mit einem Eintrag der skalare Differenzkarte verknüpft ist, in eine Klasse einer Gruppe von Klassen einklassifiziert, die zumindest eine Anomalie-Klasse und eine Nicht-Anomalie-Klasse enthält.
Verfahren nach Anspruch 10, wobei innerhalb des Klassifizierungsschrittes eine oder mehrere Regeln einer Gruppe von Regeln, umfassend – eine Maximum-Regel, – eine Minimum-Regel, – eine Produkt-Regel, – eine Summen-Regel, – eine Mehrheits-Regel, angewandt werden.
Verfahren nach Anspruch 10 oder 11, wobei Voxel des Voxel-Arrays, die bei der Klassifizierung in die Anomalie-Klasse klassifiziert wurden, zur visuellen Darstellung einer Strukturanomalie in einem Bild genutzt werden.
Verfahren zur Ansteuerung einer Bildanzeigeeinrichtung zur Anzeige einer Strukturanomalie in einem Bild der Struktur, wobei das Bild aus dreidimensionalen Bilddaten (D) der Struktur gewonnen wird und die Strukturanomalie mittels eines Verfahren nach einem der Ansprüche 1 bis 12 identifiziert und im Bild graphisch dargestellt wird.
Bildverarbeitungssystem (1) zur automatischen Analyse von Bilddaten (D) einer Struktur mit einer Bilddatenquelle (313) zur Bereitstellung von Bilddaten (D) in Form eines dreidimensionalen Voxel-Arrays und mit einem Bildanalysesystem (310), welche zur Durchführung folgender Schritte ausgebildet ist: – Durchführung einer Segmentierung des Voxel-Arrays zur Ermittlung einer Voxel-Teilmenge, – Durchführung einer Merkmalsextraktion zumindest für bestimmte Voxel der Voxel-Teilmenge zur Erzeugung einer Merkmalskarte (FM), – Erzeugung einer skalaren Differenzkarte auf Basis der Merkmalskarte (FM), – Durchführung einer Klassifizierung anhand der skalaren Differenzkarte und – Identifizieren einer Strukturanomalie in den Bilddaten (D) auf Basis eines Klassifikationsergebnisses.
Computerprogrammprodukt, welches direkt in einen Speicher eines programmierbaren Bildanalysesystems (310) für ein Bildbearbeitungssystem (1) ladbar ist, mit Programmcode-Mitteln, um alle Schritte eines Verfahrens nach einem der Ansprüche 1 bis 13 auszuführen, wenn das Programmprodukt auf dem Bildanalysesystem (310) ausgeführt wird.