Bildverarbeitung für die Medizin 2006. Algorithmen - Systeme - Anwendungen Proceedings des Workshops vom 19. - 21. März 2006 in Hamburg

Informatik aktuell Herausgeber: W. Brauer im Auftrag der Gesellschaft für Informatik (GI)

H. Handels J. Ehrhardt A. Horsch H.-P. Meinzer T. Tolxdorff (Hrsg.)

Bildverarbeitung für die Medizin 2006 Algorithmen Systeme Anwendungen Proceedings des Workshops vom 19. – 21. März 2006 in Hamburg

1 23

Herausgeber Heinz Handels Jan Ehrhardt Universitätsklinikum Hamburg-Eppendorf Institut für Medizinische Informatik Martinistraße 52, 20246 Hamburg Alexander Horsch Technische Universität München, Klinikum rechts der Isar Institut für Medizinische Statistik und Epidemiologie Ismaninger Straße 22, 81675 München Hans-Peter Meinzer Deutsches Krebsforschungszentrum Abteilung für Medizinische und Biologische Informatik / H0100 Im Neuenheimer Feld 280, 69120 Heidelberg Thomas Tolxdorff Charité – Universitätsmedizin Berlin Institut für Medizinische Informatik, Biometrie und Epidemiologie Campus Benjamin Franklin Hindenburgdamm 30, 12000 Berlin

Bibliographische Information der Deutschen Bibliothek Die Deutsche Bibliothek verzeichnet diese Publikation in der Deutschen Nationalbibliografie; detaillierte bibliografische Daten sind im Internet über http://dnb.ddb.de abrufbar.

CR Subject Classification (2001): A.0, H.3, I.4, I.5, J.3, H.3.1, I.2.10, I.3.3, I.3.5, I.3.7, I.3.8, I.6.3 ISSN 1431-472-X ISBN-10 3-540-32136-5 Springer Berlin Heidelberg New York ISBN-13 978-3-540-32136-1 Springer Berlin Heidelberg New York Dieses Werk ist urheberrechtlich geschützt. Die dadurch begründeten Rechte, insbesondere die der Übersetzung, des Nachdrucks, des Vortrags, der Entnahme von Abbildungen und Tabellen, der Funksendung, der Mikroverfilmung oder der Vervielfältigung auf anderen Wegen und der Speicherung in Datenverarbeitungsanlagen, bleiben, auch bei nur auszugsweiser Verwertung, vorbehalten. Eine Vervielfältigung dieses Werkes oder von Teilen dieses Werkes ist auch im Einzelfall nur in den Grenzen der gesetzlichen Bestimmungen des Urheberrechtsgesetzes der Bundesrepublik Deutschland vom 9. September 1965 in der jeweils geltenden Fassung zulässig. Sie ist grundsätzlich vergütungspflichtig. Zuwiderhandlungen unterliegen den Strafbestimmungen des Urheberrechtsgesetzes. Springer Berlin Heidelberg New York Springer ist ein Unternehmen von Springer Science+Business Media springer.de © Springer-Verlag Berlin Heidelberg 2006 Printed in Germany Satz: Reproduktionsfertige Vorlage vom Autor/Herausgeber Gedruckt auf säurefreiem Papier SPIN: 11616719 33/3142-543210

Veranstalter IMI

Institut f¨ ur Medizinische Informatik Universit¨ atsklinikum Hamburg-Eppendorf, Hamburg GMDS Arbeitsgruppe Medizinische Bildverarbeitung der Gesellschaft f¨ ur Medizinische Informatik, Biometrie und Epidemiologie GI Fachgruppe Imaging und Visualisierungstechniken der Gesellschaft f¨ ur Informatik DGBMT Fachgruppe Medizinische Informatik der Deutschen Gesellschaft f¨ ur Biomedizinische Technik im VDE IEEE Joint Chapter Engineering in Medicine and Biology, German Section DAGM Deutsche Arbeitsgemeinschaft f¨ ur Mustererkennung BVMI Berufsverband Medizinischer Informatiker e.V. Lokaler Veranstalter Institut f¨ ur Medizinische Informatik Universit¨ atsklinikum Hamburg-Eppendorf, Hamburg Tagungsleitung und -vorsitz Prof. Dr. Heinz Handels Institut f¨ ur Medizinische Informatik Universit¨ atsklinikum Hamburg-Eppendorf Lokale Organisation Dipl.-Ing. Martin Dalladas Dr. Jan Ehrhardt Dipl.-Inform. Matthias F¨ arber Dipl.-Geogr. Dipl.-Systemwiss. Silke Hacker Priv.-Doz. Dr. Claus-J¨ urgen Peimann Dr. Andreas Pommert Renate Reche Dipl.-Ing. Martin Riemer Dipl.-Inform. Dennis S¨ aring Institut f¨ ur Medizinische Informatik Universit¨ atsklinikum Hamburg-Eppendorf

Organisation

Verteilte BVM-Organisation Prof. Dr. Heinz Handels, Dipl.-Ing. Martin Riemer Universit¨ at Hamburg und Universit¨ at zu L¨ ubeck (Begutachtung) Prof. II Dr. Alexander Horsch, Dipl.-Phys. Andreas Enterrottacher Technische Universit¨ at M¨ unchen (Tagungsband) Prof. Dr. Hans-Peter Meinzer und Dipl.-Ing. Matthias Baumhauer Deutsches Krebsforschungszentrum Heidelberg (Anmeldung) Prof. Dr. Thomas Tolxdorﬀ, Dagmar Stiller Charit´e – Universit¨ atsmedizin Berlin (Internetpr¨asenz)

Programmkomitee Prof. Dr. Til Aach, RWTH Aachen Prof. Dr. Dr. Johannes Bernarding, Universit¨ at Magdeburg Prof. Dr. Thorsten M. Buzug, Rhein Ahr Campus Remagen Prof. Dr. Hartmut Dickhaus, Fachhochschule Heilbronn Prof. Dr. Dr. Karl-Hans Englmeier, GSF Forschungszentrum Neuherberg Prof. Dr. Bernd Fischer, Universit¨ at zu L¨ ubeck Prof. Dr. Heinz Handels, Universit¨ at Hamburg Priv.-Doz. Dr. Peter Hastreiter, Universit¨ at Erlangen-N¨ urnberg Prof. Dr. Joachim Hornegger, Universit¨ at Erlangen-N¨ urnberg Prof. Dr. Ulrich Hoppe, Technische Universit¨ at Ilmenau Prof. II Dr. Alexander Horsch, TU M¨ unchen & Universit¨at Tromsø, Norwegen Priv.-Doz. Dr. Frithjof Kruggel, MPI f¨ ur neuropsychologische Forschung Leipzig Priv.-Doz. Dr. Thomas M. Lehmann, RWTH Aachen Prof. Dr. Dr. Hans-Gerd Lipinski, Fachhochschule Dortmund Prof. Dr. Tim L¨ uth, Technische Universit¨ at M¨ unchen Prof. Dr. Hans-Peter Meinzer, Deutsches Krebsforschungszentrum Heidelberg Prof. Dr. Heinrich M¨ uller, Universit¨ at Dortmund Prof. Dr. Heinrich Niemann, Universit¨ at Erlangen-N¨ urnberg Prof. Dr. Heinrich Martin Overhoﬀ, Fachhochschule Gelsenkirchen Prof. Dr. Dietrich Paulus, Universit¨ at Koblenz-Landau Prof. Dr. Heinz-Otto Peitgen, Universit¨ at Bremen Prof. Dr. Dr. Siegfried J. P¨ oppl, Universit¨ at zu L¨ ubeck Prof. Dr. Bernhard Preim, Universit¨ at Magdeburg Dr. Karl Rohr, Universit¨ at Heidelberg Prof. Dr. Georgios Sakas, Fraunhofer Institut Darmstadt Prof. Dr. Dietmar Saupe, Universit¨ at Konstanz Prof. Dr. Rainer Schubert, UMIT Innsbruck Prof. Dr. Hans-Siegfried Stiehl, Universit¨ at Hamburg Prof. Dr. Thomas Tolxdorﬀ, Charit´e – Universit¨atsmedizin Berlin Prof. Dr. Herbert Witte, Universit¨ at Jena Dr. Thomas Wittenberg, Fraunhofer Institut Erlangen

Organisation

VII

Preistr¨ ager des BVM-Workshops 2005 in Heidelberg Die BVM-Preise zeichnen besonders hervorragende Arbeiten aus. In 2005 wurden vom Programmkomitee folgende Preistr¨ ager ausgew¨ahlt: BVM-Preis 2005 f¨ ur die beste wissenschaftliche Arbeit 1. Preis: Eldad Haber und Jan Modersitzki Beyond Mutual Information: A Simple and Robust Alternative 2. Preis: Marcus Vetter, Martin Libicher, Ivo Wolf, Mehmet Ucar, Jochen Neuhaus, Mark Hastenteufel, G¨otz Martin Richter und Hans-Peter Meinzer Navigationssystem f¨ ur die perkutane CT-gesteuerte Radiofrequenz-Ablationstherapie von Lebertumoren 3. Preis: Sven Kabus, Astrid Franz und Bernd Fischer On Elastic Image Registration with Varying Material Parameters BVM-Preis 2005 f¨ ur den besten Vortrag 1. Preis: Eldad Haber und Jan Modersitzki Beyond Mutual Information: A Simple and Robust Alternative 2. Preis: Sven-Olaf Ropers, Thomas W¨ urﬂinger, Andr´e Alexander Bell, Alfred B¨ocking und Dietrich Meyer-Albrecht Automatische mikroskopische Relokalisation von Zellkern-Ensembles mit Hilfe eines multimodalen Szenenvergleichs 3. Preis: Sven Kabus, Astrid Franz und Bernd Fischer On Elastic Image Registration with Varying Material Parameters BVM-Preis 2005 f¨ ur die beste Poster- bzw. Softwarepr¨ asentation 1. Preis: Stephan Bischof, Tobias Weyand und Leif Kobbelt Snakes on Triangle Meshes 2. Preis: Silke Bommersheim, Ulf Tiede, Eike Burmester, Martin Riemer und Heinz Handels Simulation von Ultraschallbildern f¨ ur ein virtuelles Trainingssystem f¨ ur endoskopische Longitudinal-Ultraschalluntersuchungen 3. Preis: Caroline Langer, Markus Hardwiger und Katja B¨ uhler Interaktive diﬀusionsbasierte Volumensegmentierung auf Graﬁkhardware

Vorwort

In diesem Jahr ﬁndet der Workshop Bildverarbeitung f¨ ur die Medizin zum ersten Mal in der Freien und Hansestadt Hamburg statt. Er ist in dieser Form die neunte Veranstaltung. Die Bedeutung des Themas Bildverarbeitung f¨ ur die Medizin hat u ¨ ber die Jahre deutlich zugenommen. Die Bildverarbeitung ist eine Schl¨ usseltechnologie in verschiedenen medizinischen Bereichen wie der Diagnoseunterst¨ utzung, der OP-Planung und der bildgef¨ uhrten Chirurgie. Der BVM-Workshop konnte sich durch erfolgreiche Veranstaltungen in Freiburg, Aachen, Heidelberg, M¨ unchen, L¨ ubeck, Leipzig, Erlangen und Berlin als ein zentrales interdisziplin¨ ares Forum f¨ ur die Pr¨asentation und Diskussion von Methoden, Systemen und Anwendungen im Bereich der Medizinischen Bildverarbeitung etablieren. Ziel des Workshops ist die Darstellung aktueller Forschungsergebnisse und die Vertiefung der Gespr¨ ache zwischen Wissenschaftlern, Industrie und Anwendern. Der Workshop richtet sich ausdr¨ ucklich auch an Nachwuchswissenschaftler, die u ¨ ber ihre Diplom-, Promotions- und Habilitationsprojekte berichten wollen. Der Workshop wird vom Institut f¨ ur Medizinische Informatik des Universit¨ atskliniku*ms Hamburg-Eppendorf ausgerichtet. Die Organisation ist wie in den vergangenen Jahren auf Fachkollegen aus Hamburg, M¨ unchen, Heidelberg und Berlin verteilt, so dass die Organisatoren der vergangenen Jahre ihre Erfahrungen mit einﬂießen lassen k¨ onnen. Diese Aufgabenteilung bildet nicht nur eine starke Entlastung des lokalen Tagungsausrichters, sondern f¨ uhrt auch insgesamt zu einer Eﬃzienzsteigerung. Die etablierte webbasierte Einreichung und Begutachtung der Tagungsbeitr¨ age hat sich auch diesmal wieder bew¨ ahrt. Anhand anonymisierter Bewertungen durch jeweils drei Gutachter wurden aus 130 eingereichten Beitr¨agen 92 zur Pr¨ asentation ausgew¨ ahlt: 52 Vortr¨ age und 40 Poster. Die Qualit¨at der eingereichten Arbeiten war insgesamt sehr hoch. Die besten Arbeiten werden auch im Jahr 2006 mit BVM-Preisen ausgezeichnet. Am Tag vor dem wissenschaftlichen Programm werden drei Tutorien angeboten: Prof. Dr. Bernd Fischer und Priv.-Doz. Dr. Jan Modersitzki vom Institut f¨ ur Mathematik der Universit¨ at zu L¨ ubeck werden ein Tutorium zum Thema Medizinische Bildregistrierung halten. Dieses hochaktuelle Thema der Bildverarbeitung hat in den letzten Jahren eine st¨ urmische Entwicklung genommen. Ziel dieses Tutorials ist es, allgemeine Konzepte vorzustellen, die eine inhaltliche Einordnung von modernen medizinischen Bildregistrierungsverfahren erlauben. Insbesondere sollen die g¨ angigen Verfahren u ¨bersichtlich dargestellt werden. Das zweite Tutorium tr¨ agt den Titel Visualisierung f¨ ur die bildbasierte Diagnostik und Therapieplanung. Die Referenten der Universit¨at Magdeburg, Prof. Dr.-Ing. Bernhard Preim, Dipl.-Ing. Steﬀen Oeltze und Christian Tietjen werden grundlegende Methoden der medizinischen Visualisierung vorstellen und ihre Anwendungsm¨ oglichkeiten demonstrieren. Hierbei wird erl¨autert, welche Vor-

Vorwort

und Nachteile 2D- und 3D-Visualisierungen bieten und wie beide Darstellungsarten bei komplexen Problemstellungen in der Therapieplanung integriert werden k¨onnen. Der Schwerpunkt liegt auf der direkten Visualisierung medizinischer Volumendaten (Volume Rendering). In einem dritten Tutorium wird Herr Dr. Ph. Lahorte vom Europ¨aischen Patentamt in M¨ unchen u ¨ber das Patentrecht und geistiges Eigentum in der Forschung im Bereich der Medizinischen Bildverarbeitung referieren. Ziel dieses Tutoriums ist es, eine Einf¨ uhrung in das Gebiet des geistigen Eigentums allgemein und des Patentwesens im Besonderen zu geben und deren Bedeutung f¨ ur die Forschung in der Medizinischen Bildverarbeitung aufzuzeigen. Es wird ein ¨ Uberblick u ¨ber verschiedene Formen geistigen Eigentums mit dem Schwerpunkt Patente gegeben. Die gesetzliche Bedeutung eines Patents wird ebenso diskutiert wie verschiedene Verfahren, ein Patent zu erlangen. Anhand der Bewertungen der Gutachter wurden 95 Beitr¨age f¨ ur den Workshop in Vortrags- und Postersession sowie Softwaredemonstrationen ausgew¨ahlt und in diesem Tagungsband abgedruckt. Die Internetseiten des Workshops bieten ausf¨ uhrliche Informationen u ¨ber das Programm und organisatorische Details rund um den Workshop. Sie sind abrufbar unter der Adresse: http://www.bvm-workshop.org Wie schon in den letzten Jahren, wurde der Tagungsband auch diesmal als LATEXProjekt erstellt und in dieser Form an den Verlag u ¨bergeben. Von den 92 Beitr¨ agen wurden erfreulicherweise 82 von den Autoren bereits im LATEX-Format eingereicht (das sind 89% und damit 11% mehr als im Vorjahr). Die 10 im Winword-Format abgefassten Arbeiten wurden konvertiert und nachbearbeitet. Die Vergabe von Schlagworten nahmen die Autoren selbst vor. Die Literaturverzeichnisse s¨ amtlicher Beitr¨ age wurden wieder mit BibTEX generiert. Der gesamte Erstellungsprozess erfolgte ausschließlich u ¨ber das Internet. Die Herausgeber dieser Proceedings m¨ ochten allen herzlich danken, die zum Gelingen des BVM-Workshops 2006 beigetragen haben: Den Autoren f¨ ur die rechtzeitige und formgerechte Einreichung ihrer qualitativ hochwertigen Arbeiten, dem Programmkomitee f¨ ur die gr¨ undliche Begutachtung, den Referenten der Tutorien sowie den Mitarbeitern des Instituts f¨ ur Medizinische Informatik der Universit¨ at Hamburg f¨ ur ihre tatkr¨ aftige Unterst¨ utzung bei der Organisation und Durchf¨ uhrung des Workshops. Frau Dagmar Stiller vom Institut f¨ ur Medizinische Informatik, Biometrie und Epidemiologie der Charit´e, Universit¨ atsmedizin Berlin, danken wir f¨ ur die engagierte Mithilfe bei der Erstellung und Pﬂege der Internetpr¨ asentation. Herrn Andreas Enterrottacher vom Institut f¨ ur Medizinische Statistik und Epidemiologie der Technischen Universit¨ at M¨ unchen danken wir f¨ ur die tatkr¨ aftige Mitarbeit bei der Erstellung der Workshop-Proceedings. Dem Springer-Verlag, der nun schon den neunten Tagungsband zu den BVM-Workshops herausbringt, wollen wir f¨ ur die gute Kooperation ebenfalls unseren Dank aussprechen. F¨ ur die webbasierte Durchf¨ uhrung des Reviewingprozesses geb¨ uhrt Herrn Dipl.-Ing. Martin Riemer vom Institut f¨ ur Medizinische Informatik der Universit¨ at Hamburg unser Dank.

Vorwort

F¨ ur die ﬁnanzielle Unterst¨ utzung bedanken wir uns bei den Fachgesellschaften und der Industrie. Wir w¨ unschen allen Teilnehmerinnen und Teilnehmern des Workshops BVM 2006 lehrreiche Tutorials, viele interessante Vortr¨age, Gespr¨ache an den Postern und der Industrieausstellung sowie interessante neue Kontakte zu Kolleginnen und Kollegen aus dem Bereich der Medizinischen Bildverarbeitung.

Januar 2006 Heinz Handels, Jan Ehrhardt (Hamburg) Alexander Horsch (M¨ unchen) Hans-Peter Meinzer (Heidelberg) Thomas Tolxdorﬀ (Berlin)

Inhaltsverzeichnis

Die fortlaufende Nummer am linken Seitenrand entspricht den Beitragsnummern, wie sie im endg¨ ultigen Programm des Workshops zu ﬁnden sind. Dabei steht V f¨ ur Vortrag, P f¨ ur Poster.

Bildanalyse V01 Thies C, Schmidt-Borreda M, Lehmann T: Einsatz von Klassiﬁkatoren zum Lernen von Objektbeschreibungen aus hierarchisch partitionierten Bildern . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

V02 Werner R, Ehrhardt J, Frenzel T, S¨aring D, Low D, Handels H: Rekonstruktion von 4D-CT-Daten aus r¨ aumlich-zeitlichen CTSegmentfolgen zur Analyse atmungsbedingter Organbewegungen . . .

V03 J¨ager P, Vogel S, Knepper A, Kraus T, Aach T: 3D-Erkennung, Analyse und Visualisierung pleuraler Verdickungen in CT-Daten . . .

V04 Hensel M, Pralow T, Grigat R-R: LAST Filter for Artifact-Free Noise Reduction of Fluoroscopic Sequences in Real-Time . . . . . . . . . .

V05 Hahn T, Condurache AP, Aach T, Scharfschwerdt M, Misfeld M: Automatic In-Vitro Oriﬁce Area Determination and Fluttering Analysis for Tricuspid Heart Valves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

V06 Lamecker H, Wenckebach TH, Hege H-C, Duda GN, Heller MO: Atlas-basierte 3D-Rekonstruktion des Beckens aus 2D-Projektionsbildern . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

V07 Menze BH, Kelm BM, Heck D, Lichy MP, Hamprecht FA: Machine-Based Rejection of Low-Quality Spectra and Estimation of Brain Tumor Probabilities from Magnetic Resonance Spectroscopic Images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

V08 Oehler M, Buzug TM: Maximum-Likelihood-Ansatz zur Metallartefaktreduktion bei der Computertomographie . . . . . . . . . . . .

V09 Kier C, Meyer-Wiethe K, Seidel G, Aach T: UltraschallPerfusionsbildgebung f¨ ur die Schlaganfalldiagnostik auf Basis eines Modells f¨ ur die Destruktionskinetik von Kontrastmittel . . . . . . . . . . . .

XIV

Inhaltsverzeichnis

P01 Hensel M, Lundt B, Pralow T, Grigat R-R: Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

P02 Kelm BM, Menze BH, Neﬀ T, Zechmann CM, Hamprecht FA: CLARET: A Tool for Fully Automated Evaluation of MRSI with Pattern Recognition Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

P03 S¨aring D, Ehrhardt J, Stork A, Bansmann MP, Lund GK, Handels H: Analysis of the Left Ventricle After Myocardial Infarction Combining 4D Cine-MR and 3D DE-MR Image Sequences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

P04 Lessmann B, Nattkemper TW, Huth J, Loyek C, Kessar P, Khazen M, Pointon L, Leach MO, Degenhard A: Content Based Image Retrieval for Dynamic Time Series Data . . . . . . . . . . . . . . . . . . . .

P05 Varini C, Lessmann B, Degenhard A, Hans V, Nattkemper T: Visual Exploration of Pathology Images by a Discrete Wavelet Transform Preprocessed Locally Linear Embedding . . . . . . . . . . . . . . . .

P06 Fischer B, Winkler B, Thies C, G¨ uld MO, Lehmann TM: Strukturprototypen zur Modellierung medizinischer Bildinhalte . . . .

P07 Wirjadi O, Breuel TM, Feiden W, Kim Y-J: Automated Feature Selection for the Classiﬁcation of Meningioma Cell Nuclei . . . . . . . . . .

P08 Iserhardt-Bauer S, Schoell S, Hammen T, Stefan H, Doerﬂer A, Hastreiter P: Eﬃcient Atlas-Based Analysis of the Hippocampus . . .

Segmentierung V10 Brox T, Kim Y-J, Weickert J, Feiden W: Fully-Automated Analysis of Muscle Fiber Images with Combined Region and Edge-Based Active Contours . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

V11 Gr¨ unerbl A, Fritscher K, Blauth M, Kuhn V, Schubert R: Shape-Based 3D Level Set Segmentation of the Proximal Femur in CT-Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

V12 Dornheim L, Dornheim J, Seim H, T¨onnies K: Aktive Sensoren: Kontextbasierte Filterung von Merkmalen zur modellbasierten Segmentierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

V13 Stehle TH, Ecabert O: Vergleich von Geschwindigkeitsfunktionen zur Segmentierung der Koronararterien aus CT-Volumina mit Front-Propagation-Algorithmen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

Inhaltsverzeichnis

V14 Seim H, Dornheim J, Preim U: Ein 2-Fronten-Feder-Masse-Modell zur Segmentierung von Lymphknoten in CT-Daten des Halses . . . . . 106 V15 Schreiber J, Schubert R, Kuhn V: Femur Detection in Radiographs Using Template-Based Registration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 V16 von Berg J, Lorenz C: A Statistical Geometric Model of the Heart . 116 V17 Jetzek F, Rahn C-D, Dreschler-Fischer L: Ein geometrisches Modell f¨ ur die Zellsegmentierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 V18 Br¨ uß C, Strickert M, Seiﬀert U: Towards Automatic Segmentation of Serial High-Resolution Images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 V19 Chen L, Wagenknecht G: Topology Correction for Brain Atlas Segmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 V20 Rink K, T¨orsel A-M, T¨onnies K: Segmentation of the Vascular Tree in CT Data Using Implicit Active Contours . . . . . . . . . . . . . . . . . . 136 V21 Bruijns J, Peters FJ, Berretty RPM, Barenbrug B: Shifting of the Aneurysm Necks for Enhanced Aneurysm Labelling . . . . . . . . . . . . . . . . 141 P09 Sch¨onmeyer R, Rotarska-Jagiela A, Prvulovic D, Haenschel C, Linden DEJ: Vollautomatische Segmentierung der weißen Hirnsubstanz oberhalb der Seitenventrikel aus kernspintomographischen Datens¨ atzen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 P10 Popovic A, Engelhardt M, Radermacher K: Knowledge-Based Segmentation of Calvarial Tumors in Computed Tomography Images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 P11 Zerfaß T, Mues-Hinterw¨aller S, Paulus D, Wittenberg T: Live-Wire Segmentierung f¨ ur hochaufgel¨ oste Farbbilder mit optimierter Graphensuche . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 P12 Kov´ acs T, Cattin P, Alkadhi H, Wildermuth S, Sz´ekely G: Automatic Segmentation of the Vessel Lumen from 3D CTA Images of Aortic Dissection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 P13 Heimann T, Wolf I, Meinzer H-P: Automatische Erstellung von gleichm¨ aßig verteilten Landmarken f¨ ur statistische Formmodelle . . . 166 P14 B¨ohler T, Boskamp T, M¨ uller H, Hennemuth A, Peitgen H-O: Evaluation of Active Appearance Models for Cardiac MRI . . . . . . . . . 171

XVI

Inhaltsverzeichnis

P15 G¨ob S, Maier T, Benz M, Lowitzsch S, Wittenberg T, Kullmann W, Nkenke E, Neukam FW, H¨ausler G: Automatische Segmentierung der Gewebegrenzen in 2D-Ultraschalldaten aus der Mund-, Kieferund Gesichtschirurgie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176 P16 Riegg T, Zucker U, Horsch A: Intelligente Kantendetektion in endoskopischen Ultraschallbildern mit dem Centered-CompassFilter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181 P17 W¨orz S, Rohr K: Limits on Estimating the Width of Thin Vessels in 3D Medical Images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186 P18 Wilharm S, Maier T, Benz M, Sußner G, Lowitzsch S, H¨ausler G, Greiner G: Weichgewebemodellierung durch Fl¨acheninterpolation heterogen verteilter Messdaten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191

Registrierung V22 Winter S, Brendel B, Pechlivanis I, Schmieder K: Registrierung verschiedener Knochenstrukturen in Ultraschall- und CT-Daten anhand von pr¨ a- und intraoperativen Patientendatens¨atzen . . . . . . . . 196 V23 Franz A, Carlsen IC, Renisch S: An Adaptive Irregular Grid Approach Using SIFT Features for Elastic Medical Image Registration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201 V24 W¨orz S, Rohr K: New Approximating Gaussian Elastic Body Splines for Landmark-Based Registration of Medical Images . . . . . . . 206 ˇ V25 Cech P, Andronache A, Wang L, Sz´ekely G, Cattin P: Piecewise Rigid Multimodal Spine Registration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 211 V26 Mattes J, Steingruber I, Netzer M, Fritscher K, Kopf H, Jaschke W, Schubert R: A Robust Semi-automatic Procedure for Motion Quantiﬁcation of Aortic Stent Grafts Using Point Set Registration . 216 P19 Hahn D, Wolz G, Sun Y, Sauer F, Hornegger J, Kuwert T, Xu C: Utilizing Salient Region Features for 3D Multi-modality Medical Image Registration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221 P20 Feldmann T, Bouattour S, Paulus D, Deinzer F: Kombination ¨ verschiedener Ahnlichkeitsmaße f¨ ur die 2D/3D-Registrierung von R¨ ontgenbildern mittels Demokratischer Integration . . . . . . . . . . . . . . . . 226 P21 Braumann U-D, Einenkel J, Horn L-C, Kuska J-P, L¨oﬄer M, Scherf N, Wentzensen N: Registration of Histologic Colour Images of Diﬀerent Staining . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 231

Inhaltsverzeichnis

XVII

P22 J¨ager F, Han J, Hornegger J, Kuwert T: Wissensbasierte nicht-starre Registrierung von SPECT/CT Datens¨atzen . . . . . . . . . . . 236 P23 Erbacher M, Korosoglou G, Dickhaus H: Computergest¨ utzte Auswertung koronarangiographischer Bildfolgen hinsichtlich des Myocardialen Blushgrades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 241 P24 Han J, Berkels B, Rumpf M, Hornegger J, Droske M, Fried M, Scorzin J, Schaller C: A Variational Framework for Joint Image Registration, Denoising and Edge Detection . . . . . . . . . . . . . . . . . . . . . . . 246 P25 Palm C, Vieten A, Bauer D, Pietrzyk U: Evaluierung von Registrierungsstrategien zur multimodalen 3D-Rekonstruktion von Rattenhirnschnitten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 251 P26 Ehrhardt J, S¨aring D, Handels H: Interpolation of Temporal Image Sequences by Optical Flow Based Registration . . . . . . . . . . . . . . . . . . . . . 256 P27 K¨ohn A, Drexl J, Ritter F, K¨onig M, Peitgen HO: GPU Accelerated Image Registration in Two and Three Dimensions . . . . . 261

Visualisierung V27 Baer A, Tietjen C, Spindler M, Preim B: Hardwaregest¨ utztes Stippling von medizinischen Oberﬂ¨ achenmodellen . . . . . . . . . . . . . . . . . . 266 V28 Merhof D, Sonntag M, Enders F, Nimsky C, Hastreiter P, Greiner G: Streamline Visualization of Diﬀusion Tensor Data Based on Triangle Strips . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 271 V29 Hacker S, Handels H: Repr¨ asentation und Visualisierung von 3D-Formvarianten von Organen f¨ ur die medizinische Ausbildung . . . 276 V30 Bischoﬀ S, Kobbelt L: Extracting Consistent and Manifold Interfaces from Multi-valued Volume Data Sets . . . . . . . . . . . . . . . . . . . . 281 V31 Mang A, Wagner M, M¨ uller J, Fuchs M, Buzug TM: Restoration of the Sphere-Cortex Homeomorphism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 286 P28 Oeltze S, Kuß A, Hennemuth A, K¨ uhnel C, Preim B: Integrierte Visualisierung von Anatomie und Perfusion des Myokards zur Fr¨ uherkennung der Koronaren Herzkrankheit . . . . . . . . . . . . . . . . . . . . . . 291 P29 M¨ uhler K, Bade R, Preim B: Skriptbasierte Animationen f¨ ur die Operationsplanung und Ausbildung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 296

XVIII Inhaltsverzeichnis

P30 Wegner I, Wolber P, Gebhard B, Vetter M, Meinzer H-P: Verzerrung einer virtuellen Szene zur Erweiterung der Realit¨at eines Endoskopiebildes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 301 P31 Seitel M, Vaidya V, Raghu K, Mullick R: Modeling Dynamic Aspects in Virtual Interactive Environments . . . . . . . . . . . . . . . . . . . . . . . 306 P32 M¨onch T, Bernarding J: Simulation und kombinierte Visualisierung von aktivierten Hirnarealen, Diﬀusionstensordaten und Magnetenzephalographie (MEG)-Signalverl¨aufen . . . . . . . . . . . . . . . . . . . 311

Navigation und Tracking V32 Richter D, Egger J, Straßmann G: Ein Algorithmus zur Positionsbestimmung von Patienten und Biopsienadeln mit einem Vierkamerasystem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 316 V33 Riefenstahl N, Walke M, Michaelis B, Gademann G: Multimodale Bilddatenfusion zur verbesserten Patientenlagerung und -¨ uberwachung in der Strahlentherapie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 321 V34 Baumhauer M, Richter G, Gutt C, Rassweiler J, Meinzer H-P, Vetter M: Entwicklung eines Navigationssystems f¨ ur die laparoskopische Prostatektomie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 326 V35 Kenngott H, Neuhaus J, Gutt C, Wolf I, Meinzer H-P, Vetter M: Entwicklung eines Navigationssystems f¨ ur die telemanipulatorgest¨ utzte Oesophagektomie . . . . . . . . . . . . . . . . . . . . . . . . 331 P33 Meyer B, Tietjen C, Preim B: Schichtbasierte Illustration medizinischer Volumendaten zur intraoperativen Navigation . . . . . . . 335 P34 Kahrs LA, Raczkowsky J, W¨orn H: Optische Vermessung mittels kodierten Lichts von variabel reﬂektierenden Oberﬂ¨achen zur Registrierung oder Dokumentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 340 P35 Krause UHW, Nagel M, Seibel RMM: IGS (Image Guided Surgery) – Phantomversuche und erste klinische Erfahrungen mit einem neuartigen CT-basierten Navigationssystem . . . . . . . . . . . . . . . . . . . . . . . 345 P36 Richter D, Bekkaoui F, Mostarkic Z, Straßmann G: Tetraoptisches Kamerasystem zur rahmenlosen Repostionierung und respirativen ¨ Uberwachung in der extrakraniellen Hochpr¨azisionsbestrahlung . . . . 350 P37 Gr¨oger M, Arbter K, Hirzinger G: Analysis of Colour Distributions of Anodised Titanium Clips and the Heart Surface for Tracking . . . . 355

Inhaltsverzeichnis

XIX

Visible Light V36 Tscherepanow M, Z¨ollner F, Kummert F: Segmentierung ungef¨ arbter, lebender Zellen in Hellfeld-Mikroskopbildern . . . . . . . . . . 359 V37 Yang S, K¨ohler D, Teller K, Cremer T, Eils R, Rohr K: Non-rigid Registration of 3D Microscopy Images for the Normalization of Diﬀerent Cell Nuclei . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 364 V38 Fehr J, Sauer C, Kurz H, Ronneberger O, Burkhardt H: Identiﬁkation von Zellen in intaktem Gewebe . . . . . . . . . . . . . . . . . . . . . . 369 V39 Harder N, Neumann B, Held M, Liebel U, Erﬂe H, Ellenberg J, Eils R, Rohr K: Automated Analysis of Mitotic Phenotypes in Fluorescence Microscopy Images of Human Cells . . . . . . . . . . . . . . . . . . . 374 V40 Braumann U-D, Franke H, Hengstler J, Kuska J-P, Weber M: Graph-Based Quantiﬁcation of Astrocytes . . . . . . . . . . . . . . . . . . . . . . . . . 379 P38 Gladilin E, Eils R, Rohr K: 3D Formnormalisierung von Zellkernen mit Hilfe einer elastischen Kugelabbildung . . . . . . . . . . . . . . . . . . . . . . . . . 384

Simulation und Planung V41 Gladilin E, Ivanov A, Roginsky V: Biomechanische Optimierung individueller craniofazialer Implantate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 389 V42 Tropp A, Giesel FL, Dickhaus H, Bendl R, Neﬀ T: Integration der T2*-Perfusionsmessung niedergradiger Gliome in die Strahlentherapieplanung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 394 V43 von Jan U, Sandk¨ uhler D, Kirsch L, Maas S, R¨ uhmann O, Overhoﬀ HM: Ultrasound Volume Guided Navigated Implantation of the Humeral Part of a Shoulder Prosthesis . . . . . . . . . . . . . . . . . . . . . . 399 V44 Valvoda JT, Ullrich S, Kuhlen T, Bischof CH: Interactive Biomechanical Modeling and Simulation of Realistic Human Musculature in Virtual Environments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 404 P39 Bade R, Riedel I, Schmidt L, Oldhafer KJ, Preim B: Combining Training and Computer-Assisted Planning of Oncologic Liver Surgery . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 409

Inhaltsverzeichnis

Endoskopie und minimalinvasive Chirurgie V45 Rilk ME, Winkelbach S, Wahl FM: Partikelﬁlter-basiertes Tracking chirurgischer Instrumente in Endoskopbildern . . . . . . . . . . . . . . . . . . . . . . 414 V46 Wengert C, Reeﬀ M, Cattin PC, Sz´ekely G: Fully Automatic Endoscope Calibration for Intraoperative Use . . . . . . . . . . . . . . . . . . . . . . 419 V47 Rupp S, Winter C, Wittenberg T: Camera Calibration from Fiberscopic Views with Accuracy Evaluation . . . . . . . . . . . . . . . . . . . . . . 424 V48 Winter C, Weisensel S, Rupp S, Wittenberg T: Auﬂ¨ osungssteigerung von ﬁberskopischen Bildsequenzen im Ortsraum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 429

Freie Themen V49 Baron S, Orhan G, Hornung O, Blume H, Misske J, Norozi K, Wessel A, Yelbuz TM, Heimann B: Konstruktion und Etablierung einer Klimakammer f¨ ur die Untersuchung der embryonalen Herzentwicklung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 434 V50 Weber M, Ruiter NV, Schwarzenberg G, Zapf M, M¨ uller TO: Ultraschallsimulation f¨ ur die Ultraschall-Computertomographie . . . . 439 V51 Krings M, Mahnken A, Schroer C, Patommel J, Kalender W, Glasmacher B: CT, µ-CT und µ-Tomographie (Synchrotron) der in vitro Kalziﬁzierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444 V52 Weichert F, Linder R, Landes CA, Groh A, W¨ underlich R, Wagner M: Photorealistische Generierung histologischer und histopathologischer Schnittpr¨ aparate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 449 P40 Rink E, D¨ urschmied D, Harder D, Zhou Q, Freund G, Rossknecht A, Hehrlein C: Contrast Enhanced Ultrasound Perfusion Imaging . . . . . 454

Kategorisierung der Beitr¨ age . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 459 Autorenverzeichnis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 461 Stichwortverzeichnis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 465

Einsatz von Klassiﬁkatoren zum Lernen von Objektbeschreibungen aus hierarchisch partitionierten Bildern Christian Thies, Marcel Schmidt-Borreda und Thomas Lehmann Institut f¨ ur Medizinische Informatik, RWTH Aachen, 52057 Aachen Email: [emailprotected]

Zusammenfassung. Die automatisierte explizite Extraktion von Objekten aus Bildserien, erfordert eine reproduzierbare Beschreibung der entsprechenden Bildregionen. Ein hierarchisches Partitionierungsverfahren zerlegt dazu ein Bild in seine visuell plausiblen Regionen, f¨ ur die dann ein Vektor mit beschreibenden ordinalen Merkmalen berechnet wird. Objektextraktion entspricht damit der Klassiﬁkation entsprechender Merkmalsvektoren, die vom Anwender durch markieren, trainiert werden. Da f¨ ur die Klassiﬁkation das “No-Free-Lunch-Theorem” gilt, m¨ ussen Klassiﬁkator und Merkmalsauswahl f¨ ur jede Dom¨ ane experimentell ermittelt werden. Beim Vergleich des Nearest-Neighbor Klassiﬁkators mit dem Bayes Klassiﬁkator mit Gaußschen Mischverteilungen und der Supportvektormaschine liefert letztere f¨ ur die Handknochenextraktion aus R¨ ontgenaufnahmen das beste Ergebnis.

Einleitung

In der medizinische Bildverarbeitung w¨ achst die Menge der digital erfassten Bilder, sowie die Arten von Objekten, die in ihnen gesucht werden, in einem Maße, dass eine rein manuelle Auswertung mittelfristig nicht mehr m¨oglich sein wird. Daher werden Verfahren immer wichtiger, mit denen a-priori unbekannte Objektbeschreibungen ﬂexibel und eﬃzient formulierbar sind. Hierzu ist in erster Linie eine anwenderorientierte Benutzerschnittstelle erforderlich, in der das Lernen anwendungsangepasst und mit kleinen Stichproben erfolgt. Hierarchische Partitionierungen liefern eine vollst¨andige Bildzerlegung in visuell plausible Regionen, die sich durch einen Merkmalsvektor beschreiben lassen [1]. Um eine solche Region als Objekt zu identiﬁzieren, muss der regionenbeschreibende Merkmalsvektor analysiert werden, was einer Klassiﬁkationsaufgabe entspricht. Mit diesem Ansatz l¨ asst sich die komplexe Formalisierung von a-priori unbekannten Objekten als anwenderorientiertes “Point & Click”- Training realisieren. Problem dabei ist die Auswahl des geeigneten Klassiﬁkators und dessen Integration in einen Workﬂow. In der Literatur existieren zahlreiche Klassiﬁkatoren, sowie Methoden zur Merkmalsreduktion, Sampling der Daten und Parametersuche [2]. Ferner gibt es Verfahren, die Bilder vollst¨ andig und hierarchisch in die visuell plausiblen

C. Thies, M. Schmidt-Borreda und T. Lehmann

Regionen aller Skalen partitionieren [1]. Die Auswahl des Klassiﬁkators h¨angt allerdings von der Struktur des verwendeten Merkmalsraums ab und l¨asst sich nach dem “No-Free-LunchTheorem” nicht generisch beantworten, sondern muss experimentell ermittelt werden [2]. Die mit dem hier vorgestellten Framework durchgef¨ uhrten Experimente dienen dazu, die Abh¨ angigkeit der Klassiﬁkationsergebnisse vom Merkmalsraum, den verwendeten Parametern und dem gew¨ ahlten Klassiﬁkationsparadigma f¨ ur hierarchisch partitionierte Bilder zu verstehen. Auf diese Weise wird untersucht, ob es m¨ oglich ist, Anwenderwissen durch simples Markieren relevanter Beispielregionen ohne formale Interaktion zu lernen, und dieses Wissen automatisiert zu reproduzieren. Dabei werden reale Daten der klinischen Routine verwendet.

Klassiﬁkatoren fu ¨r die Objektsuche

Als Werkzeug dient ein Framework, das den Datenﬂuss der Klassiﬁkation modelliert und so eﬃziente Experimente erm¨ oglicht. Zur Merkmalsreduktion stehen die Lineare Diskriminanz Analyse (LDA), die Hauptkomponenten Analyse (PCA) sowie das auf- und absteigende Greedy-Verfahren zur Verf¨ ugung. Als Klassiﬁkatoren wurden die Supportvektormaschine (SVM), der Bayes Klassiﬁkator mit Gaußschen Mischverteilungen (GMM) und das k-Nearest-Neighbor¨ Ahnlichkeitsmaß (KNN) gew¨ ahlt. SVM und GMM sind modellbasiert, wobei die SVM versucht, strukturelles und empirisches Risiko zu minimieren, w¨ahrend GMM zwar theoretisch durch empirische Risikominimierung immer das beste Resultat liefert, dieses aber mit Overﬁtting einhergeht. Im Vergleich dazu wird der KNN-Klassiﬁkator als modellfreier Ansatz verwendet. Die Tests werden mit unterschiedlich gesampelten Test- und Trainingsdaten durchgef¨ uhrt. Dazu steht das Bootstrapping zur Verf¨ ugung. Beim GMM kann zus¨atzlich Varianz-Pooling verwendet werden. Aus Anwendersicht ist bei der Klassiﬁkation die Zahl der korrekt identiﬁzierten Objekte entscheidend. Dies sowohl im Verh¨altnis zur Gesamtzahl aller R¨ uckgaben (Precision) als auch in Bezug zu allen im Datensatz vorhandenen Werten (Recall). Die Accuracy als Maß der Klassiﬁkationsg¨ ute ist f¨ ur die Objektsuche von untergeordneter Bedeutung. Daher wird das F-Measure als harmonisches Mittel aus Precision und Recall als Qualit¨atsmaß verwendet [3]. 2 · Recall · Precision (1) Recall + Precision Ziel der Experimente ist also die Bestimmung des Klassiﬁkators, der ein maximales F-Measure (Gl. 1) f¨ ur die gegebene Ground-Truth liefert. F-Measure =

Experimente

Die Experimente werden auf 105 zuf¨ allig ausgew¨ahlten Handradiographien aus der klinischen Routine durchgef¨ uhrt. Diese Bilder werden z.B. f¨ ur die vergleichende Auswertung von Handknochen zur Maturit¨atsbestimmung bei Heranwachsenden eingesetzt. In dieser Anwendung w¨ urde eine automatisierte Extraktion der Knochen dem befundenden Radiologen reproduzierbare Abmessungen

Einsatz von Klassiﬁkatoren zum Lernen von Objektbeschreibungen

der Objekte liefern. F¨ ur die Experimente wurden die Mittelhandknochen als Referenzobjekte gew¨ ahlt. Die Bilder werden auf eine maximale Kantendimension von 256 Pixeln unter Beibehaltung der Aspect-Ratio verkleinert. In dieser Gr¨oße sind die gesuchten Handknochen immer noch erkenn- und vermessbar, wobei der Raum der zu durchsuchenden Daten beschr¨ ankt bleibt. Die Partitionierung erfolgt mittels eines hierarchischen Bereichswachstumsverfahrens [1], das die Verschmelzung von Pixeln u ur al¨ ber Regionen hin zum gesamten Bild protokolliert. Dabei werden f¨ le Regionen 38 beschreibende Merkmale (Rundheit, Gr¨oße, Grauwertvarianz,...) berechnet. F¨ ur jedes Bild ergeben sich ca. 2.500 Regionen, was wiederum f¨ ur die gesamte Bildmenge einen 43-dimensionalen Merkmalsraum mit ca. 265.000 Elementen ergibt. Als Ground-Truth wurden manuell 372 Vektoren bestimmt, die die erkennbaren Mittelhandknochen im Datenraum repr¨asentieren. Alle anderen Vektoren fallen in die R¨ uckweisungsklasse. 3.1

Merkmalsauswahl

Im ersten Schritt wird die Abh¨ angigkeit des Ergebnisses von der Merkmalsauswahl und damit ein optimaler Merkmalssatz bestimmt. Die Tests erfolgen dabei u allige Auswahl von 10 Trainings- und 53 Testbilder aus der ¨ ber eine zuf¨ Ground-Truth von 105 Bildern, also etwa 130.000 Merkmalsvektoren und 50 Trainigsobjekte. Greedy-Verfahren. F¨ ur alle drei Klassiﬁkatoren wird jeweils das auf- und absteigende Greedy-Verfahren, zur Merkmalsauswahl angewandt. Die Parametersuche f¨ ur die SVM erfolgt dabei mittels Grid-Suche und Gaußschem Kernel (RBF), f¨ ur den GMM-Klassiﬁkator wird ein Cluster sowie Varianz-Pooling verwendet. Der KNN wird f¨ ur den ersten Nachbarn (1-NN) mit minimaler Euklidischer Norm (L2), Mahalanobis-Distanz (MAH) sowie Cosinusnorm (COS) untersucht. Korrelationsanalyse. Zum Vergleich mit den beiden Greedy-Verfahren werden die PCA und die LDA f¨ ur eine bis 38 Zieldimensionen jeweils f¨ ur alle drei Klassiﬁkatoren getestet. 3.2

Anwendungsszenario

Im zweiten Schritt wird das Anwendungsszenario simuliert. Es entspricht dem Markieren von gesuchten Objekten in einer kleinen und zuf¨alligen Auswahl von Bildern durch den Anwender. Mit den markierten Vektoren wird dann der Klassiﬁkator trainiert, und auf die verbleibende Testmenge angewandt. Es werden 5, 10, 20 und 40 Mittelhandknochen gew¨ ahlt, was 1,2,4, bzw. 8 aus 105 Bildern zum Training entspricht. Um die statistische Verteilung der Auswahl zu simulieren, werden jeweils 20 Samples mittels Bootstrapping gezogen.

C. Thies, M. Schmidt-Borreda und T. Lehmann

Tabelle 1. F-Measure, Precision und Recall f¨ ur das aufsteigende Greedy-Verfahren im Vergleich der untersuchten Klassiﬁkatoren Klassiﬁkator # Merkmale Recall Precision F-Measure 1-NN, L2 7 0,36 0,28 0,31 1-NN, COS 12 0,52 0,28 0,36 1-NN, MAH 4 0,36 0,28 0,32 GMM 12 0,58 0,37 0,45 SVM 14 0,58 0,67 0,62

4 4.1

Ergebnisse Merkmalsauswahl

Greedy-Verfahren aufsteigend. Bei der Berechnung der Merkmalsreduktion mit dem aufsteigenden Greedy-Verfahren war die SVM mit einem F-Measure von 0.62 und 14 Merkmalen der Klassiﬁkator mit der h¨ochsten erzielbaren Genauigkeit w¨ ahrend der modellfreie 1-NN Klassiﬁkator mit der L2-Norm das schlechteste Ergebnis liefert(Tab. 1).Von den 38 verf¨ ugbaren Merkmalen wurden maximal 14 ausgew¨ ahlt. Precision und Recall variieren f¨ ur die einzelnen Klassiﬁkatoren. Greedy-Verfahren absteigend. Mit dem absteigenden Greedy-Verfahren lieferte ebenfalls die SVM bei 32 Merkmalen das beste Ergebnis allerdings nur mit einem F-Measure von 0.46. F¨ ur den 1-NN ergaben sich mit COS bei 12 Merkmalen ein F-Measure von 0.36. PCA. Die zum Vergleich durchgef¨ uhret PCA lieferte mit einem F-Measure von 0.21 ihr bestes Ergebnis bei einer Zieldimension von 26 f¨ ur GMM. Dabei verbessert sich das Ergebnis mit wachsender Zieldimension. LDA. Bei der LDA liefert 1-NN mit der L2 Norm das beste Ergebnis mit einem F-Measure von 0.31 bei einer Zieldimension von 23. Dabei variieren die Werte f¨ ur die Zieldimensionen von 1 bis 38 nicht monoton zwischen 0.2 und 0.3. 4.2

Anwendungsszenario

Die Simulation des Anwendungsszenarios liefert f¨ ur die SVM bei 20 unterschiedlichen Trainings-Samples von jeweils 40 Knochen ein bestes mittleres F-Measure von 0.42 (Tab. 2). Die angegebenen Werte sind dabei jeweils eine Mittelung aus 20 Einzelklassiﬁkationen. Precision und Recall weisen f¨ ur die SVM und den GMM Klassiﬁkator hohe Diﬀerenzen auf, aber man erkennt einen monotonen Anstieg mit wachsender Anzahl Trainingsobjekte Ein Wilcoxon-Test auf dem 5% Niveau best¨ atigt außerdem das signiﬁkant bessere Abschneiden der SVM u ur den Grenzwert von ¨ ber die 20 Samples im Vergleich zu 1-NN und GMM. F¨ 80 Trainingsvektoren ergeben sich keine h¨ oheren Werte.

Diskussion

Die SVM liefert f¨ ur die gestellte Aufgabe das beste Ergebnis hierbei dient die Kombination aus hoher Generalisierungsf¨ ahigkeit bei gleichzeitiger Minimierung

Einsatz von Klassiﬁkatoren zum Lernen von Objektbeschreibungen

Tabelle 2. Ergebnisse der Klassiﬁkatoren bei wachsender Anzahl von Trainigsobjekten. # Trainingobjekte Klassiﬁkator 5 Knochen 1-NN L2-Norm GMM (1) SVM 10 Knochen 1-NN L2-Norm GMM (1) SVM 20 Knochen 1-NN L2-Norm GMM (1) SVM 40 Knochen 1-NN, L2-Norm GMM (1) SVM

Recall Precision F-Measure 0,25 0,26 0,25 0,28 0,26 0,25 0,26 0,29 0,25 0,27 0,28 0,28 0,50 0,22 0,30 0,24 0,44 0,28 0,29 0,29 0,29 0,30 0,30 0,30 0,35 0,43 0,37 0,32 0,33 0,33 0,60 0,28 0,38 0,38 0,52 0,42

des empirischen Risikos dem Ausgleich der extremen Klassenschieﬂage zwischen Objekt- und R¨ uckweisungsklasse. Die Greedy-Heuristik liefert die Merkmalsauswahl, die zum besten Ergebnis f¨ uhrt. Dies liegt an der nichtlinearen Korrelation der Merkmalsverteilungen, die von LDA und PCA nicht modelliert wird. Diese Beobachtungen lassen sich jedoch nach dem “No-Free-Lunch-Theorem” nicht ohne weiteres verallgemeinern. F¨ ur andere Aufgabenstellungen, d.h. Merkmalsr¨ aume, m¨ ussen die Experimente entsprechend wiederholt werden, wozu das vorgestellte Framework ein nachvollziehbares Werkzeug bietet. Die Erkennungsrate f¨ ur die vorgestellte Anwendung der lokalen Bildanalyse ist in jedem Fall verbesserungsf¨ ahig, sie liegt allerdings in einer Gr¨oßenordnung wie sie auch aktuelle und von der Komplexit¨ at vergleichbare Ans¨atze zum globalen inhaltsbasierten Bilddatenbankzugriﬀ liefern [4]. Die Experimente zeigen, das ein “Point & Click”-Training zur Beschreibung und Extraktion a-priori unbekannter Objekte mittels Klassiﬁkation prinzipiell m¨oglich ist. Um die Genauigkeit zu verbessern, m¨ ussen die verwendeten Merkmale und ihre Variationen jedoch weitergehend untersucht werden.

Literaturverzeichnis 1. Beier D, Thies C, G¨ uld MO, Fischer B, Kohnen M, Lehmann TM. Ein lokal¨ adaptives Ahnlichkeitsmaß als Kriterium der hierarchischen Regionenverschmelzung. In: Procs BVM; 2004. p. 100–4. 2. Duda RO, Hart PF, Stork DG. Pattern Classiﬁcation. 2nd edition, Wiley Interscience, New York; 2001. 3. van Rijsbergen CJ. Information Retrieval. 2nd edition, Butterworths, London; 1979. 4. Clough P, M¨ uller H, Sanderson M. The CLEF 2004 Cross-Language Image Retrieval Track. LNCS 2005;3491:597–613.

Rekonstruktion von 4D-CT-Daten aus r¨ aumlich-zeitlichen CT-Segmentfolgen zur Analyse atmungsbedingter Organbewegungen Ren´e Werner1 , Jan Ehrhardt1 , Thorsten Frenzel2 , Dennis S¨ aring1 , Daniel Low3 und Heinz Handels1 1 Institut f¨ ur Medizinische Informatik, Universit¨ atsklinikum Hamburg-Eppendorf, 20246 Hamburg 2 Hermann-Holthusen-Institut f¨ ur Strahlentherapie, Allgemeines Krankenhaus St. Georg, 20099 Hamburg 3 Mallinckrodt Institute of Radiology, Washington University, St. Louis, USA Email: [emailprotected]

Zusammenfassung. Atmungsbedingte Organbewegungen stellen eines der Hauptprobleme der Strahlentherapie thorakaler und abdominaler Tumoren dar. Die Entwicklung von L¨ osungsans¨ atzen bedarf der Analyse des r¨ aumlich-zeitlichen Verhaltens der strahlentherapeutisch relevanten Volumina, z.B. auf Basis von 4D-CT-Daten. Moderne CT-Scanner gestatten allerdings lediglich die simultane Aufnahme einer begrenzten Anzahl benachbarter K¨ orperschichten. Um dennoch Bewegungen gr¨ oßerer Volumina untersuchen zu k¨ onnen, werden durch wiederholtes Scannen der entsprechenden anatomischen Segmente r¨ aumlich-zeitliche CT-Segmentfolgen generiert. Es kann jedoch nicht sichergestellt werden, dass f¨ ur die Scans der unterschiedlichen anatomischen Segmente die Zeitpunkte relativ zum Atemzyklus des Patienten einander entsprechen. Werden 3DCT-Daten zu einem vorgegebenen Zeitpunkt des Atemzyklus ausschließlich aus den aufgezeichneten Datensegmenten zusammengesetzt, treten folglich Bewegungsartefakte an den Segmentgrenzen auf. Dieser Beitrag pr¨ asentiert ein Verfahren zur Rekonstruktion von (3D+t)-CT-Daten (4DCT-Daten), das zur Reduktion der Artefakte f¨ uhrt. Hierzu wird der Optische Fluss zwischen den aufgezeichneten Datensegmenten bestimmt und zur Interpolation von Datensegmenten f¨ ur vorgegebene Zeitpunkte des Atemzyklus eingesetzt. Die rekonstruierten 4D-CT-Daten bilden die Grundlage der Analyse und Visualisierung der Bewegungen der dargestellten Bildstrukturen (Lungenﬂ¨ ugel, Bronchialbaum, Lungentumoren).

Problemstellung

Atmungsbedingte Organbewegungen stellen eines der Hauptprobleme im Kontext der Strahlentherapie thorakaler und abdominaler Tumoren dar. Derzeitig werden m¨ ogliche Bewegungen des Tumors durch Addition statischer Sicherheitss¨ aume zu dem Volumen ber¨ ucksichtigt, das auf Basis eines 3D-CT-Datensatzes als tats¨ achlich bzw. sehr wahrscheinlich tumorhaltig identiﬁziert wird.

Rekonstruktion von 4D-CT-Daten aus r¨ aumlich-zeitlichen CT-Segmentfolgen

Aufgrund der großen Amplituden atmungsbedingter Organbewegungen kann ein solches Vorgehen stark zu Lasten des gesunden Gewebes gehen. Eine explizite Ber¨ ucksichtigung der atmungsbedingten Dynamik ist demnach erstrebenswert, bedarf allerdings einer detaillierten Kenntnis bzw. der Analyse des r¨aumlichzeitlichen Verhaltens der strahlentherapeutisch relevanten Volumina. Zur Durchf¨ uhrung entsprechender Analysen werden mittels eines MehrzeilenCT r¨ aumlich-zeitliche CT-Segmentfolgen aufgenommen (siehe [1]): Die derzeitige Mehrzeilen-Computertomographie gestattet lediglich die simultane Aufnahme einer begrenzten Anzahl benachbarter K¨ orperschichten. Um dennoch Bewegungen gr¨ oßerer Volumina untersuchen zu k¨ onnen, werden die das zu untersuchende Volumen abdeckenden Couchpositionen wiederholt gescannt. Zeitgleich wird mittels Spirometrie das Atemvolumen des Patienten gemessen. Hierauf basierend werden den aufgezeichneten Datensegmenten das zugeh¨orige Atemvolumen und die Atemphase (Ein- oder Ausatmung) zugeordnet. Die Gesamtheit der Datensegmente stellt somit eine r¨ aumlich-zeitliche CT-Segmentfolge dar. Gesucht ist nun ein Verfahren, das es erlaubt, auf Basis der CT-Segmentfolgen (3D+t)-CT-Daten (d.h. 4D-CT-Daten) zu rekonstruieren. Als problematisch erweist sich hierbei, dass bei Aufzeichnung der Segmentfolgen, insbesondere bei freier Atmung des Patienten, f¨ ur unterschiedliche Couchpositionen die Datensegmente zumeist zu verschiedenen Zeitpunkten des Atemzyklus erfasst werden. Setzt man die 3D-CT-Daten des gesuchten 4D-Datensatzes aus den aufgezeichneten Datensegmenten zusammen, treten Bewegungsartefakte auf, die sich durch Konturspr¨ unge an den Segmentgrenzen manifestieren (siehe Abb. 1, oben). Eine Minimierung der Bewegungsartefakte bildet die Grundlage f¨ ur eine computergest¨ utzte Analyse und Visualisierung der atmungsbedingten Organbewegungen.

Methodik

Durch Vorgabe einer sortierten Sequenz von Zeitpunkten des Atemzyklus (= ˆ Wertepaaren, je bestehend aus Atemvolumen und Atemphase) wird die zeitliche Auﬂ¨ osung des zu generierenden (3D+t)-CT-Datensatzes festgelegt. Zu jedem der Zeitpunkte ist ein korrespondierender 3D-CT-Datensatz zu rekonstruieren. 2.1

Rekonstruktion durch N¨ achster-Nachbar-Interpolation

Bislang eingesetzte Verfahren zur Rekonstruktion von (3D+t)-CT-Daten aus bei freier Atmung des Patienten aufgenommenen, r¨aumlich-zeitlichen CT-Segmentfolgen beruhen auf dem Prinzip der N¨achster-Nachbar-Interpolation: Um einen 3D-CT-Datensatz zu einem vorgegebenen Zeitpunkt des Atemzyklus des Patienten zu rekonstruieren, wird auf Basis der spirometrischen Daten zu jeder Couchposition dasjenige aufgezeichnete Datensegment ermittelt, dessen zugeordneter Zeitpunkt dem Vorgegebenen am n¨ achsten kommt. Der 3D-CT-Datensatz wird dann generiert, indem diese Datensegmente entsprechend ihrer r¨aumlichen Positionen zusammengesetzt werden. Aufgrund der dargestellten Problematik f¨ uhrt dieses Vorgehen zu den beschriebenen Bewegungsartefakten.

2.2

R. Werner et al.

Rekonstruktion durch Interpolation unter Verwendung des Optischen Flusses

Um eine Minimierung dieser Bewegungsartefakte zu erzielen, wurde ein neues Rekonstruktionsverfahren entwickelt und umgesetzt. Zur Rekonstruktion eines 3D-Datensatzes zu einem gegebenen Atemvolumen werden hierbei f¨ ur jede zu ber¨ ucksichtigende Couchposition Datensegmente interpoliert, die dem zur Rekonstruktion vorgegebenen Atemvolumen entsprechen. Hierzu werden auf Basis der spirometrischen Informationen f¨ ur die einzelnen Couchpositionen je die beiden aufgezeichneten Datensegmente ermittelt, die dem vorgegebenen Zeitpunkt am n¨ achsten kommen (zeitlich vorhergehend, zeitlich nachfolgend). Zwischen diesen wird durch Anwendung der d¨amonenbasierten Registrierung nach Thirion [2] der Optische Fluss [3] gesch¨ atzt. Das resultierende Geschwindigkeitsfeld repr¨ asentiert die auftretenden Bewegungen innerhalb des zeitlichen Intervalls, das durch die den Datensegmenten zugeordneten Zeitpunkte deﬁniert wird, und wird in Analogie zu [4] zur Erzeugung eines interpolierten Datensegmentes zum vorgegebenen Zeitpunkt herangezogen. Durch auftretende Bewegungen u ¨ ber die Segmentgrenzen hinaus k¨ onnen hierbei allerdings Registrierungsartefakte an den Segmentgrenzen auftreten. Diesen wird entgegengewirkt, indem die origin¨ar zu verarbeitenden Datensegmente um die zeitlich ann¨ahernd korrespondierenden Datensegmente benachbarter Couchpositionen erweitert werden. Die beschriebene Methodik wird auf diese erweiterten Datensegmente“ angewendet. Indem ” somit w¨ ahrend des Registrierungsprozesses Voxelkorrespondenzen u ¨ber die Grenzen der origin¨ aren Datensegmente hinaus deﬁniert werden k¨onnen, werden die Artefakte reduziert. 2.3

Evaluation der Rekonstruktionsverfahren

Die Rekonstruktionsverfahren werden evaluiert, indem die innerhalb der generierten 4D-CT-Datens¨ atze auftretenden Spr¨ unge quantiﬁziert werden. Hierzu ¨ wird das Verh¨ altnis der mittleren Ahnlichkeit der Schichten eines Segment¨ uberganges (unterste Schicht des oberen Segmentes, oberste Schicht des unteren Seg¨ mentes) und der mittleren Ahnlichkeit der innerhalb eines Segmentes benachbarten Schichten betrachtet. Gemittelt wird u uber¨ber alle entsprechenden Schicht¨ ¨ g¨ange s¨ amtlicher 3D-Datens¨ atze des betrachteten 4D-Datensatzes. Als Ahnlichkeitsmaß wird die SSD-Metrik (SSD: Summed Squared Deviation) eingesetzt. 2.4

Analyse und Visualisierung von Organbewegungen auf Basis der rekonstruierten 4D-CT-Daten

Die rekonstruierten 4D-CT-Datens¨ atze werden zur Analyse und Visualisierung atmungsbedingter Organbewegungen eingesetzt. Hierzu werden die interessierenden Strukturen (hier: Lungenﬂ¨ ugel, Bronchialbaum, Lungentumoren) segmenardaten wird dann f¨ ur die einzelnen 3D-CTtiert. Unter Verwendung der Bin¨ Datens¨ atze u uft, ob das Volumen der in der Lunge enthaltenen Luft (die ¨ berpr¨ Bestimmung des Volumens erfolgt in Anlehnung an [5]) und die spirometrisch

Rekonstruktion von 4D-CT-Daten aus r¨ aumlich-zeitlichen CT-Segmentfolgen

Abb. 1. Schnitt eines rekonstruierten 3D-CT-Datensatzes zu einem Atemvolumen von 800 ml (Ausatmung). Oben: Rekonstruktion durch N¨ achster-Nachbar-Interpolation. Unten: Rekonstruktion durch Interpolation unter Verwendung des optischen Flusses.

gemessenen Atemvolumenwerte korrelieren. Weitergehend werden organspeziﬁsche Landmarken verfolgt und deren Trajektorien analysiert. 4D-Darstellungen der Bewegungen der Strukturen werden erzeugt, indem zeitliche Sequenzen von Oberﬂ¨ achenmodellen generiert werden.

Ergebnisse

Die Rekonstruktionsverfahren wurden jeweils auf vier r¨aumlich-zeitliche CTSegmentfolgen (Lungentumor-Patienten) angewendet, aufgenommen mittels eines 12-Zeilen-CT (r¨ aumliche Auﬂ¨ osung: 0.9 x 0.9 x 1.5 mm bei 512 x 512 x 12 Voxel pro Couchposition; Datens¨ atze umfassen 16-19 Couchpositionen). F¨ ur jede Couchposition liegen 15 Datensegmente vor, die jeweils im Abstand von 0.75 s aufgezeichnet wurden. Die Datens¨ atze umfassen den Thorax sowie das obere Abdomen und bestehen je aus mehr als 2880 transversalen Schichtbildern. Die visuelle Pr¨ ufung der Datens¨ atze, die aus der Rekonstruktion durch Interpolation unter Verwendung des optischen Flusses resultieren, ergibt im Vergleich

R. Werner et al.

zu den Daten, die mittels der N¨ achster-Nachbar-Interpolation erzeugt sind, eine deutliche Reduktion der Bewegungsartefakte (siehe Abb. 1). Anhand der zur Evaluation der Verfahren eingef¨ uhrten Kenngr¨ oße werden f¨ ur die vier Patientendatens¨ atze Reduktionen der Spr¨ unge um 32.7 %, 23.3 %, 24.6 % und 27.3 % erzielt. Die Bewegungsanalyse der rekonstruierten 4D-CT-Daten erbringt f¨ ur den Zusammenhang zwischen dem spirometrisch gegebenen Atemvolumen und dem in der Lunge enthaltenen Luftvolumen f¨ ur die Patienten jeweils einen Korrelationskoeﬃzienten von 0.99. Die ermittelten Bewegungsamplituden organspeziﬁscher Landmarken stimmen mit verf¨ ugbaren Literaturwerten [6] u ¨ berein (z.B. Bifurkation und Lungenbasis: 1-2 cm). F¨ ur Lungentumoren wird eine starke Abh¨ angigkeit des Bewegungsausmaßes von der Lage des Tumors beobachtet.

Diskussion

Die Rekonstruktion von 4D-CT-Daten aus r¨ aumlich-zeitlichen CT-Segmentfolgen durch Interpolation von Datensegmenten unter Verwendung des optischen Flusses erlaubt es, zu vorgegebenen Zeitpunkten des Atemzyklus 3D-CT-Datens¨atze zu generieren. Insbesondere bei Anwendung auf Segmentfolgen, die bei freier Atmung des Patienten aufgenommen werden, zeichnet sich das Verfahren im Vergleich zu der bisher eingesetzten N¨ achster-Nachbar-Interpolation durch eine ¨ deutliche Reduktion atmungsbedingter Bewegungsartefakte an den Uberg¨ angen zwischen r¨ aumlich benachbarten Datensegmenten der 3D-CT-Datens¨atze aus. Die resultierenden 4D-CT-Daten sollen zu weiterf¨ uhrenden computergest¨ utzten Analysen atmungsbedingter Bewegungen anatomischer Strukturen verwendet werden. Der Schwerpunkt wird hierbei auf dem Auﬃnden von Korrelationen zwischen Bewegungen k¨ orperinterner Strukturen und Bewegungen der K¨orperoberﬂ¨ ache liegen, um z.B. Regionen der K¨ orperoberﬂ¨ache zu identiﬁzieren, die ein nicht invasives Tracking der Bewegungen strahlentherapeutisch relevanter Volumina erlauben.

Literaturverzeichnis 1. Low DA, Nystrom M, Kalinin E, et al. A method for the reconstruction of fourdimensional synchronized CT scans acquired during free breathing. Med Phys 2003;30(6):1254–1263. 2. Thirion JP. Image matching as a diﬀusion process: an analogy with Maxwell’s demons. Med Image Anal 1998;2(3):243–260. 3. Horn BKP, Schunck BG. Determining optical ﬂow. Artif Intell 1981;17:185–203. 4. Ehrhardt J, S¨ aring D, Handels H. Optical Flow-based Interpolation of Temporal Image Sequences. In: Procs BVM; 2006, accepted. 5. Lu W, Parikh PJ, Naqa IMEl, et al. Quantitation of the reconstruction quality of a four-dimensional computed tomography process for lung cancer patients. Med Phys 2005;32(4):890–901. 6. Langen KM, Jones DT. Organ motion and its management. Int J Radiat Oncol Biol Phys 2001;50(1):265–278.

3D-Erkennung, Analyse und Visualisierung pleuraler Verdickungen in CT-Daten Patrick J¨ ager1, Stefan Vogel3 , Achim Knepper1 , Thomas Kraus2 und Til Aach1 1

Lehrstuhl f¨ ur Bildverarbeitung, RWTH Aachen Templergraben 55, 52056 Aachen 2 Institut f¨ ur Arbeitsmedizin, Universit¨ atsklinikum der RWTH Aachen, 52057 Aachen 3 Lehrstuhl f¨ ur medizinische Informationstechnik, RWTH Aachen Templergraben 55, 52056 Aachen Email: [emailprotected] Internet: www.lfb.rwth-aachen.de/mesotheliom

Zusammenfassung. Aufbauend auf fr¨ uheren Arbeiten wurde ein Bildverarbeitungssystem entwickelt, welches pleurale Verdickungen automatisch lokalisiert und visualisiert. Es liefert reproduzierbare, quantitative Daten, die eine genauere Beobachtung der Verdickungen erm¨ oglichen als die konventionelle Befundungsmethode, und reduziert den f¨ ur die Befundung n¨ otigen Zeitaufwand. Die automatische Detektion ﬁndet innerhalb eines zweistuﬁgen Algorithmus statt, der zuerst aus allen Schichten des Datensatzes die Pleurakonturen extrahiert und darauf aufbauend die ¨ Verdickungen in den Pleurakonturen detektiert. Da die Anderung der Form einer Verdickung ein wichtiges Kriterium bei der Entscheidung ist, ob eine Verdickung entartet, wurde eine M¨ oglichkeit zur Visualisierung der Verdickungen und der Lungenﬂ¨ ugel implementiert. Diese k¨ onnen nun aus allen Perspektiven betrachtet werden. Unterschiede in den Verdickungen zweier aufeinanderfolgender Scans k¨ onnen so erkannt werden.

Problemstellung

Asbest ist ein Sammelbegriﬀ f¨ ur nat¨ urlich vorkommende, faserf¨ormige, mineralische Silikatmaterialien, die dank ihrer guten thermischen, mechanischen und isolatorischen Eigenschaften in vielen Anwendungsbereichen eingesetzt wurden. Asbest ist eine gesicherte humankanzerogene Substanz und kann Lungen- und Pleuraﬁbrosen, Lungen- und Kehlkopfkrebs sowie maligne Mesotheliome verursachen. Die Asbestfasern gelangen u ¨ber die Atemluft in die Lunge und lagern sich an der Pleura ab, wo sich gutartige und b¨osartige Ver¨anderungen ausbilden k¨ onnen. Die Latenzzeit zwischen der Asbestexposition und dem Entstehen eines Mesothelioms kann 10 bis 65 Jahre betragen [1, 2, 3]. Der H¨ohepunkt des Asbesteinsatzes wurde in Deutschland Ende der 70er Jahre erreicht. Daher geht man von einem weiteren Anstieg der Neuerkrankungen auf u ¨ ber 1000 F¨alle pro Jahr aus. Bei den Todesf¨ allen Berufserkrankter ist Asbest die weitaus h¨auﬁgste

P. J¨ ager et al.

Todesursache und stellt damit in der Arbeitsmedizin ein großes Problem dar. In anderen europ¨ aischen L¨ andern sind die Fallzahlen noch deutlich h¨oher [4]. Daher werden in den kommenden 30 Jahren in Westeuropa etwa 250.000 Mesotheliomtote prognostiziert [5].

Stand der Forschung

Eine Heilung des Mesothelioms ist bisher nur in Einzelf¨allen m¨oglich, wenn eine ¨ Fr¨ uherkennung erfolgte. Zurzeit betr¨ agt die mittlere Uberlebenszeit nach Dia¨ gnosestellung ca. 9 Monate. In dem Bestreben, die Uberlebenszeit Betroﬀener zu erh¨ ohen, kommt der Fr¨ uherkennung des malignen Tumors eine entscheidende Rolle zu. Eine Wachstumstendenz pleuraler Verdickungen im zeitlichen Verlauf stellt ein wesentliches Kriterium zur Diﬀerenzierung zwischen gutartigen und b¨ osartigen Ver¨ anderungen dar. Zur Diagnosestellung von gutartigen und b¨ osartigen Ver¨ anderungen wird die Pleura mittels CT auf Verdickungen untersucht, die dabei vom befundenden Arzt nach Gr¨oße und Form in Klassen eingeteilt werden. Die Ergebnisse werden in einen standardisierten Erfassungsbogen eingetragen. Dies ist zum einen mit sehr hohem zeitlichen Aufwand verbunden, zum anderen unterliegen die Auswertungsergebnisse aufgrund der subjektiven visuellen Beurteilung hohen Schwankungen. Intra- und Interreadervariabilit¨ atsuntersuchungen zeigen, dass eine genaue und reproduzierbare Aussage u ¨ ber die quantitative Ausdehnung einer Pleuraverdickung schwer zu treﬀen ist. Heute zur Befundung eingesetzte, kommerziell erh¨ altliche Workstations werden momentan lediglich als Werkzeug zur Visualisierung der CT-Scans eingesetzt und erleichtern es dem Befunder, die Ergebnisse in Formulare einzutragen. In [6] wird von einem System berichtet, das semiautomatisch zur eindimensionalen Vermessung von Mesotheliomen in der Lage ist. Systeme, welche pleurale Verdickungen automatisiert erfassen und vermessen, sind zurzeit noch nicht verf¨ ugbar.

Wesentlicher Fortschritt durch den Beitrag

Die Verwendung eines vollautomatischen Bildverarbeitungssystems bei der Befunderhebung hat zum einen das Ziel, die zeitliche Belastung des befundenden Arztes zu verringern. Zum anderen soll sie zu exakten, reproduzierbaren Messwerten f¨ uhren, die eine einfache, m¨ oglichst objektive Verlaufsdokumentation der pleuralen Verdickungen erm¨ oglichen. Die vom System dargestellten Ergebnisse k¨ onnen dabei jederzeit vom untersuchenden Arzt korrigiert werden. So k¨ onnen Ver¨ anderungen der Verdickungen, die auf die Bildung eines malignen Tumors hinweisen, fr¨ uher erkannt und schon im Fr¨ uhstadium geeignete Therapiemaßnahmen ergriﬀen werden, welche sowohl die Lebenserwartung als auch die Lebensqualit¨ at des Patienten erh¨ ohen.

3D-Erkennung, Analyse und Visualisierung pleuraler Verdickungen

Abb. 1. Extraktion und Visualisierung der Pleurakonturen.

(a) Alle im Bild gefundenen Konturen

(b) Gefundene Pleurakontur, Bronchialbaum herausgeschnitten

eines

Methoden

Im CT-Scan werden zuerst die Pleurakonturen extrahiert, daraufhin nach vordeﬁnierten Kriterien die Verdickungen lokalisiert und deren charakteristische Abmessungen bestimmt. Unser Verfahren arbeitet wie folgt: Der Anwender l¨adt zu Beginn einen CT-Scan zur Befundung in den Arbeitsbereich. Per Knopfdruck startet er die Detektion der Pleurakonturen. Dabei wird jedes DICOM-Bild mit einer Schwellwert-Operation (-550 HU) [7] in ein Bin¨arbild umgewandelt. In diesem Bin¨ arbild werden mittels eines Konturﬁndungsalgorithmus [8] alle geschlossenen Konturen gefunden. Diese werden nach vordeﬁnierten Kriterien sortiert, so dass am Ende die Konturen der Pleura u ¨brig bleiben. Da im Bereich des Bronchialbaumes und mediastinaler Strukturen keine Verdickungen zu erwarten sind, werden diese Bereiche im Scan nicht betrachtet. Mit Hilfe der so gefundenen Pleurakonturen lassen sich die Lungenﬂ¨ ugel in einem gesonderten Fenster in einer 3D-Ansicht visualisieren und von allen Seiten betrachten. Dabei sind die Bereiche der Verdickungen gut erkennbar (Abb. 1c). Mit einem weiteren Knopfdruck startet der Anwender den Algorithmus zum Auﬃnden der Verdickungen. Verdickungen stellen sich immer als nichtkonvexe Bereiche in der konvex modellierten Kontur der Lunge dar. Daher werden in den gefundenen Konturen der Pleura alle Punkte identiﬁziert, welche die Pleurakontur vollst¨ andig umschließen, wenn sie mit Geraden verbunden werden. Aus je zweien dieser Punkte wird nun mit Hilfe ihrer Verbindungsgerade und dem zwischen ihnen liegenden Konturst¨ uck eine neue, kleine Kontur erstellt (Abb. 2). Diese neuen Konturen werden als potentielle Verdickungen angesehen und nach drei Kriterien beurteilt: Sie werden nur weiter betrachtet, wenn sie mindestens ein Pixel umschließen, sich mindestens eine Kontur im selben Bereich des Schichtbildes dar¨ uber oder darunter beﬁndet und sich in ihrem Innern ein bestimmter Anteil an Pixeln beﬁndet, deren Hounsﬁeld-Werte denen von Tumorgewebe (gr¨ oßer 40 HU) entsprechen. Ist eine dieser Bedingungen nicht erf¨ ullt, wird die Kontur verworfen. Die gefundenen Verdickungen (Abb. 3) werden mit

P. J¨ ager et al. Abb. 2. Erstellen der Diﬀerenzkonturen

(a) Knotenpunkte der konvexen H¨ ulle

¨ (b) Uberlagerung von Kontur und konvexer H¨ ulle

der

den sich im selben Bereich beﬁndenden Konturen aus den Schichten dar¨ uber und darunter vervollst¨ andigt. Zur Bestimmung des idealtypischen Verlaufs der Pleura wird ein thin plate model“ verwendet, das mit Hilfe einiger Pixel außerhalb ” des Beginns der Nichtkonvexit¨ aten eingepasst wird. Dadurch wird zum einen die zur Volumenbestimmung n¨ otige R¨ uckwand“ der Verdickung realistisch nachge” bildet und zum anderen der Anfangs- und Endpunkt der Verdickungen in jeder Schicht unter Hinzunahme der 3D-Information optimal bestimmt [9]. Die so gefundenen Verdickungen k¨ onnen ¨ ahnlich wie die Lungen in einem eigenen Fenster in 3D visualisiert und von allen Seiten betrachtet werden. Mittels eines Justage-Verfahrens [7] werden charakteristische Punkte des K¨orpers und ihre Lage relativ zum CT ermittelt. Dies erm¨oglicht eine Koregistrierung der Verdickungen und damit die automatisierte Erstellung einer Verlaufsdokumentation f¨ ur jede Verdickung. Die Ergebnisse der Verarbeitung (Konturen, charakteristische Werte der Verdickungen) werden bei deutlich erh¨ohter Sensitivit¨at des Untersuchungsverfahrens gespeichert und stehen so jederzeit zum Vergleich mit neuen Aufnahmen des Patienten zur Verf¨ ugung.

Ergebnisse

Die Algorithmen zum Auﬃnden der Pleurakonturen, der Verdickungen und die Funktionalit¨ at zur 3D-Visualisierung sowohl der Lungen als auch der Verdickungen funktionierten auf dem zur Verf¨ ugung stehenden Bildmaterial zuverl¨assig. Als n¨ achster Schritt steht eine Validierung der Ergebnisse des Systems in der Praxis im Vergleich zu den Ergebnissen erfahrener Befunder an, ebenso eine Beurteilung der Anwendbarkeit durch Befunder im t¨aglichen Betrieb.

Ausblick

Um das System im klinischen Routinebetrieb einsetzen zu k¨onnen, muss eine Validierung der gelieferten Ergebnisse stattﬁnden. Um dem Ziel eines vollautoma-

3D-Erkennung, Analyse und Visualisierung pleuraler Verdickungen

Abb. 3. Verdickungen - rot: HU > 10, blau: HU < 10

(a) Spindelf¨ ormige ckung

Verdi-

(b) Tafelbergf¨ ormige Verdickung

tischen Systems zur Erfassung und Vermessung pleuraler Verdickungen n¨aher zu kommen, ist neben einer Integration der entwickelten Algorithmen in die schon bestehende Benutzeroberﬂ¨ ache eine Automatisierung der Koregistrierung und die Implementierung einer Verlaufsdokumentation n¨otig.

Literaturverzeichnis 1. Kraus T, Raithel HJ. Fr¨ uhdiagnostik asbeststaubverursachter Erkrankungen. Hauptverband der gewerbl. Berufsgenossenschaften, Sankt Augustin; 1998. 2. R¨ uhle KH. Pleurale Erkrankungen: Diagnostik und Therapie. Kohlhammer, Stuttgart; 1997. 3. Sugarbaker DJ, Norberto JJ, Bueno R. Current Therapy of Mesothelioma. Cancer Control Journal 1997;4(4):326–334. 4. Deriot G, Godefroy JP. Rapport d’information fait au nom de la mission commune d’information sur le bilan et les cons´equences de la contamination par l’amiante. Senat de la republique francaise; 2005. 5. Dami˜ ao. Bericht u ur eine Richtlinie des Europ¨ aischen Parlaments ¨ ber den Vorschlag f¨ ¨ und des Rates zur Anderung der Richtlinie 83/477/EWG des Rates u ¨ber den Schutz der Arbeitnehmer gegen Gef¨ ahrdung durch Asbest am Arbeitsplatz. Europ¨ aisches Parlament; 2002. 6. Armato SG, Oxnard GR, et al. Evaluation of Semiautomated Measurements of Mesothelioma Tumor Thickness on CT Scans. Academic Radiology 2005;12(10):1301– 1309. 7. Vogel S. Semiautomatische Segmentierung, quantitative Vermessung und Verlaufsdokumentation von Pleuramesotheliomen in Spiral-CT-Bildsequenzen. Lehrstuhl f¨ ur Messtechnik und Bildverarbeitung, RWTH Aachen; 2003. 8. Suzuki S, Abe K. Topological Structural Analysis of Digital Binary Images by Border Following. CVGIP 1985;30(1):32–46. 9. J¨ ager P. 3D-Analyse und Visualisierung pleuraler Verdickungen in Spiral-CT-Daten. Lehrstuhl f¨ ur Bildverarbeitung, RWTH Aachen; 1997.

LAST Filter for Artifact-Free Noise Reduction of Fluoroscopic Sequences in Real-Time Marc Hensel1 , Thomas Pralow2 and Rolf-Rainer Grigat1 1 2

Vision Systems, TU Hamburg-Harburg, 21079 Hamburg General X-Ray, Philips Medical Systems, 22335 Hamburg Email: [emailprotected]

Abstract. We present a spatio-temporal ﬁlter for real-time noise reduction of strongly corrupted X-ray image sequences. It possesses eﬃcient noise reduction while, at the same time, preventing typical artifacts of state-of-the-art methods. Decisive for these features are, in particular, innovative motion detection as well as noise-adaptive ﬁlter parametrization. Motion detection based on twofold signed binarization proved to be a powerful method for pixelwise separation of motion and strong noise. Drawbacks of threshold determination by Euler curve analysis as applied previously were eliminated by integration of signal-dependent noise estimation.

Introduction

Real-time noise reduction, enhancement, and visualization of X-ray image sequences is the basis of several established clinical applications. For instance, in angiography, medical instruments are navigated and the distribution processes of contrast agents are observed. As medical staﬀ and patients are exposed over a long period, very low X-ray doses are applied in these examinations. Consequently, images exhibit severe noise. The advance of digital technologies and powerful computers enables us to apply image processing methods of increasing complexity. State-of-the-art in real-time noise reduction is spatio-temporal ﬁltering. Temporal ﬁltering is well suited for noise reduction of static image regions. Motion, however, generates perturbing artifacts. Hence, moving structures should be detected and temporal ﬁltering reduced at these structures. Motion detection is complicated signiﬁcantly by severe signal-dependent noise in the image data at hand. In last year’s BVM workshop our groups have presented motion detection based on independent binarization of positive and negative values in temporal diﬀerence images [1]. In this paper we report on the development status of the enhanced noise reduction method, which we named – for lack of better ideas – LAST (Local Adaptive Spatio-Temporal ) ﬁlter. Following a brief description of the ﬁlter, we focus on advances in motion detection and ﬁlter parametrization. In particular, we replaced blockwise processing based on Euler curves by a pixelwise motion probability that is motivated by an appropriate noise model.

LAST Filter for Artifact-Free Noise Reduction of Fluoroscopic Sequences

Fig. 1. Temporal recursive ﬁlter structure (α + β + γ = 1)

State-of-the-Art

Typically, temporal recursive ﬁlters [2] are applied for real-time noise reduction of X-ray image sequences exhibiting severe noise [3, 4]. In static image regions, temporal ﬁltering signiﬁcantly improves image quality and preserves ﬁne structures far better than spatial ﬁlters of comparable complexity. However, motion generates distracting artifacts (motion shadows). Motion artifacts can be prevented by detecting motion and reducing the strength of temporal ﬁltering in favor of spatial ﬁltering at these locations. As motion generates strong signal in temporal diﬀerence images, motion detection is commonly based on the current and the last ﬁltered image [5]. However, pixelbased motion detection without consideration of the neighborhood yields artifacts as frayed edges or salt-and-pepper noise in very noisy images. Aufrichtig and Wilson [4] assume motion artifacts in non-cardiac angiography to be caused predominantly by long, thin objects like catheters or guide-wires rather than organic structures. They use templates of small line segments to identify these objects and control spatio-temporal ﬁltering. However, at least in angiography, large vessels brought out by contrast agents and, in general, bone structures are likely to produce artifacts as well. Recently, we have presented a method to detect strong noise and motion by independent binarization and postprocessing of positive and negative gray value diﬀerences [1]. For binarization, images are divided into blocks and a threshold for each block is determined in a noise-adaptive way by Euler curve analysis. Problems are, in particular, blocking artifacts and the complexity of Euler curve generation and analysis.

Methods

The ﬁlter is composed of a recursive spatio-temporal ﬁlter structure (Fig. 1) and a control unit determining local ﬁlter parameters. Latter includes a shutter, noise estimation, motion detection (twofold signed binarization and binary postprocessing) and automatic parametrization. The shutter determines the region of X-ray (ROX). Outside this region no motion is expected. Thus, maximal temporal ﬁltering is applied. The integrated noise estimator is subject to a further paper presented by our groups [6].

M. Hensel, T. Pralow and R.-R. Grigat

Main features of the proposed motion detection and ﬁlter control are: 1. Twofold Signed Binarization: Independent binarization of positive and negative temporal gray value diﬀerences generates two binary images. 2. Postprocessing: By means of morphological operations, both binary images are divided into two masks representing motion and strong noise, each. Positive and negative motion masks are combined to one binary image. The proceeding for the noise masks is analogous. 3. Local adaptive control of spatial and temporal ﬁlter strengths via motion and noise masks as well as a motion probability based on local noise estimation. 3.1

Twofold Signed Binarization and Postprocessing

By independent binarization of positive and negative gray value diﬀerences, motion forms connected objects in the binary images while strong noise above the threshold does not form objects of signiﬁcant size. Consequently, motion and noise can be separated into binary motion and noise masks bmotion and bnoise with bmotion ∩ bnoise = ∅ by simple morphological operations [1]. For each pixel, thresholds are determined based on an estimation of the local signal-dependent additive noise [6]. For this purpose, gray value diﬀerences g∆ (x, t) = g(x, t)−gst (x, t−1) of the current (g) and the last spatio-temporal ﬁltered (gst ) image are normed to the signal-dependent standard deviation σηe (x, t) of the estimated noise ηe : |g∆ (x, t)| q(x, t) = (1) σηe (x, t) Absolute diﬀerences above three times the estimated noise level (q > 3) are interpreted as being caused by structure, i.e. motion. Further, by applying a sigmoid function, the normed diﬀerences are mapped to a continuous motion probability: c2 prob(q) = c1 + (2) 1 + exp c3c−q 4 Extremely noisy pixels wrongly classiﬁed as motion are not problematic as these are detected and collected in the noise mask in postprocessing. 3.2

Control

The overall goal of ﬁlter parametrization is to maximize the temporal ﬁlter strength, i.e. parameter β, while, at the same time, preventing motion artifacts. Pixels identiﬁed as representing motion are predominantly ﬁltered spatially, noise pixels are mainly ﬁltered temporally. The spatio-temporal weighting w( · ) of the remainder is carried out on basis of the calculated motion probability: ⎧ : bmotion (x, t) = 1 ⎨ βmotion (3) β(x, t) = βnoise : bnoise (x, t) = 1 ⎩ w prob q(x, t) : else

LAST Filter for Artifact-Free Noise Reduction of Fluoroscopic Sequences

Table 1. Coeﬃcients a1 and a2 for correlated αβγ-ﬁltering Filter Binomial 3 × 3 Binomial 5 × 5 Binomial 7 × 7

a2 0.6406 0.7935 0.8556

a1 1.5 1.7188 1.8047

The remaining parameters α and γ = 1 − (α + β) are optimized in such way to maximize noise reduction performance for given β. For observed medical X-ray images g, the true signal s being corrupted by additive zero-mean uncorrelated Gaussian noise η with signal-dependent variance, i.e. g = s + η(s), is a valid assumption [6]. The signal is modeled to be approximately hom*ogeneous in a small neighborhood. Furthermore, in practice, the spatial ﬁlter is a (linear) binomial ﬁlter h of small size (3 × 3 to 7 × 7 pixels). Assume the neighborhood of a pixel to be corrupted with noise of standard deviation σ0 . Then the spatio-temporal ﬁltered image exhibits noise of strength

a2 2 α + a1 α + (β − 1) · σ0 σst (t) = β 2t + (β 2t − 1) (4) β−1 with: a1 = 2(1 − h0,0 ) a2 =

(5) [1 − δ(i)δ(j)] h2i,j + (1 − h0,0 )2

(6)

i=−K j=−K

For static regions (t → ∞) and constant β < 1, the optimum parameters are given by a1 α= (1 − β) (7) 2a2 and γ = 1 − (α + β). Table 1 shows parameters a1 and a2 for typical ﬁlters used.

Results

Diverse clinical sequences as well as a X-ray image sequence showing a dynamic test object were used in the evaluation. For quantiﬁcation in terms of signal-tonoise-ratio (SNR), only the ROX was evaluated as the clinical irrelevant shuttered regions falsiﬁed the results by approximately 5 to 10 dB. State-of-the-art ﬁltering without noise estimation or motion detection generated perceivable motion artifacts even for low temporal ﬁlter portions of about 5% (β = 0.05). At the same time, static image regions were hardly denoised with these parameter settings. Contrarily, even for strong noise reduction settings no temporal or spatial artifacts were observed when applying the presented

M. Hensel, T. Pralow and R.-R. Grigat

Fig. 2. Original image detail (left), motion artifacts due to temporal averaging at dark structures (middle), processed with proposed method (right)

method. Due to strong temporal ﬁltering of broadly static image regions, clearly improved image quality and good temporal stability was notable in large image regions (Fig. 2). Objective evaluation of frames exhibiting about 18.5 dB SNR resulted in SNR improvement (SNRI) of about 9.1 dB.

Discussion

The presented method exhibits eﬃcient noise reduction while, at the same time, preventing artifacts. Motion artifacts, typical for temporal averaging, are prevented by motion detection based on processing positive and negative instead of absolute gray value diﬀerences. Artifacts due to blockwise binarization were eliminated by integration of a motion probability based on signal-dependent noise estimation. By this, improved adaptation of the method to signal-dependent noise has been achieved and the former discrete ﬁlter parametrization could be replaced by a continuous parametrization.

References 1. Hensel M, Wiesner G, Kuhrmann B, et al. Motion Detection for Adaptive SpatioTemporal Filtering of Medical X-Ray Image Sequences. In: Procs BVM; 2005. p. 68–72. 2. Lagendijk RL, van Roosmalen PMB, Biemond J, et al. Video Enhancement and Restoration. In: Bovik A, editor. Handbook of Image and Video Processing. 2nd ed. Elsevier Academic Press; 2005. p. 275–294. 3. Sanchez-Marin FJ, Srinivas Y, Jabri KN, Wilson DL. Quantitative Image Quality Analysis of a Nonlinear Spatio-Temporal Filter. IEEE Trans Image Proc 2001;10(2):288–295. 4. Aufrichtig R, Wilson DL. X-Ray Fluoroscopy Spatio-Temporal Filtering with Object Detection. IEEE Trans Med Imaging 1995;14(4):733–746. 5. Konrad J. Motion Detection and Estimation. In: Bovik A, editor. Handbook of Image and Video Processing. 2nd ed. Elsevier Academic Press; 2005. p. 253–274. 6. Hensel M, Lundt B, Pralow T, Grigat RR. Robust and Fast Estimation of SignalDependent Noise in Medical X-Ray Image Sequences. In: Procs BVM; 2006.

Automatic In-Vitro Oriﬁce Area Determination and Fluttering Analysis for Tricuspid Heart Valves Tobias Hahn1 , Alexandru Paul Condurache1 , Til Aach2 , Michael Scharfschwerdt3 and Martin Misfeld3 1 Institute for Signal Processing, University of Luebeck Institute of Imaging and Computer Vision, RWTH Aachen University 3 Department of Cardiac Surgery, University Clinic of Schleswig-Holstein, Luebeck Email: [emailprotected] 2

Abstract. Patients suﬀering from a tricuspid heart valve deﬁciency are often treated by replacing the valve with artiﬁcial or biological implants. In case of biological implants, the use of porcine heart valves is common. To supply patients with the best such xenograft implants, quality assessment and inspection methods for biological heart valves are necessary. We describe initial approaches to automatically testing tricuspid heart valves in-vitro in an artiﬁcial circulation system. We use two criteria to judge the quality of the valve, i.e. the evolution of the oriﬁce area during a heart cycle and the ﬂuttering of the valve leaﬂets in the stream.

Introduction

Replacing a deﬁcient tricuspid heart valve with a xenograft implant is a common procedure in cardiac surgery [1]. Porcine heart valves (see Figure 1) are favoured for this procedure because of their similarity to the human heart valve and the easy access to these valves due to food production. With the high number of available porcine valves quality assessment becomes a vital task to provide the surgeons with the best material only. Major quality criteria of these valves are the oriﬁce area, the opening and closing behaviour and the ﬂuttering of the leaﬂets in the blood stream. The determination of the oriﬁce area currently is done - if at all - manually by the researcher. An assessment of the ﬂuttering behaviour is provided only by visual observation. In our work we test and compare diﬀerent algorithms that automatically determine the oriﬁce area. We have implemented three threshold based methods for this task. We ﬁrst segment the oriﬁce and then compute its area by counting all object pixels. Using the segmented oriﬁce we perform a quantiﬁcation of the leaﬂet ﬂuttering. Section 2 of this paper explains the test apparatus and the state of the art. In section 3 we describe our algorithms. Results are shown in section 4, and section 5 closes with conclusions.

The Test Apparatus

Figure 2 a) depicts the test environment as described in [2]. From a reservoir (1) a transparent ﬂuid is transported through a disc valve(4) by a piston pump

T. Hahn et al.

Fig. 1. Frames 47 and 124 from heart valve sequences A and B respectively showing boundary outlines, projection lines and anchor points

(2), which is driven by a waveform adapted cam plate (3). After passing an input compliance (5) the ﬂuid is pressed through the inspected heart valve (11) into a visualization chamber (7) located in another ﬂuid reservoir (6). Pressure sensors (10) are installed below and above the heart valve. Passing an aortic compliance (9) the ﬂuid reaches a height variable column and ﬂows back to the ﬁrst reservoir. The heart valve is illuminated by light sources (12) outside the ﬂuid tank (6) made from perspex. A high-speed video camera (13) takes images of the heart valve. The digitized images are stored on a PC (14). Images of the heart valve are taken with 500 fps and an interlaced resolution of 480x420 gray pixels. An image sequence of one valve cycle has a duration of half a second leading to 250 frames. Figure 1 shows two heart valve images. Determing the oriﬁce area is done manually for each frame so far. The researcher has to select a set of points at the leaﬂet’s tips and the leaﬂet’s connections to deﬁne six triangles. Their summed areas give an approximation to the oriﬁce area. An evaluation of the leaﬂet ﬂuttering is currently done only qualitively by visual inspection. The oriﬁce approximation neglects the rather curved boundary of the leaﬂets. The methods we propose here are intended to give the researcher a more precise and reproducible way to calculate the area. The ﬂuttering analysis of the leaﬂets oﬀers a new way to quality inspection and is more reproducible than just human inspection.

Methods and Algorithms

Algorithm Overview In Figure 2 b) an overview of our algorithm is shown. In the preprocessing part we perform Tophat ﬁltering, equalize illumination variations and by user interaction obtain once for the whole sequence three anchor points situated at the intersections between the leaﬂets and the valve’s ﬁxation ring (see Figure 1). The oriﬁce area is determined by one of three threshold methods. These methods include a user-deﬁned threshold, a threshold automatically selected with the method proposed by Otsu in [3] and a threshold derived from the Finite mixture models algorithm proposed by Figueiredo and Jain in [4]. For each image we obtain the oriﬁce area by binarizing it with the threshold. To quantify the ﬂuttering of the heart valve leaﬂets we ﬁrst obtain the oriﬁce area

Automatic In-Vitro Oriﬁce Area Determination and Fluttering Analysis

Fig. 2. Scheme of test apparatus (a) and overiew of algorithm (b)

boundary by subtracting the eroded area segmentation from the result itself. As the tricuspid heart valve consists of three leaﬂets we split up the boundary at those contour points closest to the anchor points. This results in three segments each of which belonging to one leaﬂet. For each segment we perform a Fourier analysis of the contour points. The contour points coordinates are either treated as complex numbers (2-D data) or are projected to the connection lines between the splitting points leading to 1-D data and also oﬀering better control over the sampling rate. In Figure 1 the oriﬁce boundary outlines, the splitting points and the projection lines are depicted for both valves. Oriﬁce Area Determination Using Thresholds User-deﬁned threshold: The threshold used for binarization can be directly speciﬁed by the user as a certain percentile of the histogram of the heart valve image. Although we were able to obtain good results, this method has the downside that one has to pick the right threshold by trial and error. Otsu threshold: The threshold can be automatically selected using the method from Otsu [3]. It is histogram-based and assumes the histogram to be bimodal. It selects the threshold by maximizing the separability between the two resulting classes, i.e. dark ’oriﬁce area’ and bright ’leaﬂet area’. Finite mixture models threshold: In our case the histogram turned out not to be bimodal, thus a method is needed that can separate multiple classes, i.e. generate multiple thresholds. Knowing that the oriﬁce area is the darkest structure in our images we can select the threshold that separates the histogram peak with the smallest mean from the rest. We use the Finite Mixture Models [4] for this purpose. For a histogram with n samples consisting of multiple distributions the FMM algorithm returns estimates for the number k of distributions in the histogram, the mean values µi and the variances σi2 . Initialized with a random set of m Gaussian distributions having the probabilities αi the algorithm ﬁts the estimates iteratively to the histogram data y using an expectation maximization method and a quality criterion for de-

T. Hahn et al.

Fig. 3. Oriﬁce areas for ﬁve sequences (a) and ﬂuttering energies for sequences A (solid line) and B (dashed line) (b)

terming the number of distributions. For each distribution pi and sample yj we calculate the probability wij that this sample was generated by this distri−1(

yj −µi 2

)

. In each iteration step t the mean valbution as pi (yj ) = αi σ √1 2π e 2 σi i ues and variances are updated for each distribution according to µi (t + 1) = n −1 n j j −1 n 2 2 j ( nj=1 wij ) j=1 yj wi and σi (t + 1) = ( j=1 wi ) j=1 (yj − µi (t + 1)) wi . Furthermore the distribution probabilities α are evaluated and eventually set to 0 if the current distribution does not ﬁt onto the data. To determine the number of distributions hidden in the histogram in each iteration step the criterion knz (N +1) knz n i − log(p(Y |µ, σ 2 )) is L(µ, σ 2 , α, Y ) = N2 i:αi >0 (log( nα 12 )) + 2 log 12 + 2 calculated. With Y being the complete set of samples, N being the number of parameters deﬁning the distributions (2 for 1-dimensional Gauss distributions) and knz being the number of non-zero (i.e. α > 0) distributions this criterion reaches its minimum at that iteration where the number k of distributions ﬁts optimal to the histogram data. Assuming that our histogram consists of at least two classes and that the darkest class, i.e. the distribution with the lowest mean, represents the oriﬁce, we obtain the threshold by calculating the intersection point of the two distributions with the lowest mean values. Fluttering Analysis We apply a Fourier transform to the oriﬁce boundary obtained from thresholding to analyze the ﬂuttering behaviour of the heart valve leaﬂets. High-pass ﬁltration suppresses the frequencies corresponding to the leaﬂets’ curved shape showing only frequencies characteristic for the leaﬂet ﬂuttering. Then we compute the energy in the spectrum thus obtained. Only images where the oriﬁce area is above a certain minimum value are considered.

Results

Figure 3 a) shows the area over time plot for ﬁve of our test sequences. The plot exhibits the typical plateau-like form for heart valve oriﬁce areas. Beginning in the closed state the heart valves open rapidly, remain open for a while and then

Automatic In-Vitro Oriﬁce Area Determination and Fluttering Analysis

Fig. 4. Projected boudary (a) and magnitude spectra (b) for frame 124 of sequence B

close again. The plateau is slightly falling to the right side. It is also exhibits some variations due to the heart valves’ leaﬂets ﬂuttering in the blood stream thus altering the oriﬁce area to some extent. Figure 4 a) shows the contour coordinates projected to 1D and freed from the mean for each of the three heart valve leaﬂets. This data was taken from our test sequence B, frame 124. Figure 4 b) shows the corresponding energy of the data. We used high-pass ﬁltration to cut oﬀ the frequencies belonging to the leaﬂets’ natural curvature. Figure 3 b) shows the ﬂuttering energies for sequences A (solid) and B (dashed). The strong ﬂuttering of B is clearly visible in the plot.

Conclusions and Outlook

We tested three threshold based approaches to determine the oriﬁce area of heart valves automatically. Although these were only ﬁrst tests the methods showed good results. We were able to derive typical oriﬁce plots from the given data. This provides the researcher with a tool to perform the so far tedious work automatically. The ﬂuttering analysis we introduced here is a novel approach giving information about the behaviour of the leaﬂets ﬂoating in the blood stream. Further work on this subject may include snakes as an alternative approach to determine the oriﬁce area. The selection of the anchor points may be automated using a Hough transform.

References 1. Vongpatanasin W, Hillis LD, Lange RA. Prosthetic Heart Valves. MP 1996;355(6):407–416. 2. Scharfschwerdt M, Misfeld M, Sievers HH. The Inﬂuence of a Nonlinear Resistance Element upon In-Vitro Aortic Pressure Tracings and Aortic Valve Motions. ASAIO 2004;. 3. Otsu N. A Threshold Selection Method from Gray-Level Histograms. IEEE TSMC 1979;9(1):62–66. 4. Figueiredo MAT, Jain AK. Unsupervised Learning of Finite Mixture Models. IEEE TPAMI 2002;24(3):381–396.

Atlas-basierte 3D-Rekonstruktion des Beckens aus 2D-Projektionsbildern Hans Lamecker1, Thomas H. Wenckebach1 , Hans-Christian Hege1 , Georg N. Duda2 und Markus O. Heller2 1

Visualisierung und Datenanalyse, Zuse Institut Berlin, 14195 Berlin Centrum f¨ ur Muskuloskeletale Chirurgie, Charit´e-Berlin, 13353 Berlin Email: [emailprotected]

Zusammenfassung. H¨ auﬁg liegt einer computergest¨ utzten Operationsplanung lediglich ein Satz von wenigen 2D-R¨ ontgenbildern zugrunde. Dennoch ist es ein Anliegen, auf Basis solcher Daten R¨ uckschl¨ usse auf die dreidimensionale Anatomie des Patienten zu ziehen. In dieser Arbeit wird ein Verfahren vorgestellt, das mithilfe eines statistischen 3DFormmodells (SFM) des Beckens die geometrisch sowie toplogisch komplexe 3D-Form aus wenigen R¨ ontgenbildern rekonstruiert. Dies geschieht durch eine Optimierung, welche den Abstand der Silhouette des Modells in den Projektionsebenen zur Silhouette des Objekts in den R¨ ontgenbildern minimiert. Das Verfahren wird an 23 synthetisch erzeugten Datens¨ atzen validiert.

Einleitung

Der k¨ unstliche Gelenkersatz der H¨ ufte ist ein Standardverfahren zur Behandlung degenerativer Gelenkerkrankungen. Allein in Deutschland werden j¨ahrlich ca. 120.000 H¨ uftendoprothesen implantiert. Aufgrund der Alterung der Bev¨olkerung ist in Zukunft mit einem starken Anstieg der Zahl der Implantationen zu rechnen. Obwohl bekannt ist, dass die Belastungen des H¨ uftgelenkes eine wesentliche Rolle f¨ ur die langfristige Funktion und den Erfolg des k¨ unstlichen Gelenkersatzes spielen, stehen dem Chirurgen bei der Planung des Eingriﬀes derzeit keine validen Informationen zur erwarteten Gelenkbelastung zur Verf¨ ugung. Die computergest¨ utzte Planung soll daher helfen, den Erfolg der Behandlung zu verbessern. Dazu sollen die auftretenden Kr¨ afte und Belastungen des Gelenkes vor und nach dem Eingriﬀ mit Hilfe validierter biomechanischer Modelle berechnet werden, um eine m¨ oglichst optimale biomechanische Rekonstruktion des H¨ uftgelenkes f¨ ur jeden Patienten gew¨ ahrleisten zu k¨ onnen. Die Grundlage f¨ ur die Operationsplanung sind in der Regel R¨ontgenbilder. W¨ahrend aus einer Becken¨ ubersichtsaufnahme (koronale Projektion) bereits wesentliche Daten f¨ ur eine Anpassung biomechanischer Modelldaten zur Berechnung der Belastungsbedingungen der H¨ ufte gewonnen werden k¨onnen [1], erfordert eine genauere Belastungsanalyse die Kenntnis der 3D-Geometrie der Anatomie, insbesondere der Knochen und Muskeln. Die Herausforderung besteht somit darin, aus R¨ ontgenbildern eine 3D-Rekonstruktion der anatomischen Objekte als Grundlage f¨ ur die nachfolgende Berechnung der Belastungen zu bestimmen.

Atlas-basierte 3D-Rekonstruktion des Beckens aus 2D-Projektionsbildern

Abb. 1. Registrierung eines 3D-Formmodells mit zwei R¨ ontgenbildern.

Ein Großteil der Literatur zum Thema 2D/3D-Registrierung (siehe [2] und enthaltene Referenzen) geht davon aus, dass ein pr¨a-operativer 3D-Datensatz des Patienten vorliegt, der dann mit intra- oder post-operativ gewonnene Daten registriert werden muss. In unserem Fall ist jedoch kein 3D-Referenzdatensatz vorhanden. Der Einsatz eines SFM erscheint daher sinnvoll f¨ ur die Rekonstruktion einer unbekannten 3D-Anatomie aus 2D-R¨ ontgenbildern. Verschiedene Arbeiten verfolgen einen ¨ ahnlichen Ansatz: Fleute et al. [3] verwenden ein SFM des distal gelegenen Teils des Femur zur Registrierung mit R¨ ontgenbildern aus einem C-Arm. Dazu werden die Abst¨ande zwischen den Konturen der SFM-Oberﬂ¨ ache und einer diskreten Anzahl von Projektionsstrahlen der R¨ ontgenbilder mittels eines ICP-Algorithmus minimiert. Benameur et al. [4] benutzen ein SFM der Wirbelknochen zur Registrierung und Segmentierung von R¨ ontgenbildern. Die Registrierung besteht in der Minimierung eines Kantenpotentials, das den Abstand der proji*zierten Konturen von den Konturen der R¨ ontgenbilder misst. Yang et al. [5] bauen auf ein hybrides Formmodell zur Rekonstruktion von Femuren aus R¨ ontgenbildern. Als Abstandsmaß dient die Korrelation zwischen simulierten Dickenbildern des Formmodells und den R¨ ontgenbildern.

Methoden

Das verwendete SFM wird aus einem Satz von individuellen Trainingsdaten erzeugt. Die Herausforderung dabei liegt in der Identiﬁzierung von korrespondierenden Punkten auf den Trainingsoberﬂ¨ achen. Hierzu wurde das von Lamecker et al. [6, 7] beschriebene Verfahren verwendet. Als Ergebnis lassen sich alle Trainingsformen in einem gemeinsamen Vektorraum darstellen. Eine Hauptkomponentenzerlegung (PCA) dieses Ensembles liefert ein lineares Modell f¨ ur die zuge lassene Formvariabilit¨ at: S(b, T ) = T (v + nk=1 bk pk ), wobei v das Formmittel, pk die Eigenmoden der Formvariabilit¨ at, bk die Formgewichte und T eine aﬃne Transformation darstellen (vgl. Abb. 2, links).

H. Lamecker et al.

Abb. 2. Links: Statistisches Formmodell des Beckens aus einer Trainingsmenge, Mitte: Kantenerkennung in R¨ ontgenbildern, Rechts: Silhouette des Beckens in R¨ ontgenbildern.

F¨ ur eine gegebene Kamerakalibrierung K (Lage und Ort der R¨ontgenquelle relativ zu den Projektionsebenen) und eine gegebene Instanz des Formmodells S(b, T ) wird die Silhouette des Modells in den Projektionsebenen berechnet und als 2D-Bild gerastert. Die Silhouette des Objektes in den R¨ontgenbildern wird einmalig interaktiv bestimmt: Nach einer automatischen Kantenbestimmung (Canny-Filter) werden unbrauchbare Konturen entfernt bzw. zur Silhouette vervollst¨ andigt (vgl. Abb. 2, mitte und rechts). Gegeben zwei Silhouetten s und s , sei d(x, s ) = minx ∈s x − x der Abstand zwischen einem Punkt x ∈ s und der Silhouette s . Das zu minimierende ¨ Ahnlichkeitsmaß ist der symmetrische quadratische Abstand

d = x∈s d(x, s )2 dx + x∈s d(x, s)2 dx. Zur eﬃzienten Bestimmung des Abstandes zwischen den Silhouetten werden sowohl f¨ ur die Modell- als auch die R¨ ontgensilhouetten Distanzfelder berechnet. Die Berechnung von d besteht dann in der Abtastung und Summierung der Distanzen der Punkte einer Silhouette zur jeweils anderen Silhouette. Das Distanzfeld der Modellsilhouette muss dabei st¨ andig aktualisiert werden. Die gesuchte Rekonstruktion ist die L¨ osung des Optimierungsproblems (T ∗ , b∗ , K ∗ ) = arg min(T,b,K) d. Zur Minimierung wird ein Suchverfahren verwendet, das einem Gradientenabstieg mit Liniensuche ¨ahnelt. Die Suchrichtung wird nicht mittels einer numerischen Approximation der partiellen Ableitungen ¨ berechnet, sondern mit einem Verfahren, welches das Ahnlichkeitsmaß an den Grenzen eines sehr großen Intervalls auswertet. Iterativ werden die Intervalle verkleinert, wodurch eine quasi-globale Optimierung erzielt wird. Zur Vermeidung von lokalen Minima und zur Steigerung der Eﬃzienz wird die Registrierung in einer Datenpyramide durchgef¨ uhrt, d.h. die Silhouetten werden in verschiedenen Auﬂ¨ osungen in der Optimierung ber¨ ucksichtigt [8].

Ergebnisse

Zur Validierung wurden 23 CT-Datens¨ atze des Abdomens ohne Knochendefekte verwendet (Schichtabstand 5mm). Zu allen Datens¨atzen lagen manuelle Segmentierungen des Beckens als Goldstandard zur Auswertung vor. Das SFM des Beckenknochens wurde ebenfalls aus diesen Daten erzeugt [6]. Aus den Datens¨ atzen wurden mittels eines Volume Renderings k¨ unstliche R¨ontgenbilder be-

Atlas-basierte 3D-Rekonstruktion des Beckens aus 2D-Projektionsbildern

Tabelle 1. Quantitative Ergebnisse der verschiedenen Experimente. Experiment Leave-All-In (KO) Leave-All-In (KO-SA) Leave-One-Out (KO) Leave-One-Out (KO-SA) Leave-One-Out (DFM)

Mittelwert [mm] 1, 5 ± 0, 5 1, 3 ± 0, 5 2, 6 ± 0, 4 2, 4 ± 0, 4 2, 0 ± 0, 2

Median [mm] 1, 2 ± 0, 3 1, 1 ± 0, 5 2, 1 ± 0, 3 2, 0 ± 0, 3 1, 6 ± 0, 2

Maximum [mm] 8, 8 ± 3, 3 7, 9 ± 2, 6 17, 6 ± 6, 5 14, 9 ± 3, 1 13, 3 ± 2, 6

rechnet. Diese sind hier allerdings lediglich f¨ ur die Deﬁnition geeigneter Projektionsebenen erforderlich. Die Silhouetten dieser R¨ontgenbilder wurden direkt aus den manuellen Segmentierungen berechnet. Ziel der Validierung ist, die erreich¨ bare Genauigkeit des Verfahrens mit dem vorgeschlagenen Ahnlichkeitsmaß zu untersuchen. F¨ ur die konventionelle Operationsplanung werden bisher außerdem bereits kalibrierte R¨ ontgenbilder unter standardisierten Aufnahmebedingungen erstellt [9]. Daher wurden in der Optimierung hier nur die Formgewichte b optimiert (T, K sind bekannt). Bei Unkenntnis von T, K ist mit einer Verschlechterung der Ergebnisse zu rechnen. Die Experimente wurden f¨ ur den Fall einer R¨ ontgenaufnahme (koronal, KO) und zweier Aufnahmen (koronal und sagittal, KO-SA) durchgef¨ uhrt (vgl. Abb. 1). Als Fehler wurde der mittlere symmetrische Fl¨ achenabstand zwischen Optimierungsergebnis und Goldstandard gemessen. Den geringsten Wert f¨ ur den Fehler erreicht man durch direktes Fl¨ achenmatching (DFM) des SFM an den Goldstandard in 3D (Minimierung des mittleren quadratischen Fehlers, vgl. [6]). Als Test der Optimierungsstrategie wurde das Verfahren mit einem SFM durchgef¨ uhrt, das die Formen der zu rekonstruierenden R¨ontgenbilder enth¨alt (Leave¨ All-In-Test). Hier wird optimale Ubereinstimmung erwartet. Im Leave-One-OutTest ist der zu rekonstruierende Datensatz nicht im SFM enthalten. Dies stellt den realistischen Fall dar. Die Ergebnisse sind in Tabelle 1 zusammengefasst.

Diskussion

Ziel der Arbeit ist die Validierung eines neues Verfahrens zur 3D-Rekonstruktion des Beckens aus 2D-Projektionsbildern anhand synthetischer Datens¨atze. Die Mittelwerte der gemessenen Fehler liegen mit einer Gr¨oßenordnung von wenigen Millimetern in einem Bereich der erwarten l¨ asst, dass evtl. bereits aus einer einzigen Projektion eine f¨ ur die Berechnung der Belastungsbedingungen hinreichende Absch¨ atzung der 3D-Geometrie m¨ oglich ist [1]. Der Leave-One-Out-Test zeigt weiter, dass der gr¨ oßte Anteil des Fehlers in der Unvollst¨andigkeit des SFM liegt. Dieses gilt es daher in Zukunft um weitere Trainingsdaten zu erweitern. F¨ ur die klinische Anwendung bleibt zu kl¨ aren, wie das Verfahren durch degenerative Ver¨ anderungen der H¨ uftknochen beeinﬂusst wird. Hierzu muss in Zukunft untersucht werden, inwiefern die Methode bei unvollst¨andigen Silhouetten eine ausreichende Genauigkeit liefert. Der Einsatz eines Formmodells normaler

H. Lamecker et al.

Becken ist jedoch auch f¨ ur eine solche Situation sinnvoll, da das Ziel der Rekonstruktion vor allem diejenigen Bereiche sind, in denen wichtige Muskeln ansetzen. Solche Bereiche sind h¨ auﬁg nicht degenerativ ver¨andert. Eine m¨ ogliche Steigerung der Genauigkeit besteht darin, mehr Information aus ¨ den Bilddaten und dem Modell in das Ahnlichkeitsmaß zu integrieren. Versuche mit weiteren aus der Kantendetektion berechneten Konturen zu arbeiten, wie Benameur et al. [4] es vorschlagen, schlugen fehl, da die Konturberechnung im Allgemeinen nicht robust genug ist. R¨ ontgenbilder messen die Dichte des Knochens entlang eines Projektionsstrahls. Allerdings brachten Versuche, in denen ausschließlich Dickenbilder des SFM mit den R¨ontgenbildern registriert wurden (¨ahnlich [5]), unbefriedigende Ergebnisse, da der Inhom*ogenit¨at des Knochens im Dickenbild nicht Rechnung getragen wird. Zuk¨ unftig sollen daher Silhouettenund Dickeninformation im Registrierungsprozess kombiniert werden. Zusammenfassend muss man davon ausgehen, dass mit dem hier vorgeschlagenen Ansatz eine f¨ ur die Absch¨ atzung der Belastungen des Gelenkes hinreichend genaue 3D-Rekonstruktion des Beckens erzielt werden kann. Dies soll auch an klinischen Daten evaluiert werden. F¨ ur die vollst¨andige Analyse der Belastungen im Rahmen des H¨ uftgelenkersatzes wird das Verfahren in einem n¨achsten Schritt auf das Femur erweitert werden. Danksagung. Hans Lamecker wird gef¨ ordert durch das DFG Forschungszentrum Matheon “Mathematik f¨ ur Schl¨ usseltechnologien” in Berlin. Die Autoren danken Olaf Etzmuss f¨ ur erste Untersuchungen und Implementierungen.

Literaturverzeichnis 1. Heller MO, Bergmann G, Deuretzbacher G, et al. Musculo-skeletal loading conditions at the hip during walking and stair climbing. J Biomech 2001;34(7):883–893. 2. van de Kraats EB, Penney GP, et al. Standardized Evaluation Methodology for 2D-3D Registration. IEEE Trans Med Imaging 2005;24(9):1177–1190. 3. Fleute Markus, Lavall´ee St´ephane. Nonrigid 3-D/2-D Registration of Images Using Statistical Models. In: Procs MICCAI. vol. 1679; 1999. p. 138–147. 4. Benameur S, Mignotte M, Parent S, et al. 3D/2D Registration and segmentation of scoliotic vertebrae using statistical models. Computerized Medical Imaging and Graphics 2003;27(5):321–337. 5. Tang ThomasSY, Ellis RandyE. 2D/3D Deformable Registration Using a Hybrid Atlas. In: Procs MICCAI. vol. 3750; 2005. p. 223–230. 6. Lamecker Hans, Seebaß Martin, Hege HansChristian, Deuﬂhard Peter. A 3D Statistical Shape Model of the Pelvic Bone For Segmentation. In: Procs SPIE Medical Imaging: Image Processing. vol. 5370; 2004. p. 1341–1351. 7. Lamecker Hans, Lange Thomas, Seebaß Martin. Erzeugung statistischer 3DFormmodelle zur Segmentierung medizinischer Bilddaten. In: Procs BVM; 2003. p. 398–403. 8. Wenckebach ThomasH. Volumetrische Registrierung zur medizinischen Bildanalyse. Master’s thesis. Institut f¨ ur Informatik, Humboldt-Universit¨ at zu Berlin; 2005. 9. The B, Diercks RL, Stewart RE, et al. Digital correction of magniﬁcation in pelvic x-rays for preoperative planning of hip joint replacements: theoretical development and clinical results of a new protocol. Med Phys 2005;32(8):2580–2589.

Machine-Based Rejection of Low-Quality Spectra and Estimation of Brain Tumor Probabilities from Magnetic Resonance Spectroscopic Images Bjoern H. Menze1 , B. Michael Kelm1 , Daniel Heck1 , Matthias P. Lichy2,3 and Fred A. Hamprecht1 1

Interdisziplin¨ ares Zentrum f¨ ur Wissenschaftliches Rechnen (IWR), Universit¨ at Heidelberg, 69120 Heidelberg 2 Abteilung Radiologie, Deutsches Krebsforschungszentrum (dkfz), 69120 Heidelberg 3 currently: Diagnostische Radiologie, Universit¨ at T¨ ubingen, 72076 T¨ ubingen

Abstract. Magnetic resonance spectroscopic images (MRSI) carry spatially resolved information about the in vivo metabolism, however, their evaluation is diﬃcult. Problems arise especially from artifacts and noise, yielding non-evaluable signals in many voxels. We propose a two-step approach to the processing of MRSI. In the ﬁrst step a non-linear classiﬁer is employed in every voxel to determine whether the spectral signal is evaluable, and if so, the tumor probability is computed in the second step. Thus, the quality control is strictly separated from the diagnostic evaluation of the spectrum. For an assessment of the proposed approach we consider MRSI-based brain tumor detection and localization and a tumor probability mapping by pattern recognition. In a quantitative comparison against the standard operator-controlled processing our interaction-free approach shows similar to superior performance.

Introduction

Magnetic resonance spectroscopy (MRS) allows the study of physiological and patho-physiological processes in the living tissue. Thus, characteristic changes in the relative concentration of a number of metabolites can be detected by MRS. Magnetic resonance spectroscopic imaging (MRSI) also allows for a spatial mapping of their occurrence. – One of the major applications of MRSI in clinical routine is in the diagnosis of tumor, in particular in the brain, but also in the breast or prostate. MRSI in Tumor Diagnosis The standard approach in the analysis of MR spectra is a ﬁtting of resonanceline shaped template functions to the spectral pattern. The parameters of these (amplitude, line width, area) provide the basis for the following diagnostic interpretation of the spectrum, e.g. using the value of the resonance line integral

B.H. Menze et al.

of Choline relative to the integral of Creatine as an indicator for the malignity of a tissue change. Unfortunately, the robustness of this resonance line “quantitation” is limited in the presence of noise or artifacts within the spectra, and a visual control by the operator is always required. This eﬀort hampers the routine use of MRSI, but is compulsory to assess the diagnostic value of an MR spectrum. An alternative approach, which is primarily suited for highly resolved (short echo time) spectra, bypasses the quantitation and uses algorithms from pattern recognition and machine learning. Tissue type and state of disease are directly inferred on the spectral pattern, without a model based reduction of the spectral information. Pattern recognition on the full spectrum, e.g. by a projection of the spectral image to a “nosologic” or diagnostic map, is more robust than the ﬁtting of adaptive resonance line functions [1]. Nevertheless, this method also fails in the presence of spectral artifacts or at low signal-to-noise ratio. So, an inspection of the spectra is also required here and is limiting a fully automated analysis also for this approach. Published examples for the use of MRSI in tumor diagnostic can be found for both analysis approaches, e.g. in the non-linear classiﬁcation of brain tumor [2] (EU-project INTERPRET), or in the quantitation of MR spectroscopic images in the detection of recurrent tumors after radiation therapy [3]. The most important obstacle for an application of these automated routines is the lack of a speciﬁc assessment about the conﬁdence in the diagnostic analysis of a spectrum. While a number of algorithms are available for the diagnostic interpretation of MR spectra [2, 3, 4], quality control and the evaluation of conﬁdences still demand a large amount of operator time and bind the resources of the human expert. Two Step Analysis Using the example of the detection and localization of brain tumor, we propose in the following a two-step procedure which is able to overcome some of the mentioned limitations in the automated analysis of MR spectra. In the ﬁrst step, we demonstrate the applicability of pattern recognition in the classiﬁcation of evaluable and non-evaluable (“nice vs. noise”) spectra. In the following diagnostic analysis, we use regression methods from chemometrics for pattern recognition and the generation of probabilistic tumor-maps with conﬁdence assessment. Overall, we present an automated processing of MR spectroscopic images which yields a data product similar to any other standard imaging modality. Most important in comparison to the standard processing of MRSI, it is not a collection of localized spectra that need to be interpreted laboriously, but the mapping of the tumor-related information within the spectral image. Results of this approach will be compared against the outcome of a operatorcontrolled analysis. In order to allow a direct application, the algorithms can be integrated into the radiation therapy planning software at the German Cancer Research Center (dkfz) [5].

Machine-Based Rejection of Low-Quality Spectra

Methods

Data The data was acquired on a standard MR scanner (Siemens, Magnetom Vision, 1.5T, 1 H-MRS, TE 135 ms) at the dkfz. Spectral images of 31 patients with tumor diagnosis (astrocytoma, meningioma, glioma, metastases – conﬁrmed during follow-up) were available at a spatial resolution of 322 and under a spatial resolution of approx. 1 cm voxel side-length (ﬁg. 1a). Within this data 269 spectra could be labeled as “tumor”, “tumor border” or “normal state” (69/49/151, ﬁg. 1d). – Spectral images of another 36 patients were classiﬁed completely by visual inspection according to the signal quality and resulted in approx. 37000 spectra (10% “nice”) for the training of the “nice vs. noise” classiﬁer. Classiﬁer – Quality Assessment While the signal-to-noise ratio can be assessed easily by linear methods, the separation of the various artifacts demands for the training of a non-linear classiﬁer (“randomForest” decision tree ensemble [1]). RandomForests extend the idea of bagging tree ensembles by seeking optimal splits for each tree node on random subspaces of the data. Outcome of a randomForest classiﬁcation is either the binary class membership, or the normalized number of votes from the tree ensemble, which allows for a subjective readjustment of the rejection threshold.– Diﬀerent feature representations of the data (diﬀerent normalization strategies, data binning) were tested and optimized according to the classiﬁcation error. As the “noise” class of the training data only contained the most frequent artifacts, it was extended by uniformly distributed random samples of the spanned feature space of all spectra. This nonparametric outlier detection ensures a proper signal classiﬁcation even in the presence of previously unknown artifacts. Classiﬁer – Tumor Detection The character of the binary “normal state vs. tumor” classiﬁcation problem allows the application of linear models on the data. In order to overcome the sparseness of the training problem (256 spectral channels (P ) in the spectral region of Choline, Creatine and NAA, 269 (N ) highly collinear spectra), the spectral regions relevant for the classiﬁcation were identiﬁed in a ﬁrst feature selection by help of the randomForest Gini importance measure (ﬁg. 1d), the accumulated Gini decrease for all nodes of all trees within the ensemble. A partial least squares regression (PLS) [1] on these variables was used for the discrimination of the two tissue types. To allow an assessment of the tumor probability, a logistic model was build on the score of the PLS subspace (ﬁg. 1c). The leave-‘one-patient’-out cross-validated PLS score was compared by means of the receiver-operator-characteristic (ROC) against the results of diﬀerent standard approaches [5], all based on resonance line quantitation (line integral ratios and linear combinations thereof, ﬁg. 1b).

B.H. Menze et al.

Results

Areas with high signal-to-noise ratio were reliably identiﬁed, virtually all artifacts being present in the data set were separated. An inspection of the leave-’one patient’-out shows that the cross-validated test error of 1.04 % is most probably due to ambivalent training labels. The classiﬁer proposed to reject 60 low quality spectra (16/10/34) from the data set for the tumor classiﬁcation. Within that data set, the randomForest Gini importance identiﬁed the spectral regions of the metabolites choline and NAA to be relevant for the classiﬁcation. Thus a subset of Psubset = 14 out of the 256 spectral channels was used for training and evaluation of the binomial PLS model (ﬁg. 1d). The resulting tumor classiﬁer reached an accuracy of 95.5% in the diﬀerentiation between normal state and both tumor groups, normal and tumor were separated without error (100% accuracy). These values, however, reduce to 92.5% and 93.3% when the 60 low quality spectra are included. Comparing the ROC of this approach with other methods as described in [5], we ﬁnd that the automated pattern recognition performs comparable to the outcome of a quantitation by a human operator (ﬁg. 1b with noise) or better (ﬁg. 1b without noise). All studied methods tend to show a higher reliability when low quality spectra are separated by the ﬁrst step classiﬁer.

Discussion

The proposed application of a speciﬁc, non-linear classiﬁer for the the separation of non-reliable spectra is independent from the diagnostic question and the examined organ. Although we observe optimal results in conjunction with our pattern recognition approach, such a quality assessment prior to the diagnostic analysis shows to be beneﬁcial also in conjunction with traditional resonance line quantitation. So far, only information theoretic criteria (e.g. Cramer-Rao-bounds) were used to assess the quality of a spectrum. Explicit and easily available knowledge about frequent artifacts which is the base of our noise classiﬁer are disregarded by these traditional quality measures. In the speciﬁc task of brain tumor detection and localization, we can show that the fully automated pattern recognition operates as reliably as the operatorsupervised resonance line quantitation. The high degree of automation of both steps, the trained quality control and the chemometrical tumor recognition, allows for the application of further image processing and the exploitation of spatial relations within the MR spectroscopic image. Promising extensions might be the automated segmentation of tumor regions or the combination of MRSI nosologic images with other diagnostic imaging modalities. The operator-free data processing is a prerequisite for the clinical application of spectral images with high spatial resolution and therefore thousands of spectra in each acquisition.

Machine-Based Rejection of Low-Quality Spectra

1.0 0.8 0.6 0.4

True positive rate

pred cho/cr cho/naa cho/cr + a*cho/naa

0.0

0.2

0.4

0.6

pred cho/cr cho/naa cho/cr + a*cho/naa

0.2

True positive rate

0.8

1.0

Fig. 1. a) Tumor probability map; left: PRESS box indicated, right: gold standard (T (t) - tumor (-border), N - normal), regions which were classiﬁed as low quality spectra were disregarded for display. b) Receiver-operator-characteristic of the presented PLS based pattern recognition and other methods as described in [5], including(left) and excluding (right) low quality spectra. c) Transformation of class densities and binomial tumor-probability in red-green color scale. (here: normal vs. both tumor groups) d): Average spectral patterns of healthy and tumorous tissue.

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.8

1.0

(b) tumor

healthy

0.6

normal

Cho

NAA Cr

3 2

tumor border tumor

0.5

2.2 ppm

randomForest Gini importance

Median Quartiles Outlier

healthy

10 12 14

3.3 ppm tumor tumor−border

tumor probability

class densities

(a) 4

0.4

False positive rate

100

PLS score

(c)

(d)

References 1. Menze BH, Lichy MP, Bachert P, et al. Optimal classiﬁcation of Lonag Echo Time In Vivo Magnetic Resonance Spectra in the Detection of Recurrent Brain Tumors. NMR in Biomedicine submitted;. 2. Devos A, Lukas L, Suykens JA, et al. Classiﬁcation of brain tumor using short echo time 1 H MR spectra. J Magn Reson 2004;p. 164–175. 3. Lichy MP, Pathow C, Schultz-Ertner D, et al. Follow-up of gliomas after radiotherapy: 1 H MR spectroscopic imaging for increasing diagnostic accuracy. Neuroradiology 2005;p. 1–12. 4. Kelm BM, Menze BH, Neﬀ T, et al. CLARET: a tool for fully automated evaluation of MRSI with pattern recognition methods. In: Procs BVM; 2006. . 5. Neﬀ T, Bachert P, Lichy MP, et al. Einsatz der Protonten-MR-spektroskopischen Bildgebung zur Zielvolumenbestimmung in der Strahlentherapieplanung bei glialen Hirntumoren. In: Procs BVM; 2004. p. 145–149.

Maximum-Likelihood-Ansatz zur Metallartefaktreduktion bei der Computertomographie May Oehler und Thorsten M. Buzug Department of Mathematics and Technology, RheinAhrCampus Remagen, 53424 Remagen, Deutschland Email:{moehler,buzug}@rheinahrcampus.de

Zusammenfassung. Aufgrund der Aufh¨ artung des R¨ ontgenspektrums treten besonders bei Metallen Artefakte in der CT-Bildrekonstruktion auf. Mathematisch liegt die Ursache in der Inkonsistenz der Daten des Radonraums. Zur Unterdr¨ uckung der Bildfehler werden die inkonsistenten Daten mit polynomialer Interpolation u uckt und anschließend ¨berbr¨ mit einem gewichteten Maximum-Likelihood-Expectation-MaximizationAlgorithmus rekonstruiert. Dieses Verfahren wird als λ-MLEM Algorithmus bezeichnet.

Einleitung

Die Computertomographie (CT) z¨ ahlt heute zu einem der wichtigsten bildgebenden Verfahren in der Medizin. Ein Problem bei der Beurteilung von CT-Bildern stellen Metallartefakte dar. Metalle in Form von Zahnf¨ ullungen, wie Amalgam oder Gold, sowie H¨ uftprothesen aus Stahl absorbieren auch die hochenergetische R¨ontgenstrahlung st¨ arker als das umliegende Gewebe und f¨ uhren zu falschen Schw¨ achungswerten. Im rekonstruierten Bild f¨ uhrt dies h¨auﬁg zu dunklen Streifen zwischen den Metallobjekten und zu nadelstrahlf¨ormigen, von den Objekten ausgehenden Streifen, die das umliegende Gewebe u ¨ berlagern und hierdurch die diagnostische Beurteilung erschweren. Zur¨ uckzuf¨ uhren sind diese Probleme auf ein niedriges SNR im Bereich des Metallschattens, ein erh¨ohtes Verh¨altnis von Streustrahlung zu Prim¨ arstrahlung und auf Strahlaufh¨artungsartefakte [1]. Um die diagnostische Beurteilung der Bilder, sowie die Planung von Operationen im Kieferbereich oder nach H¨ uftimplantatoperationen zu erleichtern, besteht ein Interesse an der Reduktion der Artefakte.

Material und Methoden

Als Standardrekonstruktion wird in der Computertomographie auf Grund ihrer Geschwindigkeit die geﬁlterte R¨ uckprojektion (FBP) eingesetzt. Dieses Verfahren stellt jedoch keine ad¨ aquate Methode im Umgang mit den inkonsistenten Daten im Radonraum dar. Einen besseren Ansatz hierf¨ ur stellt das 1984 eingef¨ uhrte

Maximum-Likelihood-Ansatz zur Metallartefaktreduktion

Maximum-Likelihood-Expectation-Maximization-Verfahren (MLEM) f¨ ur die Computertomographie dar [2]. Der Vorteil des MLEM-Verfahrens im Vergleich zur Rekonstruktion mit der geﬁlterten R¨ uckprojektion liegt darin, dass w¨ahrend der Rekonstruktion in geeigneter Weise Einﬂuss auf die aufgenommenen Rohdaten genommen werden kann. Weiterhin hat das MLEM-Verfahren das Potenzial die zur Aufnahme von CT-Bildern ben¨ otigte Dosis zu reduzieren, denn es ist in der Lage, Bilder aus Rohdaten zu rekonstruieren, die eine schlechte Statistik aufweisen und nur aus einer geringen Anzahl an Projektionen bestehen. Aus diesem Grund stellt es in der Nuklearmedizin das Standardverfahren zur Rekonstruktion dar. Zentrale Vorteile gegen¨ uber der FBP sind weiterhin das Fehlen der streifenf¨ ormigen Artefakte im Bild und ein exakteres betonen der Organ und K¨ orperkonturen [3]. Zun¨ achst werden die Rohdaten eines Torsophantoms der Firma CIRS Inc. mit dem Scanner Tomoscan M/EG der Firma Philips aufgenommen (Abb. 1a). Das Torsophantom wird mit zwei Stahlmarkern versehen. Die Vergleichsgrundlage bilden jeweils die rekonstruierten Bilder der Rohdaten dieses Torsophantoms ohne Stahlmarker in derselben Schnittlage. Den ersten Schritt zur Reduktion der entstehenden Metallartefakte stellt ein einfaches Reparaturverfahren des Rohdatenraumes dar (Abb. 1b). Hierzu werden zun¨achst die Metallmarker in einem vorl¨ auﬁgen, mit der FBP rekonstruierten Bild mit Hilfe eines Schwellwertes detektiert. Durch Vorw¨ artsprojektion der segmentierten Marker erh¨alt man ein reines Metallmarker-Sinogramm, das den Verlauf der Metallprojektionen im Rohdatenraum liefert. Alle Werte innerhalb dieses Sinogramms, die verschieden von Null sind, werden als inkonsistente Daten bezeichnet. In der Literatur [4], [5] und [6] ﬁnden sich verschiedene Ans¨ atze, die inkonsistenten Daten zu u ucken. In ¨ berbr¨ dieser Arbeit werden die inkonsistenten Daten mit polynomialer Interpolation innerhalb einer Projektion unter einem festen Winkel u uckt. ¨ berbr¨ Damit wird die unter dem Marker liegende Weichteilabsorption gesch¨atzt. Abbildung 1c zeigt das Ergebnis der polynomialen Interpolation (gepunktete schwarze Kurve) im Vergleich zur linearen Interpolation (schwarze StrichPunkt-Kurve), der Projektion, gemessen mit Metallmarkern (graue Kurve) und der Ground Truth, gemessen ohne Metallmarker (schwarze Kurve) anhand eines vergr¨ oßerten Proﬁl-Ausschnittes aus Abb. 1b unter einem Winkel von γ = 120◦ . Hierin zeigt sich, dass mit der polynomialen Interpolation eine gute N¨aherung an die Ground Truth gelingt. Anschließend werden die so gesch¨atzten Weichteildaten mit dem MLEM-Algorithmus rekonstruiert. Da es sich bei diesem Reparaturschema des Radonraums um einen rein pragmatischen Ansatz handelt, dem ein schwaches physikalisches Model zu Grunde liegt, kann an dieser Stelle nicht erwartet werden, dass die interpolierten Daten jenen entsprechen, die tats¨ achlich ohne Metallmarker gemessen werden. Aus diesem Grund werden im Folgenden die jeweiligen Projektionen, die durch einen Metallmarker verlaufen, sowie die entsprechenden Zeilen der dazugeh¨origen Systemmatrix. mit einem Vertrauensparameter λ gewichtet. Dieses neue Verfahren wird λ-MLEM-Algorithmus genannt.

M. Oehler und T.M. Buzug

Abb. 1. (a) Experimenteller Philips Tomoscan M/EG mit CIRS Torsophantom markiert mit zwei Stahlmarkern; (b) Sinogramm der aufgenommen Daten des Torsophantoms markiert mit zwei Stahlmarkern; (c) Vergr¨ oßerter Ausschnitt der Projektion im Bereich der Metallprojektion unter dem Winkel γ = 120◦ einschließlich der Ergebnisse der linearen und polynomialen Interpolationen sowie der Ground Truth (gemessen ohne Stahlmarker).

¨ Die Herleitung des neuen λ-MLEM-Algorithmus basiert auf zwei Anderungen der urspr¨ unglichen Formel [2]. Im ersten Schritt werden alle Zeilen der Systemmatrix A, die zu Projektionen geh¨ oren, die durch ein Metallobjekt verlaufen, mit dem Vertrauensparameter 0 ≤ λ ≤ 1 gewichtet, so dass die neue Log-LikelihoodFunktion f¨ ur die Transmissions-Computertomographie wie folgt lautet

l (f ∗ ) = ln (L (f ∗ )) =

⎛ ⎝−˜ ni

i=1

⎞ λi aij fj∗ − n0 e−

∗ j=1 λi aij fj

⎠ + const., (1)

j=1

˜ i der wobei n0 der Anzahl der Photonen, die die R¨ontgenr¨ohre verlassen und n modiﬁzierten Anzahl der detektierten Photonen entspricht. Der urspr¨ ungliche MLEM-Algorithmus ist ein einfacher Gradienten-Algorithmus [7], der auch nach Modiﬁkation in folgender Form dargestellt werden kann (2) f ∗(n+1) = f ∗(n) + D f ∗(n) grad(l (f ∗ )); hierbei entspricht D der Diagonalmatrix f¨ ur die Transmissions-Computertomographie ∗(n) fr ∗(n) D f , (3) = diag M ˜ i λi air i=1 n so dass Gleichung (2) wie folgt geschrieben werden kann fr∗(n+1) =

fr∗(n)

∗(n)

+ M

∂l (f ∗ ) ∂fr∗

˜ i λi air i=1 n M ∗(n) N ∗(n) fr n0 λi air e− j=1 λi aij fj M ˜ i λi air i=1 i=1 n

−

M i=1

n ˜ i λi air

(4) .

Maximum-Likelihood-Ansatz zur Metallartefaktreduktion Abb. 2. (a) FBP-Rekonstruktion, (b) MLEM-Rekonstruktion; (c) Rekonstruktion mit λi = 0.5 (innerhalb der Metallprojektionen).

λ-MLEM-

¨ ¨ Der zweite Anderungsschritt beinhaltet die Uberlegung, dass die Anzahl der gemessenen R¨ ontgenquanten im Detektor proportional zur Intensit¨at der R¨ ontgenstrahlung ist. Aus diesem Grund muss die Projektionssumme pi =

aij fj ,

(5)

j=1

das heißt, die gemessenen Projektionswerte pi im Sinogram ebenfalls ad¨aquat gewichtet werden. Die angepasste Anzahl der R¨ontgenquanten ist damit n ˜i = n0 e−λi pi . Somit lautet die neue Fix-Punkt-Iteration fr∗(n+1)

fr∗(n)

∗(n)

λi air e− j=1 λi aij fj M −λi pi i=1 λi air e

i=1

(6)

In Abb. 2a ist das Ergebnis der FBP-Rekonstruktion dargestellt. Abbildung 2b zeigt das Ergebnis der Rekonstruktion mit dem klassischen MLEMAlgorithmus [7] basierend auf der modiﬁzierten MLEM-Gleichung (6). Diese ist n¨ amlich identisch mit der urspr¨ unglichen MLEM-Rekonstruktionsvorschrift bei ur alle Projektionen i ∈ {1, . . . , M }. Innerhalb der beiden einer Wahl von λi ≡ 1 f¨ unterschiedlich rekonstruierten Bilder 2a/b sind die Metallartefakte gut zu erkennen. Dies beruht auf Inkonsistenzen innerhalb der unbehandelten Rohdaten, die durch Aufh¨ artung hervorgerufen werden. Im Vergleich hierzu zeigt Abb. 2c das Ergebnis der zun¨ achst mit polynomialer Interpolation reparierten Rohdaten und anschließender Rekonstruktion mit dem λ-MLEM-Algorithmus. Die mit dem λ-MLEM-Algorithmus rekonstruierten Bilder wirken im Vergleich zu den mit dem FBP-Algorithmus berechneten Bildern verrauschter. Daher ist die Erweiterung der Rekonstruktionsvorschrift des λ-MLEM-Algorithmus durch einen Regularisierungsterm erforderlich.

Ergebnisse und Diskussion

Als Beurteilungskriterium f¨ ur die G¨ ute der Metall-Artefakt-Unterdr¨ uckung wurde der Korrelationskoeﬃzient zwischen den Bildern jeweils mit und ohne Stahlmarker - also der Ground Truth - berechnet. Das beste Ergebnis wurde mit dem

M. Oehler und T.M. Buzug

λ-MLEM-Algorithmus mit einer Wahl von λ = 0.5 erzielt (Abb. 2c). Es zeigt sich, dass die Metallartefakte im Vergleich zur Rekonstruktion mit der StandardFBP gut reduziert werden k¨ onnen. Derzeit ist der λ-MLEM-Algorithmus in otigt bei 500 Iterationen auf einem Intel MATLABTM implementiert und ben¨ Pentium 4 Prozessor mit 3 GHz ca. 20 Minuten (Radonraum: 256x180 Daten; Bildraum: 256x256 Pixel). Auf Grund dieser langen Rechenzeit planen wir eine Implementation auf dedizierter Hardware, die Unterteilung der Projektionen in Ordered Subsets und die Einf¨ uhrung von Relaxationsmethoden. Weiterhin werden die entwickelten Verfahren mit anderen Algorithmen zur Metallartefaktreduktion verglichen. Hierbei stellen Methoden, die die polychromatischen Eigenschaften der R¨ ontgenstrahlung in der Computertomographie ber¨ ucksichtigen, einen besonderen Schwerpunkt dar.

Literaturverzeichnis 1. Kachelrieß M. Reduktion von Metallartefakten in der R¨ ontgen-ComputerTomographie. Ph.D. thesis. Friedrich-Alexander-Universit¨ at Erlangen- N¨ urnberg; 1998. 2. Lange K, Carson R. EM reconstruction algorithms for emission and transmission tomography. J Comput Assist Tomog 1984;8:306–316. 3. Jank J, Backfriede W, Bergmann H, Kletter K. Analyse des Konvergenzverhaltens von Rekonstruktionsalgorithmen anhand lokaler und globaler Parameter. Med Phys 2001;11:246–254. 4. Glover GH, Pelc NJ. An algorithm for the reduction of metal clip artifacts in CT reconstructions. Med Phys 1981;8:799–807. 5. Kalender WA, Hebel R, Ebersberger J. Reduction of CT artifacts caused by metallic implants. Radiology 1987;164:576–577. 6. Lewis RM, Bates RHT. Image reconstruction from projections III: Projection completion methods (theory). Optik 1978;50:189–204. 7. Lange K, Bahn M, Little R. A Theoretical Study of Some Maximum Likelihood Algorithms for Emission and Transmission Tomography. IEEE Trans Med Imaging 1987;6:106–114.

Ultraschall-Perfusionsbildgebung fu ¨r die Schlaganfalldiagnostik auf Basis eines Modells fu ¨r die Destruktionskinetik von Kontrastmittel Christian Kier1 , Karsten Meyer-Wiethe2 , G¨ unter Seidel2 3 und Til Aach 2

1 Institut f¨ ur Signalverarbeitung, Universit¨ at zu L¨ ubeck, 23538 L¨ ubeck Klinik f¨ ur Neurologie, Universit¨ atsklinikum S-H, Campus L¨ ubeck, 23538 L¨ ubeck 3 Lehrstuhl f¨ ur Bildverarbeitung, RWTH Aachen, 52056 Aachen Email: [emailprotected]

Zusammenfassung. Zur Diagnose von zerebrovaskul¨ aren Erkrankungen mittels Ultraschall ist neben den Methoden der Bolus- und Wiederanreicherungskinetik die Methode der Destruktionskinetik bekannt, bei der u orung von Kontrastmittel perfusionsbezoge¨ber die schnelle Zerst¨ ne Parameter bestimmt und farbkodiert als Parameterbilder dargestellt werden. Das dabei bislang verwendete Modell hat den Nachteil, diese Parameter nur durch ein Least-Squares-Fitting bestimmen zu k¨ onnen. Eine Weiterentwicklung des Modells f¨ uhrt zu der M¨ oglichkeit, den sogenannten Perfusionskoeﬃzienten direkt aus den Daten berechnen zu k¨ onnen. Desweiteren wird eine verbesserte Sch¨ atzung des daf¨ ur notwendigen Destruktionskoeﬃzienten dargelegt. Eine Versuchsstudie, in der die erhaltenen Ergebnisse mit MRT-Aufnahmen verglichen werden, liefert deutliche Hinweise darauf, daß das neue Modell eine exaktere Bestimmung der Perfusion erm¨ oglicht.

Einleitung

Die zuverl¨ assige und schnelle Messung der zerebralen Mikrozirkulation ist entscheidend f¨ ur die Diagnose und Behandlung akuter zerebrovaskul¨arer Erkrankungen. Derzeit wird diese Messung vorwiegend mit Hilfe zeit- und kostenintensiver tomographischer Verfahren (CT/MRT) durchgef¨ uhrt. Die Messstationen sind dar¨ uber hinaus ortsgebunden, so daß Transfer und eventuelle Wartezeiten zum Teil erhebliche Belastungen des Patienten mit sich bringen. Mit der transkraniellen Sonograﬁe lassen sich kosteng¨ unstig und schnell Informationen u ¨ ber die lokale Hirndurchblutung am Bett des Patienten gewinnen. Dabei werden Ultraschall-Bildsequenzen aufgezeichnet, die den Verlauf der Konzentration eines Kontrastmittels im Blut festhalten, welches durch Ultraschallpulse in regelm¨ aßigen Abst¨ anden zerst¨ ort wird. Aus diesem Verlauf lassen sich perfusionsbezogene Parameter bestimmen, die farbkodiert als Parameterbilder dargestellt werden und R¨ uckschl¨ usse auf die Perfusion in einzelnen Bild- und damit Gehirnbereichen erlauben [1]. Bekannt sind eine bolusbasierte Methode (BHI [2]) sowie zwei auf konstanter Kontrastmittelinfusion basierende Methoden der Wiederanreicherungs- (RHI [3]) und Destruktionskinetik (DHI [4]).

C. Kier et al.

Harmonic Imaging und Destruktionskinetik

Um das bedingt durch die starke Schallreﬂektion am Sch¨adelknochen geringe Signal-Rausch-Verh¨ altnis von transkraniellem Ultraschall zu verbessern, macht sich Harmonic Imaging das nichtlineare Verhalten des Kontrastmittels zu Nutze. Ultraschall-Kontrastmittel besteht im Wesentlichen aus Mikrobl¨aschen, die bei Beschallung den Schall nichtlinear reﬂektieren und abh¨angig vom Schalldruck dabei zerst¨ ort werden. W¨ ahrend der Puls mit einer Mittenfrequenz von 1,8MHz gesendet wird, werden im Echo auch die entstehenden harmonischen Schwingungen erfasst, die haupts¨ achlich auf das im Blut vorhandene Kontrastmittel zur¨ uckzuf¨ uhren sind und somit vorrangig die Perfusion abbilden. Die Diminution Harmonic Imaging (DHI) Methode beruht auf konstanter intraven¨ oser Kontrastmittelinfusion, die zu einer konstanten Kontrastmittelkonzentration im Blutkreislauf f¨ uhrt. Abbildung 1 (li.) zeigt exemplarisch das erste Bild einer DHI-Sequenz mit anfangs hohem Kontrastmittelniveau. Durch eine Folge von Ultraschallpulsen (1,67-6,67 Hz) wird die Kontrastmittelkonzentration vermindert, da durch den bei dieser Methode verwendeten Schalldruck die Mikrobl¨ aschen im Kontrastmittel degenerieren. Nach ca. 10 Pulsen wird ein Gleichgewichtszustand g erreicht, in dem bei jedem Puls genau soviel Kontrastmittel zerst¨ ort wird, wie im folgenden Interpulsintervall in die Bildebene einstr¨omt. Abbildung 1 (re.) verdeutlicht diesen Prozeß. Zuﬂuß, Abﬂuß und Zerst¨orung des Kontrastmittels k¨onnen f¨ ur jeden Zeitschritt und jede Volumeneinheit der Bildebene durch folgendes rekursives parametrisches Modell [5] beschrieben werden: c(n + 1) = c(n) · d · e−p· ∆t + c(0) · 1 − e−p · ∆t . Abf luss

(1)

Zuf luss

Dabei ist n die Anzahl der Ultraschallpulse, c(n) die Konzentration, ∆t die Zeit zwischen zwei Pulsen, d der Destruktionskoeﬃzient und p der Perfusionskoeﬃzient. Um die Konzentration im n-ten Schritt zu berechnen, kann folgende geschlossene Form verwendet werden: xn−1 − 1 xn−1 + y · , mit x−1 x = d · e−p · ∆t

c(n) = c(0) ·

(2)

y = 1 − e−p · ∆t .

Bisher wird das Modell mit Least-Squares Methoden an die Bilddaten der Ultraschallsequenz geﬁttet und so der Perfusionskoeﬃzient p errechnet, der letztlich Aufschluß u ¨ ber die Perfusion geben soll. Genauso wird der ortsabh¨angige Destruktionskoeﬃzient d bestimmt, der die Destruktionswirkung des Ultraschallpulses ortsaufgel¨ ost wiedergibt [5]. Durch das Fitting k¨onnen diese Parameter bislang nur mit eingeschr¨ ankter Genauigkeit gesch¨atzt werden [6]. Desweiteren muß die Tiefenabh¨ angigkeit der Destruktionswirkung des Ultraschallpulses auf das Kontrastmittel ber¨ ucksichtigt werden.

Ultraschall-Perfusionsbildgebung f¨ ur die Schlaganfalldiagnostik

Abb. 1. Exemplarisches Ultraschallbild einer DHI-Sequenz (links). Tats¨ achliche und beobachtete DHI-Kontrastmittelkonzentration abh¨ angig von US-Pulsen (rechts). Beobachtete Konzentration TatsÃ¤chliche Konzentration Gleichgewichtszustand g

Kontrastmittelkonzentration (normalisiert)

0.8

Destruktion 0.4

Messung

0.2

Perfusion

0.6

4 5 6 Ultraschallpuls

Weiterentwicklung des mathematischen Modells

Das sich einstellende Gleichgewicht zwischen Zuﬂuß, Abﬂuß und Destruktion resultiert in einem konstanten Intensit¨ atswert, der als Grenzwert g des zugrunde liegenden Modells (mit n → ∞) interpretiert werden kann. Die Grenzwertbildung verdeutlicht die Abh¨ angigkeit von g von der Ausgangskonzentration c(0), vom Perfusionskoeﬃzienten p sowie vom Destruktionskoeﬃzienten d: xn−1 − 1 n→∞ n→∞ n→∞ x−1 −y e−p · ∆t − c(0) xn−1 − 1 p>0 = c(0) · = = c(0) · lim y · n→∞ x−1 x−1 d · e−p · ∆t − c(0) (2)

g = lim c(n) = lim c(0) · xn−1 + lim c(0) · y ·

(3)

Zur Berechnung von p ist es zun¨ achst notwendig, g und d zu berechnen. Die Ausgangskonzentration c(0) wird aus dem ersten Bild der DHI-Sequenz bestimmt. Zur Berechnung von d wird davon ausgegangen, daß die Destruktion unmittelbar nach dem Ultraschallpuls auftritt. Im folgenden Intervall bis zum n¨achsten Puls kommen dann nur noch die Flußeﬀekte zum Tragen. Wenn also das Intervall hinreichend kurz gew¨ ahlt wird, kann d als Koeﬃzient von zwei aufeinander folgenden Werten betrachtet werden. Zur Berechnung werden deshalb von der Kurve mit der h¨ ochsten Aufnahmefrequenz (6,67 Hz) die ersten beiden Werte herangezogen, da sich das Kontrastmittelniveau hier noch nahe an der S¨ attigung beﬁndet und somit die Flußeﬀekte noch nicht so stark sind, die Destruktion vom Absolutwert jedoch am h¨ ochsten ausf¨allt: lim

∆t→0

c(1) (2) 1 lim c(0) · d · e−p · ∆t + c(0) · 1 − e−p · ∆t = d = c(0) c(0) ∆t→0

(4)

Nachdem der Grenzwert als Mittelwert der letzten drei Intensit¨atswerte gebildet wurde, kann p nach Umformung von (3) aus d, Grenzwert sowie Anfangskonzentration c(0) berechnet werden:

C. Kier et al.

Abb. 2. Parameterbilder des Perfusionskoeﬃzienten p (li.) und Destruktionskoeﬃzienzen d (re.). Je heller die Intensit¨ at, desto st¨ arker die Perfusion bzw. desto schw¨ acher die Destruktion. hom*ogen graue (li.) und schwarze (re.) Bereiche sind ung¨ ultige Pixel.

1 · ln p= ∆t

g · d − c(0) g − c(0)

(5)

Nach der Berechnung des Perfusionskoeﬃzienten werden die Werte normalisiert und farbkodiert in einem Parameterbild dargestellt. Dabei werden die Bildpunkte, bei denen der Grenzwert g > α · c(0) ist (mit α frei w¨ahlbar, i.A. 1) als ung¨ ultig markiert, da hier die Modellannahme nicht zutriﬀt. Abbildung 2 zeigt Parameterbilder von p und d. Als weitere Parameter wurden aus der Zeit-Intensit¨ atskurve die maximale Amplitude und die Halbwertszeit des Kontrastmittels errechnet und ebenfalls als farbkodiertes Parameterbild dargestellt. Alle Berechnungen erfolgen pixelbezogen und erlauben somit eine pixel- oder regionenbasierte semi-quantitative Auswertung.

Versuchsstudie an 8 Patienten

In einer laufenden Studie wurden von bisher 8 Patienten mit Schlaganfallsymptomen jeweils DHI-Sequenzen sowie MRT-Aufnahmen erstellt. In den MRTBildern wurden solche Bereiche als Perfusionsdefekte klassiﬁziert, in denen die Durchblutung mehr als 4 Sekunden verz¨ ogert auftrat. In den DHI-Parameterbildern wurden Perfusionsdefekte von klinischen Experten ohne Kenntnis der MRTBefunde von Hand markiert. Die Gr¨ oße der Perfusionsdefekte wurde mittels des Spearman Ranking Verfahrens zueinander in Beziehung gesetzt. Dabei zeigt sich, daß von den vier Parametern (p, d, Amplitude, Halbwertszeit) allein der Perfusionskoeﬃzient eine signiﬁkante Korrelation (r = 0.8, Signiﬁkanz p = 0.017) mit den aus MRT-Aufnahmen bestimmten Perfusionsdefekten aufweist.

Diskussion

Im Gegensatz zu den derzeit verwendeten Verfahren erm¨oglicht die vorgestellte Methode die kosteng¨ unstige und schnelle Diagnostik der Hirnperfusion aufgrund von Ultraschallsequenzen direkt am Patientenbett unter minimaler Bela-

Ultraschall-Perfusionsbildgebung f¨ ur die Schlaganfalldiagnostik

stung des Patienten. Hervorzuheben ist die kurze Dauer der Untersuchung mittels der DHI-Methode, die neben verminderter Patientenbelastung zu stabileren Untersuchungsbedingungen im Vergleich zu anderen Harmonic Imaging Methoden (Bolus- und Wiederanreicherungskinetik) f¨ uhrt und der Methode damit das gr¨oßte Potential verschaﬀt [4, 2, 5]. Die gezeigten Ergebnisse sind aufgrund des Umfangs der Studie nur vorl¨auﬁg, zeigen aber, daß unsere Weiterentwicklung des mathematischen Modells deutlich zur Verbesserung der DHI-Methode gef¨ uhrt hat und eine exaktere Berechnung des Perfusionskoeﬃzienten p erlaubt. Ebenfalls zur Verbesserung der Genauigkeit tr¨ agt die tiefenabh¨ angige Sch¨ atzung des Destruktionskoeﬃzienten d bei, da somit die Tiefenabh¨ angigkeit von p reduziert wird. Eine Schw¨ache der Methode ist derzeit noch die untersucherabh¨ angige Bestimmung der Perfusionsdefekte in den Parameterbildern, die allerdings durch ein standardisiertes und somit automatisches Verfahren behoben werden kann. Desweiteren kann davon ausgegangen werden, daß mit unserer Weiterentwicklung die DHI-Perfusionsbildgebung auch an anderen per Ultraschall erreichbaren Organen verbessert werden kann.

Literaturverzeichnis 1. Metzler Volker, Seidel G¨ unter, Wiesmann Martin, Meyer Karsten, Aach Til. Perfusion Harmonic Imaging of the Human Brain. In: Ultrasonic Imaging and Signal Processing. vol. 5035 of Procs SPIE. San Diego, CA; 2003. p. 337–348. 2. Kier Christian, Toth Daniel, Meyer-Wiethe Karsten, Schindler Angela, Cang¨ ur Hakan, Seidel G¨ unter, et al. Cerebral Perfusion Imaging with Bolus Harmonic Imaging. In: Ultrasonic Imaging and Signal Processing. vol. 5750 of Procs SPIE. San Diego, CA; 2005. p. 437–446. 3. Metzler V, Seidel G, Toth D, Claassen L, Aach T. Schnelle Messung der lokalen Hirnperfusion zur Diagnoseunterst¨ utzung bei zerebrovaskul¨ aren Erkrankungen. In: Procs BVM; 2001. p. 280–284. 4. Meyer Karsten, Seidel G¨ unter. Transcranial contrast diminution imaging of the human brain – a pilot study on healthy volunteers. Ultrasound Med Biol 2002;28:1433– 1437. 5. Eyding J, Wilkening W, Reckhardt M, et al. Contrast burst depletion imaging (CODIM): a new imaging procedure and analysis method for semiquantitative ultrasonic perfusion imaging. Stroke 2003;34(1):77–83. 6. Eyding Jens, Wilkening Wilko, Reckhardt Markus, et al. Reliability of semiquantitative ultrasonic perfusion imaging of the brain. J Neuroimaging 2004;14:143–149.

Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences Marc Hensel1 , Bernd Lundt2 , Thomas Pralow2 and Rolf-Rainer Grigat1 1 2

Vision Systems, TU Hamburg-Harburg, 21079 Hamburg General X-Ray, Philips Medical Systems, 22335 Hamburg Email: [emailprotected]

Abstract. We present a practice-oriented, i.e. fast and robust, estimator for strong signal-dependent noise in medical low-dose X-ray images. Structure estimation by median ﬁltering has shown to be superior to linear binomial ﬁltering. Falsiﬁcations due to remaining structure in the estimated noise image are signiﬁcantly reduced by iterative outlier removal.

Introduction

Medical X-ray image sequences as applied, e.g., in angiography exhibit severe signal-dependent noise. This is founded on low X-ray doses used for radiation protection reasons. Noise reduction as well as contrast enhancement methods can be signiﬁcantly improved by consideration of the signal-dependent noise level. However, in noise estimation several intricacies arise. First of all, the application requires real-time processing while sophisticated noise estimators exhibit high computing time (complexity). Also, there is the problem of structure diversity. Noise estimation implies signal estimation (object structures). However, clinical sequences contain large regions of low contrast and smooth transitions (e.g. organs) as well as ﬁne structures of high contrast (e.g. catheters, vascular trees). Finally, existing noise models are not universally applicable for all observed X-ray images. Existing noise models are based on the Poisson distribution as noise in lowdose X-ray images is dominated by statistical variability of the X-ray quanta and, thus, by additive Poisson-distributed quanta noise. However, all too often such distribution cannot be observed in medical sequences. Flat-panel detectors might consist of multiple detector arrays exhibiting signiﬁcantly diﬀering properties (e.g. sensitivity) and, hence, require corrective preprocessing. Moreover, device-speciﬁc logarithmic mappings of digital X-ray images with the purpose to establish linear dependency of gray-values from the thickness of imaged objects are applied. These mappings invert the monotony of signal-dependent noise curves. Even worse, mapping parameters not matching the detector or camera type can totally degrade the noise curve. As concluding example, preprocessing by Wiener ﬁltering leads to a sharp bend in noise curves that cannot be described by the known models.

Robust and Fast Estimation of Signal-Dependent Noise

State-of-the-Art

State-of-the-art noise models for low-dose X-ray images are based on the assumption of the data being dominated by additive Poisson-distributed quanta noise. In Poisson-distributed data, mean µ (i.e. the expected uncorrupted signal inten√ sity s(x, y)) equals variance σ 2 . Thus, noise standard deviation σ = µ increases with the square root of the expected signal strength. Most often, and leading to a diﬀering model, the exponential attenuation of X-rays passing through matter is compensated by logarithmic mapping. Hereby, linear dependency of image intensity and traversed object thickness is established. At the same time, logarithmic mapping fundamentally changes the noise characteristic. The resulting noise standard deviation is monotonically decreasing over large signal ranges and can be modeled by an exponential function. However, both models are devicespeciﬁc and preprocessing might yield observed noise characteristics that are not adequately described by these models. There are a multitude of noise estimation methods [1]. Many of these are not suited for real-time X-ray image processing, because of the huge amount of data, the complexity of the method, or because they are not suited for signal-dependent noise. Examples are methods based on wavelet transformation or Maximum Likelihood estimation [2, 3]. Noise estimators for real-time applications are mostly based on signal estimation by low-pass ﬁltering. The high-pass portion, i.e. the data ﬁltered out, is interpreted as noise and analyzed to yield the noise estimation [1, 4]. Naturally, the resulting noise estimation can be considerably falsiﬁed by structure, as edges contain signiﬁcant high-frequency components.

Methods

In the context of real-time noise estimation, three major questions arise: What are the noise properties? Which ﬁlter is suited to approximate the signal? Finally, how can the inﬂuence of signal on the estimated noise image be further reduced? 3.1

Noise Model

The physical process underlying noise in low-dose X-ray images is Poissondistributed quanta noise. For a – in practice fulﬁlled – ”suﬃciently large”number of quanta contributing per pixel (about 20), the discrete distribution can be approximated by the continuous normal (Gaussian) distribution with same mean and variance. Our measurements of X-rayed brass wedges using varying preprocessing steps veriﬁed that the observed signal-dependent noise can be described fairly well by additive zero-mean normal-distributed noise η with signaldependent standard deviation ση (s), i.e. g(x, y) = s(x, y) + η s(x, y) (1) 2 η 1 1 exp − pdf(η; s) = √ (2) 2 σ 2π · ση (s) η (s) with signal s and the observed gray-valued image g.

M. Hensel et al. Table 1. Noise estimation in hom*ogeneous images using 3 × 3 standard ﬁlters Filter Median Binomial

3.2

σse /σ0 0.4078 0.3750

σηe /σ0 0.9720 0.8004

ηe [%] 2.85 19.96

Signal Estimation Filter

With estimated signal se , the estimated noise component at location (x, y) is given by ηe (x, y) = g(x, y)−se (x, y). For real-time signal estimation in very noisy images, 3 × 3 median and (linear) binomial ﬁlters h are taken into account [1]. In hom*ogeneous signals s(x, y) = s0 corrupted by noise σ0 , linear ﬁltering se = h∗g yields an estimated noise image ηe with standard deviation ση2e according to: ⎡ ⎤ ση2e = ⎣ [1 − δ(i)δ(j)] h2i,j + (1 − h0,0 )2 ⎦ · σ02 (3) i

The performance of median ﬁltering has been analyzed by experiments. For both, median and binomial ﬁltering, Table 1 shows the noise levels in se and ηe as well as the relative error ηe = |σηe − σ0 |/σ0 of the resulting noise estimation. Regarding se , noise reduction performance of the median ﬁlter is only slightly below the binomial ﬁlter. However, most remarkably, noise in ηe reﬂects σ0 by far better than in case of the binomial ﬁlter. Furthermore, median ﬁltering is superior in structure preservation and, thus, in limiting the inﬂuence of signal structure on ηe . 3.3

Intra-frame Noise Estimation

The proposed method is organized in ﬁve steps: 1. Signal estimation se (x, y) by 3 × 3 median ﬁltering and creation of the corresponding noise image ηe (x, y) = g(x, y) − se (x, y) 2. Subdivision of the dynamic range of the signal in intervals Si of equal size 3. Estimation of noise σηe (i) in the intervals S(i) as standard deviation of all noise values ηe (x, y) whose corresponding signal values se (x, y) are in Si 4. Iterative outlier removal 5. Interpolation with consideration of the statistical basis The quality of the initial noise estimation in S(i) %−1/2 $ 2 σηe (i) = E (Ni − E{Ni }) with & & Ni = ηe (x, y) && se (x, y) ∈ Si

(4) (5)

is improved by iterative outlier removal: In the calculation of the noise standard deviation σηe (i) for an interval S(i), absolute noise values above three times the

Robust and Fast Estimation of Signal-Dependent Noise

standard deviation are with high probability caused by high-frequency components originating from structure and not from noise. Thus, these outliers are removed from Ni (t) iteratively and the standard deviation is updated from the remaining values in Ni (t + 1): & & Ni (t + 1) = ηe (x, y)&& ηe (x, y) ∈ Ni (t) ∧ |ηe (x, y)| ≤ 3σηe (i) (6) The noise estimations σηe (i) for intervals Si are mapped to a signal-dependent estimation σηe (s) by interpolation of σηe (i) of those intervals, where a ”suﬃciently high”number (> 100) of pixel values have contributed to the standard deviation. Clinical X-ray images most often contain smooth transitions. Hence, normally most intervals are considered in the interpolation for the number of intervals (e.g. 128) and image data used by our groups. 3.4

Incremental Inter-frame Noise Estimation

For applications with very restrictive real-time requirements, the computing time for noise estimation can be controlled adaptively by incremental inter-frame noise estimation. Founded on the expression σ 2 (X) = E(X 2 ) − E(X)2 , the basic idea is to incrementally improve the estimation by additional consideration of some pixels of each new frame in a sequence. Number and location of the pixels is arbitrary, however, locations should diﬀer from image to image in such way that all image coordinates have the same contribution in the long run. Incremental adding can be applied as long as the imaging parameters are constant and, hence, frames exhibit equivalent noise characteristics.

Results

Signal-dependent noise in low-dose X-ray images is described well by additive zero-mean Gaussian noise with signal-dependent standard deviation. The maximum absolute diﬀerence of normal and Poisson distributions of equal mean and variance is solely 0.010 units (normed to signal values: 0.10%) for 10 quanta per pixel and 0.005 units (0.02%) for 20 quanta per pixel. The measured noise curves σηe (s) are quasi-continuous, but not necessarily monotonic. The precision of the noise estimation has been evaluated using uniform, artiﬁcial, and clinical test images corrupted by noise with known characteristic. The clinical image showing a vascular tree has been acquired with comparatively high dose and, thus, contains low noise. This was further reduced by nonlinear diﬀusion ﬁltering prior to evaluation. On all test data, signal estimation by median ﬁltering clearly outperformed estimation based on binomial ﬁltering. In the absence of image structure, the relative error of the estimated noise was about 2.85% for median and 19.96% for binomial ﬁltering (Tab. 1). The presence of artiﬁcial structure falsiﬁed the noise estimation and yielded relative errors of 52.18% (median) and 81.69% (binomial). In the presented method, these errors were signiﬁcantly reduced, in

M. Hensel et al.

Fig. 1. Estimated noise histograms for binomial ﬁltering of clinical image. Ideal case: 1-D histograms for constant se are Gaussian-distributed. Left: All pixels contribute. Structure produces errors with predominantly negative/positive sign for small/large se (MSE: 8.96). Right: Only 37.16% of pixels with lowest entries (≤ 1% max. value) in edge image (Prewitt) contribute. Strong disturbances due to structure remain (MSE: 7.54). Proposed method: The disturbances are largely removed (MSE: 3.34). (For better visibility the data is normed to 1 for each se and the display clipped to 0.05.)

particular by structure iteration, to about 0.05% (median) and 17.38% (binomial), respectively. For the clinical image, relative estimation errors were about 2.67% and 18.91%, respectively. State-of-the-art gradient-based methods [1] did not yield comparable reduction of the inﬂuence of structure for the clinical test image (see Figure 1 and mean squared errors of normed data therein).

Discussion

The presented method is utmost robust in several aspects. Most prominent are its universality and suppression of object structures. The method adapts to the noise characteristics and is suited for exposures and ﬂuoroscopic sequences, likewise. Moreover, typical falsiﬁcations in state-of-the-art noise estimators due to highfrequency components of structure are signiﬁcantly reduced by iterative outlier removal. Finally, it has shown that signal estimation using median ﬁltering yields good results for images corrupted by severe noise (e.g. ﬂuoroscopic images).

References 1. Olsen SI. Estimation of noise in images: An evaluation. CVGIP: Graphical Models and Image Processing 1993;55(4):319–323. 2. Yuan X, Buckles BP. Subband Noise Estimation for Adaptive Wavelet Shrinkage. In: Proc IEEE 17th Int Conf Pattern Recognition. vol. 4; 2004. p. 885–888. 3. Sijbers J, den Dekker AJ. Maximum Likelihood Estimation of Signal Amplitude and Noise Variance From MR Data. Magn Reson Med 2004;51(3):586–594. 4. Amer A, Dubois E. Fast and Reliable Structure-Oriented Video Noise Estimation. IEEE Trans Circuits and Systems for Video Technology 2005;15(1):113–118.

CLARET: A Tool for Fully Automated Evaluation of MRSI with Pattern Recognition Methods B.M. Kelm1 , B.H. Menze1 , T. Neﬀ2 , C.M. Zechmann2 and F.A. Hamprecht1 1

Interdisciplinary Center for Scientiﬁc Computing (IWR), University of Heidelberg, 69120 Heidelberg 2 German Cancer Research Center (dkfz), 69120 Heidelberg Email: [emailprotected]

Abstract. Magnetic Resonance Spectroscopic Imaging (MRSI) measures relative concentrations of metabolites in vivo and can thus be used for the diagnosis of certain tumors. We introduce the program CLARET that makes MRSI accessible for clinical routine use. Instead of embarking on an error-prone quantiﬁcation of metabolites that requires manual checking of the results in many voxels, the program uses pattern recognition methods to directly compute tumor probability maps. Furthermore, non-evaluable signals are identiﬁed and masked out. The user can thus save time and concentrate on suspicious regions only.

Introduction

With magnetic resonance spectroscopic imaging (MRSI), speciﬁc metabolites can be detected and spectrally resolved in vivo. Since the concentration ratios of certain metabolites change characteristically in pathologic tissue, MRSI is in principle very well suited for the detection and localization of tumors. A major challenge in MRSI, however, lies in the postprocessing and evaluation of the acquired spectral volumes. In general, two approaches to the evaluation of MRSI have been considered in the literature. The ﬁrst approach initially performs a quantiﬁcation based on physical signal models [1]. Algortihms such as HLSVD, VARPRO, AMARES, LCModel, QUEST, etc. allow to incorporate diﬀerent kinds of prior knowledge in order to obtain robustness towards noise. Such or similar algorithms are often provided by commerical vendors (e.g. Siemens – PRISMA, Philips – LCModel) and are also available in the popular Java solution jMRUI. However, in order to avoid wrong diagnoses because of erroneous quantiﬁcation results, one has to manually verify the line ﬁt in every single voxel. Merely because of time constraints, this approach is not viable in clinical routine, in particular at higher spatial resolutions. Furthermore, manual evaluation of MRSI data requires tissue speciﬁc domain knowledge not shared by all radiologists. Finally, such a manual evaluation lacks objectivity and reproducibility.

B.M. Kelm et al.

The second approach to the evaluation of MRSI data applies methods from machine learning and pattern recognition [2]. From this point of view, quantiﬁcation is only one way of feature extraction and dimension reduction. Instead of using error-prone quantiﬁcation, pattern recognition methods can as well use any other, possibly more robust, representation of the spectral information. Several techniques have been proposed in the literature among which are e.g. artiﬁcial neural networks, support vector machines, subspace methods, generalized linear models, etc. Recent results [3, 4] show that pattern recognition methods which do not rely on an explicit quantiﬁcation step can be superior for a given diagnostic problem. Furthermore, pattern recognition also allows to identify spectra which are dominated by artifacts or which are just too noisy to be processed reliably [5]. Despite their considerable potential, pattern recognition methods for the evaluation of MRSI are not applied in clinical routine. One important reason is that commercial vendors do not oﬀer or support such evaluation strategies yet. Another reason is that no software tools are available which oﬀer pattern recognition methods for the evaluation of MRSI in a user-friendly way. In this paper we introduce the CLARET tool (C SI-based Localization And Robust E stimation of T umor probability) for the diagnostic evaluation of MRSI data. CLARET implements powerful pattern recognition methods for an automatic evaluation of MRSI volumes. Furthermore, it provides convenient tools and a user-friendly interface for time-eﬃcient analyses.

Design Principles and Use Cases

The evaluation of MRSI volumes with CLARET is designed for utmost user friendliness. After selecting an MRSI volume and a suitable MR image volume (usually T2 -weighted) from the DICOM data set, CLARET can be initiated to evaluate either individual slices or the whole loaded volume at once. The results are displayed in transparent probability maps superposed onto slices through the MRI volume (Fig. 1). One can easily switch between tumor probability estimates and their 2σ conﬁdence intervals. In addition, voxels which cannot be evaluated are masked out. In the subsequent diagnosis the user can therefore concentrate only on regions marked suspicious in the probability map. In case of doubt, the original spectral signal is easily accessible and conspicuities in the T2 -image can also be scrutinized. Finally the extracted probability map can be stored in a ﬁle together with the analyzed MRI/MRSI volumes for later reference or it can be exported for use in the radiation planning software VIRTUOS [6]. CLARET has explicitly been designed for the application of pattern recognition methods. Therefore, it can also be used for the construction of training data sets. The automatic display of the respective spectral signals together with ﬁtted model spectra and quantiﬁcation results upon selection of an MRSI voxel allow for a semi-manual evaluation of the spectral data. The results from such a voxel-wise evaluation can easily be entered per mouse click in the probability map and stored as training data set. Since manual labels can also be entered

CLARET: A Tool for Fully Automated Evaluation of MRSI

Table 1. Distribution of labels in the prostate data set (76 slices from 24 patients). quality \ class non-evaluable poor good total

healthy – 721 1665 2386

undecided – 437 629 1066

tumor – 284 452 736

total 15268 1442 2746 19456

after an automatic evaluation CLARET is also suitable for the correction of classiﬁcation errors and is ready for active learning.

Experimental Results

Spectral recordings from an ongoing clinical prostate MRSI study (IMAPS, [7]) from 24 patients with 16 × 16 × 16 voxels and 512 spectral channels have been included in the assessment. The data has been acquired on a common MR scanner (MAGNETOM Symphony, SIEMENS Medical Systems, Erlangen) with endorectal coil. Altogether 76 MRSI slices have been labeled with CLARET with respect to signal quality and diagnostic classiﬁcation (cf. Table 1). In this way it was possible to label about 20,000 voxels in a relatively short amount of time and to store them as training data set in a machine-readable format. The analysis of the obtained data has shown that pattern recognition methods which do not employ an explicit quantiﬁcation step for feature extraction can yield better performance. For the discrimination between tumorous and healthy tissue all considered pattern recognition methods reached a median cross-validated area under ROC (receiver operator characteristic) of at least .995 and therefore performed at least as good as the best quantiﬁcation approaches (cf. Fig. 2, left). The right part of Figure 2 illustrates the class separation that can be obtained for the example of a binomial PLS (partial least squares, ”binom PLS (m)”). The noise classiﬁer discriminates evaluable from non-evaluable spectra with an estimated out-of-bag error of 4.9%. With the current choice of pulse sequence and ﬁeld of view, approximately three quarters of the MRSI voxels have to be discarded by the noise classiﬁer as they are either heavily distorted by artifacts or contain noise only.

Discussion

For the ﬁrst time a software tool is available which generates pathophysiologic probability maps from MRSI data fully automatically. CLARET is currently employed in a prostate study at the German cancer research center Heidelberg (dkfz). The graphical user interface and integrated workﬂow allow for an eﬃcient evaluation of MRSI. Direct import of DICOM data from the MR scanner and the subsequent fully automatic evaluation by means of powerful pattern recognition

B.M. Kelm et al.

Fig. 1. The CLARET GUI can be used to evaluate MRSI eﬃciently. In routine use, the program automatically computes and displays tumor probability maps and conﬁdence intervalls on top of morphologic MR images. The program also allows for a point-andclick display of spectral raw data, it can perform quantiﬁcation, and it may be used for the manual labeling or the semi-manual reﬁnement of training data sets.

healthy

tumor 4

density

10 12

(Cho+Cr)/Ci (v) Ci/(Cho+Cr+Ci) (v) gauss (v) binom (v) RF (v) SVM−rbf (v) (Cho+Cr)/Ci (a) RF (a) SVM−rbf (a) (Cho+Cr)/Ci (q) RF (q) SVM−rbf (q) binom (m) binom PCR (m) binom PLS (m) binom GPLS (m) binom PSR (m)

undecided

0.85

0.90

0.95 AUC

1.00

0.0

0.2

0.4

0.6

0.8

1.0

Fig. 2. Left: Comparison of the tumor prediction capability of various approaches: (v)VARPRO, (a)-AMARES(jMRUI), (q)-QUEST(jMRUI) quantiﬁcation versus spectral (m)agnitude pattern recognition approaches. The area under curve (AUC) of the receiver operator characteristic has been 8-fold cross-validated. Right: Density estimates of the cross-validated tumor probability estimates.

CLARET: A Tool for Fully Automated Evaluation of MRSI

algorithms make its use simple. An application of CLARET for radiation therapy planning is projected and enabled by the integration into the software platform VIRTUOS [6]. The mentioned IMAPS prostate tumor study shows the potential of pattern recognition methods for the diagnostic evaluation of MRSI and highlights the necessity to provide suitable software. CLARET prototypically demonstrates the possibilities of a pattern recognition based MRSI evaluation. Here CLARET clearly contrasts with other MRSI evaluation tools such as jMRUI, LCModel or PRISMA which concentrate on the quantiﬁcation of spectral data. Although these programs can also be used to compute and visualize color maps, only relative metabolite concentrations or ratios thereof are displayed. In contrast, CLARET is tailored towards the generation and visualization of pathophysiologic maps that report an explicit estimate for tumor probability in every voxel.

References 1. Vanhamme L, Sundin T, Van Hecke P, et al. MR spectroscopy quantitation: a review of time-domain methods. NMR Biomed 2001;14:233–246. 2. Hagberg G. From magnetic resonance spectroscopy to classiﬁcation of tumors. A review of pattern recognition methods. NMR Biomed 1998;11:148–156. 3. Kelm BM, Menze BH, Zechmann CM, et al. Automated Estimation of Tumor Probability in Prostate MRSI: Pattern Recognition vs. Quantiﬁcation. Magn Reson Med submitted;. 4. Menze BH, Lichy MP, Bachert P, et al. Optimal Classiﬁcation of Long Echo Time In Vivo Magnetic Resonance Spectra in the Detection of Recurrent Brain Tumors. NMR Biomed submitted;. 5. Menze BH, Kelm BM, Heck D, et al. Machine based rejection of low-quality spectra, and estimation of tumor probabilities from MRSI. In: Procs BVM; 2006. . 6. Bendl R, Pross J, Hoess A, et al. VIRTUOS - A Program for Virtual Radiotherapy Simulation and Veriﬁcation. In: Proc of 11th Int. Conf. on The Use of Computers in Radiation Therapy. A. R. Hounsell u. a. Manchester: North Western Med. Physics Dept.; 1994. p. 226–227. 7. Scheenen T, Weiland E, Futterer J, et al. Preliminary Results of IMAPS: An International Multi-Centre Assessment of Prostate MR Spectroscpoy. In: Procs 13th ISMRM, Miami. Springer; 2005. .

Analysis of the Left Ventricle After Myocardial Infarction Combining 4D Cine-MR and 3D DE-MR Image Sequences D. S¨ aring1 , J. Ehrhardt1 , A. Stork2 , M.P. Bansmann2 , G.K. Lund3 and H. Handels1 1 Department of Medical Informatics, Department of Diagnostic and Interventional Radiology, 3 Department of Cardiology/Angiology, University Medical Center Hamburg-Eppendorf, 20246 Hamburg, Germany Email: [emailprotected] 2

Abstract. Spatial-temporal MR image sequences of the heart contain information about shape and motion changes and pathological structures after myocardial infarction. In this paper a software system called HeAT for the quantitative analysis of 4D MR image sequences of infarct patients is presented. HeAT supports interactive segmentation of anatomical and pathological structures. Registration of Cine- and DE-MR image data is applied to enable their combined evaluation during the analysis process. Partitioning of the myocardium in segments enables the analysis with high local resolution. Corresponding segments are generated and used for inter/intra patient comparison. Quantitative parameters were extracted and visualized. Parameters like endocard movement in the infarcted area of 6 infarct patients were computed in HeAT. Parameters in the infarct area show the expected dysfunctional characteristics. Based on theses parameters passive endocardial movement and myocardial areas with decreased contraction were identiﬁed.

Introduction

Myocardial infarction is one of the most common diseases in Germany. The infarcted area loses its ability to contract. This is attributed to left ventricle (LV) remodeling. LV remodeling is an important element in the progression of cardiac insuﬃciency characterized by wall thinning and chamber dilation for example [1]. It is a clinical problem to get quantitative parameters characterizing LV remodeling. New MR imaging techniques provide a non-invasive evaluation of anatomical and functional (local and global) parameters of the beating heart. In this paper spatial-temporal 4D Cine-MR and spatial delayed enhancement (DE)MR image sequences are combined to extract quantitative image parameters after myocardial infarction (Fig. 1). In this paper a software system called HeAT for the quantitative analysis of 4D MR image sequences of infarct patients is presented. Recently several software tools for quantitative analysis of spatial-temporal cardiac MR data were

Analysis of the Left Ventricle After Myocardial Infarction

Fig. 1. Two frames of a Cine-MR image sequence in endiastolic (a) and endsystolic (b) phase. Infarct area can be identiﬁed in the corresponding DE-MR image (c).

(a)

(b)

(c)

presented. The computation of global and local parameters in Cine-MR is supported rarely. Typically local parameters are computed and visualized using the 17 segment model. These parameters are limited for detailed evaluation of LV remodeling. Recently studies documented the mechanism of LV remodeling. Especially motion and shape changes of the myocardium in the infarct border zone were analyzed. Here the border zone was veriﬁed in Cine-MR manually using e.g. contrast echocardiography or delayed enhancement data for visual assistant. Generally this manual analysis is user-dependent and limited for inter/intra patient comparison. Other studies have underscored the impact of microvascular obstruction (MO) on the development of LV remodeling [2]. Due to the facts above the inﬂuence of MO is still poorly understood.

Methods

Computation of quantitative parameters in cardiac MR data requires localization and combination of the anatomical and pathological structures. Typically over 200 spatial-temporal Cine-MR and 7-12 spatial DE-MR images are acquired for each patient during breath-hold sequences. Segmentation and registration are necessary before the automatic algorithms for parameter extraction can be applied to evaluate the degree of myocardial infarction. Manual segmentation is used to outline myocardial structures in Cine-MR. Here identiﬁcation of pathological structures is not possible, but in DE-MR these structures can be recognized. Therefore a data-driven method is applied in DEMR images for segmenting pathological structures using user-deﬁned seed points within the healthy myocardial area. Based on their intensities mean x¯intensity and standard deviation σintensity are calculated. Due to the hom*ogeneous intensity values of healthy myocardial tissue, voxels in the myocardium with intensities higher than (1) xthreshold = x¯intensity + 2 · σintensity are regarded as necrotic (Fig. 2).

D. S¨ aring et al.

Fig. 2. Visualization of the results of direct contour overlay before (a) and after (b) rigid registration. Segmentation of the infarct area using xthreshold is presented in (c).

(a)

(b)

(c)

The MO is identiﬁed as voxels with signal intensity lower than x¯intensity − 2 · σintensity in surrounding hyperenhanced myocardium [3]. The direct transfer of these contours to Cine-MR image sequences for their combined evaluation leads to two problems. Alignment: Cine-MR and DE-MR imaging require diﬀerent settings. These different settings lead to diﬀerent ﬁeld of views or diﬀerent patient positions. A rigid registration method is applied to align the datasets (Fig. 2). Correspondence: In the same spatial position one DE-MR and several temporal Cine-MR images exist. Combination of both requires the deﬁnition of correspondence. This problem is addressed using Mutual Information [4] as similarity measure. Thus, information of both MR datasets can be combined transferring the contours of the registered DE-MR to the corresponding Cine-MR image sequence. Furthermore all Cine-MR frames of each slice are registered using non-linear demon based registration [5]. Applying the resulting transformation to the corresponding DE-MR slice a moving DE-MR can be computed. Detailed cardiac analysis requires a high local resolution. Based on the anatomical contours myocardial segments are generated and numbered [6]. Segments of diﬀerent phases and diﬀerent patients but with the same segment number can be identiﬁed as corresponding. Thereafter local quantitative parameters like size, length and transmurality are computed automatically for each myocardial segment. Combination of end-systolic and end-diastolic contours enables the calculation of parameters characterizing endocardial motion and increase of the wall thickness. Also global parameters e.g. infarct size are calculated. We implemented HeAT utilizing ITK and VTK with a problem-oriented adaptation.

Results

HeAT was evaluated using Cine-MR and DE-MR datasets of six patients (#2−7) after myocardial infarction and one patient (#1) with a healthy heart. Two

Analysis of the Left Ventricle After Myocardial Infarction

Table 1. For datasets #1 − #4 size of HS and IS relative to the myocardial area [%] and x ¯ ± σ of EM and DT in mid-ventricular slice are represented. Additionally for #7 the results for microvascular obstruction (MO) area are shown patient #1 #2 #3 #4 #7

section HS HS IS HS IS HS IS HS IS MO

[%] 100, 0 77, 1 22, 9 83, 7 16, 3 94, 3 5, 7 76, 2 23, 8 2, 8

EM[mm] 10, 0 ± 2, 8 7, 2 ± 2, 1 1, 9 ± 1, 1 5, 6 ± 1, 9 4, 5 ± 0, 7 7, 4 ± 2, 2 6, 4 ± 0, 8 12, 3 ± 1, 9 4, 1 ± 2, 2 2, 9 ± 1, 3

DT[mm] 5, 3 ± 1, 5 2, 7 ± 1, 6 1, 0 ± 1, 0 5, 6 ± 2, 2 3, 4 ± 1, 3 3, 7 ± 2, 1 3, 1 ± 0, 6 7, 4 ± 1, 6 3, 4 ± 2, 1 1, 4 ± 0, 9

physicans deﬁned an apical, a midventricular and a basal slice for each patient and outlined the anatomical and pathological structures. During the analysis the mean x ¯ and standard deviation σ of endocardial movement (EM) and diﬀerence of thickness (DT) of healthy (HS) and infarcted (IS) section were compared (Tab. 1). The quantitative parameters of #1 describe the movement of a healthy heart. #2 and #5 reﬂects the expected degree of myocardial infarction. In the infarction area EM and DT are decreased (5 mm) due to the loss of its contractile function. In contrast, in #3 and #6 the diﬀerences of EM are low (< 1, 1 mm) and the diﬀerences of DT are high ( > 2, 2 mm ). This result can be interpreted as a passive myocardial contraction. In #4 the changes of EM and DT values are not signiﬁcant ( < 1, 0 mm ). Here, the infarct area is small and heart function is inﬂuenced minimally. In #7 MO was discovered. In the basal and mid-ventricular slice the MO values of EM and DT are high decreased and indicate the center of the infarct (Fig. 3). In apical slice characteristics of passive movement can be observed. The computation of quantitative parameters using HeAT is useful to get objective and observer-independent criteria for the characterization of LV remodeling.

Discussion

The automatic parameter extraction and their visualization proofed to be helpful for the evaluation of patient image data after myocardial infarction. The visualization of the analysis results reﬂect signiﬁcant changes in shape and motion in the patient data considered. The registration of Cine- and MR image sequences enables the combined analysis of motion and shape in infarct regions and healthy tissue. According to LV remodeling the diagnostic impact of the generated DE-MR sequences at diﬀerent phases of the heart cycle has to be evaluated in future. The manual segmentation was time-consuming and its re-

D. S¨ aring et al.

Fig. 3. Representation of the quantitative parameters for all myocardial segments. The absolute values mm of endocardial movement (circles) and the diﬀerence in wall thickness between end-diastolic and end-systolic phase (triangles) are displayed as well as the transmurality % of the infarct area (diamonds) and the MO (boxes).

sults were user-dependent. Therefore automatic or semi-automatic methods have to be developed to improve the segmentation process. In the next step a database of a population of infarct patients will be established with base-line (acute phase) and follow-up (approx. 6 month past acute) evaluations concerning left ventricular remodeling. This can be useful for predicting the course of disease of an acute infarct patient.

References 1. Zimmerman SD, Criscione J, Covell JW. Remodeling in myocardium adjacent to an infarction in the pig left ventricle. Am J Physiol Heart Circ Physiol Dec 2004;(287):H2697 – H2704. 2. Wu KC, Zerhouni EA, Judd RM, et al. Prognostic signiﬁcance of microvascular obstruction by magnetic resonance imaging in patients with acute myocardial infarction. Circulation 1998;(97):765–772. 3. Gerber BL, Garot J, et al. Accuracy of Contrast-Enhanced Magnetic Resonance Imaging in Predicting Improvement of Regional Myocardial Function in Patients After Acute Myocardial Infarction. Circulation Aug 2002;(106):1083–1089. 4. Wells W, Viola P, et al. Multi-modal volume registration by maximization of mutual information. Medical Image Analysis 1996;1(1):35–51. 5. Thirion JeanPhilippe. Non-Rigid Matching Using Demons. In: Proc. Int. Conf. Computer Vision and Pattern Recognition (CVPR ’96). Washington, DC, USA: IEEE Computer Society; 1996. p. 245. 6. S¨ aring D, Ehrhardt J, Stork A, Bannsmann MP, Lund GK, Handels H. HeAT Heart Analysis Tool. In: 50. Jahrestagung der Deutschen Gesellschaft f¨ ur Medizinische Informatik, Biometrie und Epidemiologie. Freiburg; 2005. p. 28–30.

Content Based Image Retrieval for Dynamic Time Series Data Birgit Lessmann1,2 , Tim W Nattkemper2 , Johannes Huth2 , Christian Loyek2 , Preminda Kessar3 , Michael Khazen3 , Linda Pointon3 , Martin O. Leach3 and Andreas Degenhard1 1

Theoretical Physics, Physiscs Department, University of Bielefeld, 33615 Bielefeld Applied Neuroinformatics, Technical Faculty, University of Bielefeld, 33615 Bielefeld 3 Clinical MR Research Group, Royal Marsden Hospital, Institute of Cancer Research, Sutton, Surrey, UK Email: [emailprotected]

Abstract. Content based image retrieval (CBIR) systems in the ﬁeld of medical image analysis are an active ﬁeld of research. They allow the user to compare a given case with others in order to assist in the diagnostic process. In this work a CBIR system is described working on datasets which are both time- and space-dependent. Diﬀerent possible feature sets are investigated, in order to explore how these datasets are optimally represented in the corresponding database.

Introduction

In clinical practice, the amount of digitised image data acquired with respect to diagnostic decision making is steadily increasing. In particular three or four dimensional acquisition strategies give rise to large amount of data to be processed. One example is the acquisition of a series of images utilising dynamic contrast enhanced magnetic resonance imaging (DCE-MRI), regarded as a powerful technique for the early detection of breast cancer [1]. In DCE-MRI the uptake of a contrast agent is monitored over a selected time interval to access characteristics of the uptake behaviour in tissue [2]. This allows for detecting abnormalities and, additionally, classifying the detected lesions by their uptake characteristics. Malignant tumours are typically characterised by a signal intensity curve showing a strong and fast uptake often followed by a washout. Contrarily, benign tumours typically possess a slight or moderate uptake behaviour of the contrast agent [2]. The acquired image series contains both, the temporal and the spatial information, resulting in a high dimensional data set. Software tools designed to assist the radiologist in his or her daily clinical work are therefore an active ﬁeld of research. One particular class of assistant tools are the content based image retrieval (CBIR) systems based on the concept of query by example (QBE), i.e. searching database entries similar to a given query. Here we present a CBIR system for the interactive analysis of DCE-MRI images.

2 2.1

B. Lessmann et al.

Methods Content Based Image Retrieval

CBIR systems are based on the concept of query by example (QBE). Given an example case as input, the system output are those cases which are, by a chosen measure, most similar to the query case. CBIR systems for various purposes are currently under development. One type of example are retrieval systems for histopathological applications, usually based on texture, or modality independent retrieval systems for automatic categorisation of image data [3]. 2.2

Feature Vector Computation

We evaluate datasets acquired with diﬀerent acquisition protocols, an additional challenge often unavoidable in daily clinical work. We consider a database containing DCE-MRI data with four and also ﬁve post-contrast images. We compare methods that diﬀer in the treatment of these varying datasets with respect to feature vector computation. During the retrieval step the feature vectors of the diﬀerent cases are compared. To achieve comparable feature vectors for our heterogenous database we apply diﬀerent methodologies. An established method to describe signal enhancement in DCE-MRI is the use of radiological metafeatures [4], which do not depend on the dimension of the time-series. However, meta features, in general, strongly depend on the particular application, motivating the search for feature vectors applicable in a more generalised context. In this work, we propose and compare possible substitutions for meta features. In the following paragraphs the sets of feature vectors analysed in this work are described in detail. Truncation In a ﬁrst approach we truncate a six dimensional time series by omitting the last time-point, so there is the same minimum number of timepoints for each sequence. Extrapolation In a second approach the ﬁve dimensional time series are extrapolated to six dimensional ones. Wavelet features Based on a Discrete Wavelet Transform (DWT) a more sophisticated approach is considered. Wavelet analysis is a method of data analysis to access the localised and scale dependent information in signals [5]. With respect to time series analysis, wavelet-based methods are quite established [6]. We have wavelet-transformed the ﬁve or six-dimensional time series and then truncated the achieved series of coeﬃcients to the length of ﬁve. I.e. in case of a six-dimensional time series one coeﬃcient is omitted. In contrast the truncation in the time domain, the information of the last time-point, if available, is still partly included in some of the remaining coeﬃcients. Here the Daubechies-4wavelet [5] is employed for our computations.

Content Based Image Retrieval for Dynamic Time Series Data

Radiological Meta-Features Another type of feature vectors is based on established kinetic meta features used in MR radiology. We have evaluated the CBIRsystem using the features percentage enhancement, washout ratio, steepest slope and time-to-peak enhancement as described in [4].

Database

Our analysis includes 27 cases, i.e. 27 time series of DCE-MR images. Each single image in this time series is acquired within an acquisition time of 90 sec. Two images are acquired before and four or ﬁve images after injection of the contrast agent. Only one pre-contrast image was included in the further analysis. In 10 cases only four post-contrast images are available resulting in a shorter time series of images. 14 cases are histopathologically classiﬁed as malignant and 13 as benign. Segmentation of tumour tissue was performed within a preprocessing step by applying the dynamic pyramid linking algorithm in 3d. This usually results in a decomposition of a tumour into several parts. The average time series computed for each of the latter are then stored as entries in the database,normalised with the pre-contrast value. The database comprises a total number of approximately 3000 entries.

Results

A statistical evaluation of the diﬀerent feature sets is performed. All segments of each tumour for each case are selected as a query example and the entries corresponding to the remaining cases as database content. For each given query example the number of retrieved malignant and benign segments are counted. For a query example classiﬁed as malignant the retrieved malignant segments are labelled as true positive (TP) and the retrieved benign ones as false positive (FP) and vice versa. From the values of TP and FP the precision of each query is computed. The four diﬀerent methodologies are then compared to each other in a statistical evaluation. The results for the best ten matches are shown in ﬁgure 1. For each tumour an average precision curve is computed. Figure 1 presents the mean areas under the precision curve. In the averaging process the results for the diﬀerent segments are weighted regarding the size of the segment. Since this analysis is based on only retrieving the best ten matches, recall values cannot be computed. Since the database comprises 1678 feature vectors describing malignant tissue and 1073 feature vectors describing benign tissue, the a priori probabilities of retrieving malignant or benign cases are 0.61 and 0.39 respectively. Regarding these values the results of malignant tissue are signiﬁcantly better then those of benign tissue. In fact, the results of the benign lesions are quite near to the a priori probability. As can be seen in ﬁgure 1 the retrieval results for the truncation and extrapolation procedure are quite similar. The wavelet based features show a very good performance in retrieving malignant tissue and benign cases. The meta features show a performance similar to those of the wavelet features in retrieving

B. Lessmann et al.

Fig. 1. Mean areas under the precision curve for the for diﬀerent methods. Left: malignant cases, Right: benign cases 1 truncated time series extrapolated time series wavelet features meta features

mean area under precision curve

0.8

0.6

0.4

0.2

0 malignant cases

benign cases

the malignant cases, but provide a performance inferior to the others in case of benign entries. However, regarding the statistical ﬂuctuations measured none of the feature sets can deﬁnitely be determined as the best suited one. A further observation is, that the standard deviations of the benign cases are much higher then those of the malignant cases. This is assumed to be related to the following observation. Starting from a database containing the segments of 25 malignant and 17 benign lesions, we count the number of lesions providing retrieval results exceeding the a priori probabilities of retrieving malignant or benign cases. The results are shown in Table 1. In most of the malignant cases the averaged area under the precision curve exceed the a priori probability. In case of benign lesions the retrieval results seem to be quite randomly distributed. This leads to the conclusion that malignant tissue presumably shows a more characteristic uptake kinetic than benign tissue. In [7] it was shown, that malignant tumours often contain regions showing a benign uptake characteristic. Thus a lot of segments in the database derived from malignant tumours are diﬃcult to be distinguished from the benign ones.

Discussion and Conclusion

We proposed and evaluated a CBIR system for DCE-MR images capable of handling inhom*ogeneous data sets. Diﬀerent techniques are proposed to accomplish this task, including a wavelet based scheme. In a statistical analysis we showed that the feature vector sets achieve quite comparable results. However, among the approaches to be compared with the meta features the wavelet based approach shows the smallest standard deviation with respect to malignant cases. This observation corresponds to the fact that, regarding this approach, all lesions except one achieve a retrieval result exceeding the a priori probability. Furthermore, an advantage of the wavelet based approach is that it does not require an a priori decision for truncating a selected time point or extrapolating an additional time point, to achieve hom*ogeneity necessary for database retrieval. In particular the latter observation is desirable in clinical practice. The observation that the retrieval results of benign tissue only slightly exceed the a priori values

Content Based Image Retrieval for Dynamic Time Series Data

Table 1. The number of lesions reaching a retrieval results higher than the a priori probability. method

malignant cases benign cases total: 25 total: 17 truncation 21 8 extrapolation 21 9 wavelet features 24 9 meta features 24 6

is assumed to be related to tumour inhom*ogeneity. As shown in [7] malignant tumours also contain tissue showing benign characteristics. Thus a signiﬁcant amount of tumour segments in the database are classiﬁed as malignant while showing benign signal enhancement. Thus, we conclude that the representation of a tumour in a database through the result of a segmentation procedure is insuﬃcient for retrieval tasks. Further analysis with respect to tumour features will be necessary to ﬁnd a representation scheme satisfying the requirements within the clinical practice. Acknowledgements. The authors would like to thank Oliver Schulz, MaxPlanck-Institute Dortmund, for providing software resources. The work was funded by the University of Bielefeld and the federal state Nordrhein-Westfalen, Germany.

References 1. UK MRI Breast Screening Study Advisory Group. Magnetic resonance imaging screening in women at genetic risk of breast cancer: imaging and analysis protocol for the UK multicentre study. Magn Res Imag 2000;18:765–776. 2. Heywang-Koebrunner SH, Beck R. Contrast-Enhanced MRI of the Breast. Berlin: Springer-Verlag; 1996. 3. Lehmann TM, Gueld MO, Deselaers T, Keysers D, Schubert H, Spitzer K, et al. Automatic Categorization of medical images for content-based retrieval and datamining. Comp Med ImagGraph 2005;29:143–155. 4. Szabo BK, Aspelin P, Wiberg MKristoﬀerson, Bone B. Analysis of kinetic and morphologic diagnostic criteria. Acta Radiologica 2003;44:379–386. 5. Daubechies I. Ten Lectures on Wavelets. SIAM: CBMS-NFS Series Appl. Math.; 1991. 6. Percival DB, Walden AT. Wavelet Methods for Time Series Analysis. Cambridge, UK: Cambridge University Press; 2000. 7. Twellmann Thorsten, Lichte Oliver, Nattkemper TimW. An Adaptive Tissue Characterisation Network for Model-Free Visualisation of Dynamic ContrastEnhanced Magnetic Resonance Image Data. IEEE Transactions on Medical Imaging 2005;24(10):1256–1266.

Visual Exploration of Pathology Images by a Discrete Wavelet Transform Preprocessed Locally Linear Embedding Claudio Varini1,2 , Birgit Lessmann1,2 , Andreas Degenhard1 , Volkmar Hans3 and Tim Nattkemper2 1 Theoretical Physics, University of Bielefeld, Bielefeld, Germany Applied Neuroinformatics Group, University of Bielefeld, Bielefeld, Germany, 3 Institute of Neuropathology, Evangelischesn Krankenhaus, Bielefeld, Germany Email: [emailprotected] 2

Abstract. The information content of large collections of histopathological images can be explored utilizing computer-based techniques that can help the user to explore the similarity between diﬀerent brain tumor types. To visually inspect the degree of similarity between diﬀerent tumors, we propose a combined approach based on the Discrete Wavelet Transform (DWT) and Locally Linear Embedding (LLE). The former is employed as a preprocessing utility, the latter achieves the dimensional reduction required for visualization.

Introduction

The automatic exploration of large collections of histopathological images is an actual ﬁeld of research in computer science applied to medicine. Computer-based techniques can allow the user to ﬁnd relations and similarities in the data. In particular the visual exploration of a dataset containing diﬀerent types of tumor can provide valuable information concerning classiﬁcation. Techniques for dimensional data reduction oﬀer appealing characteristics for the visual exploration of collections of histopathological images and of highdimensional data in general. These techniques compute a projection of the data into a lower-dimensional space while best preserving the information content. The visualization of the low-dimensional space can reveal hidden structures in the data, e. g. the presence of meaningful patterns. A recently proposed algorithm for dimensional data reduction is Locally Linear Embedding (LLE) [1]. It attempts to reduce the dimension of the data while preserving the relationships between neighboring data points. i. e. neighboring data points in the highdimensional space are projected by LLE into neighboring data points in the dimensional reduced space. In other words, the closeness of two data points in the dimensional reduced space reﬂects the similarity of the two respective data points in the high-dimensional space. This property of neighborhood preservation in LLE can be useful for the visual exploration of a collection of histopathological images. Speciﬁcally, the

Visual Exploration of Pathology Images

image space I = N , where N is the number of pixels of an image and each image represents a data point belonging to this space, can be reduced down to two or three dimensions in order to visualize the total information content. Images with similar structural characteristics represent neighboring data points in the image space that are expected to be mapped onto nearby points in the dimensional reduced space. In this way the similarity between several diﬀerent images can be visualized at once in a customized display. On the other hand, the image space is extremely high-dimensional (104 −106 ) while the number of images, i. e. the number of points, is typically limited (102 − 103 ). This means the image space is actually very sparse and it may probe problematic for LLE to compute a faithful projection of such space. Indeed, for this purpose LLE must ﬁnd a global hidden structure in the data and this in turn requires the data to be suﬃciently dense. Ideally, the number of data points should grow exponentially with the dimension of the data space and this fact is known as ”curse of dimensionality”in the machine learning literature [2]. To reduce the sparsity of the image space there are typically two possibilities, either to increase the number of data points (i. e. images) or to reduce the dimension of the image space itself. While the ﬁrst possibility is often not achievable because of practical reasons, e. g. the limited number of patients and images, the second possibility can be done by appropriate techniques of signal processing such as the Discrete Wavelet Transformation (DWT). This is a technique for signal processing used to access the localized and scale-dependent information in signals or images. By computing wavelet based image features we are able to generate a lower-dimensional feature vector that encodes the image information. In this work we present the visual exploration of a dataset of histopathological images of diﬀerent brain tumors based on a combination of DWT and LLE. Each image is preprocessed by DWT and transformed into a feature vector. This DWT pre-processing step is used to reduce the dimension of the image space, thereby reducing its sparsity. These sets of feature vectors represent the LLE input that is reduced down to two dimensions. The dimensionally reduced dataset is ﬁnally visualized with customized colors encoding the type of tumor. This allows the user to explore the similarity among the images of the dataset.

2 2.1

Methods Locally Linear Embedding

Given V D-dimensional vectors {Xi } as input, the LLE algorithm comprises three steps: Step 1: it consists in assigning each data point Xi a predetermined number n of neighbors, typically according to the Euclidean distance. Step 2: by minimizing the following error function Ψ (W ) =

V i=1

|Xi −

n j=1

Wij Xj |2

with

n j=1

Wij = 1

(1)

C. Varini et al.

one computes the weights {Wij } that combined with the neighbors best approximate each data point Xi . Step 3: the weights are used to map each data point Xi into a d-dimensional vector Yi with d < D, such that the following error function is minimized: Φ(Y ) =

|Yi −

i=1

2.2

Wij Yj |2

j=1

with

V V 1 Yi YiT = I , Yi = 0 (2) V i=1 i=1

Discrete Wavelet Transform (DWT)

The DWT analysis enables the assessment of localized and scale-dependent information in signals and images [3]. A signal f is decomposed into a basis of shifted and dilated versions of a mother wavelet ψ [4] f (x) = dj,k ψj,k (x) , with ψj,k (x) = 2−j/2 ψ(2−j x − k) . (3) (j,k)

The index j indicates the dilation or scaling step while k refers to translation or shifting. The wavelet coeﬃcients dj,k are given by the scalar product dj,k = f (x), ψj,k (x) or dj,k = f (x), ψ˜j,k (x) in the case of biorthogonal wavelets with the dual wavelet ψ˜ [4]. An eﬃcient calculation of these coeﬃcients is accomplished by the Fast Wavelet Transform (FWT), an algorithm allowing the coeﬃcients to be calculated in a stepwise manner. To perform a FWT a scaling function φ(x) is required such that [4] √ √ h(k)φ(2x − k) and ψ(x) = 2 g(k)φ(2x − k) . (4) φ(x) = 2 k

The coeﬃcients h(k) and g(k) are termed Filter coeﬃcients. On the ﬁrst scale the signal is decomposed into its details and the remaining signal, i.e. the approximation, reﬂecting the particular scale. The details are described by the wavelet coeﬃcients of this scale while the approximation is represented by scaling coeﬃcients corresponding to the scaling function. The procedure can be iterated by a further decomposition of the approximation into details and approximation of the next coarser scale. In two dimensions the DWT can be applied to each dimension separately, resulting in wavelet coeﬃcients describing the horizontal, the vertical and the diagonal details.

Data and Experiments

The experimental dataset contains 84 histopathological images from 21 diﬀerent cases of meningioma WHO grade I, a benign brain tumor. Each case relates to one of four meningioma subtypes, namely meningotheliomatous, ﬁbroblastic, psammomatous and transitional, the latter showing intermediate features of the characteristics of the former classes. Each pathological image comprises 1300 × 1050 pixels at 400x magniﬁcation. Dividing each image into 16 subimages with 256 × 256 pixels, we derive a data set containing 1344 subimages.

Visual Exploration of Pathology Images

The LLE is performed on two diﬀerent input data: the ﬁrst one is the set of all subimages where each subimage is treated as a point in a 256 × 256dimensional data space; the second one is the set of texture describing features based on DWT. These features are computed as follows. Firstly, an average absolute coeﬃcient is computed for each color channel and each of the ﬁnest six scales: |dj,o (kx , ky )|, j = 1, .., 6, and o = o1 , o2 , o3 . (5) f1 (j) = o kx ,ky

Here j is again the scale index, while index o indicates the orientation in the image. As explained in section 2.2 The indices o1 and o2 indicate coeﬃcients in vertical or horizontal direction, while o3 indicates the diagonal details. Secondly, we introduce a feature set f2 describing a preferred orientation in the image.

f2 (j) =

1 ( |dj,o1 (kx , ky )|− |dj,o2 (kx , ky )| + c |dj,o3 (kx , ky )|).(6) f1 (j) kx ,ky

kx ,ky

with j = 1, .., 6. These computations result in two types of features for each of the ﬁnest six scales and each of the three color channels, leading to a set of 36 features for each subimage. Thus, the input of LLE comprises 1344 36dimensional feature vectors. Both input data are projected into a two-dimensional space and each data point is colored according to its respective histopathological diagnosis. The overlap among the clusters in the two embeddings has been quantiﬁed by computing for each data point the percentage of its nearest neighbors which belong to the same class. The average value computed over the entire dataset quantiﬁes the degree of overlap among the various clusters. The value is comprised between 0 and 1. The closer to 1, the less the clusters overlap.

Results

The LLE projections of the subimages and of the DWL coeﬃcients are shown in Fig. 1. Direct application of LLE to the subimages results in overlapping clusters lacking a clear structure. On the contrary, the DWT preprocessing of the images allows LLE to detect the presence of four diﬀerent clusters. The meningotheliomatous, ﬁbroblastic and psammomatous data points are clearly localized and quite well separated in the projected space. This is in agreement with the histopathological characteristics of the various tissue types. Speciﬁcally, psammomatous tissue is characterized by round calciﬁcations that diﬀer signiﬁcantly from the round structures in meningothelomatous tissue and the elongated structures (cells and cell nuclei) in ﬁbroblastic tissue (see the tissue samples in Fig. 1. At the same time, the transitional subtype cluster largely overlaps with meningotheliomatous and ﬁbroblastic, in agreement with the fact that transitional meningiomas show features incorporating characteristics from the

C. Varini et al.

Fig. 1. Two-dimensional projections obtained by LLE of the histopathological dataset without (left) and with (right) DWT processing.

0.2

LLE

0.1

0.08

0.15

LLE + DWT

0.06 0.1

0.04 0.05

0.02 0

−0.02

−0.05

−0.04 −0.1

transitional fibroblastic psamomatous meningothelomatous

−0.15

−0.2 −0.06

−0.04

−0.02

0.02

0.04

0.06

transitional fibroblastic psamomatous meningothelomatous

−0.06

−0.08

0.08

0.1

0.12

0.14

−0.1 −0.06

−0.04

−0.02

0.02

0.04

0.06

0.08

0.1

0.12

other classes. The matching between the data visualization and the prelabelling of the data ensures that signiﬁcant information has been preserved by the DWT preprocessing approach. The overlap between the clusters has been quantiﬁed as described above and the experimental values are 0.387 for the LLE projection and 0.567 for the LLE projection with the DWT preprocessing.

Discussion

Our results clearly show that LLE beneﬁts from the DWT-based preprocessing. Speciﬁcally, applying LLE directly to the 256 × 256-dimensional space of the subimages proves problematic because of the sparsity of the space itself. By computing wavelet based features, the information of each subimage can be compacted into a lower dimensional feature vector. In this way the dimension of the data is signiﬁcantly reduced and, in turn, the resulting feature space is more dense. This allows LLE to detect clinical relevant data structures in a more eﬃcient manner.

References 1. Roweis ST, Saul LK. Nonlinear Dimensionality Reduction by Locally Linear Embedding. Science 2000;290:2323–2326. 2. Bellmann R. Adaptive Control Processes: A guided Tour. Princeton University Press; 1961. 3. Mallat S, Zhong S. Characterization of Signals from Multiscale Edges. Transactions on Pattern Analysis and Machine Intelligence 1992;14(7):710–732. 4. Daubechies I. Ten Lectures on Wavelets. CBMS-NFS Series Appl. Math., SIAM; 1991.

Strukturprototypen zur Modellierung medizinischer Bildinhalte Benedikt Fischer, Benjamin Winkler, Christian Thies, Mark O. G¨ uld und Thomas M. Lehmann Institut f¨ ur Medizinische Informatik, RWTH Aachen, 52057 Aachen Email: bﬁ[emailprotected]

Zusammenfassung. Medizinische Bildinhalte sind h¨ auﬁg aus mehreren Bildregionen zusammengesetzt, die zueinander in bestimmten Verh¨ altnissen stehen. Diese lassen sich u ¨ ber Wahrscheinlichkeitsverteilungen regionaler und relationaler Attribute zu einem Strukturprototypen verallgemeinern und mit Szenen unbekannten Inhalts vergleichen. Die semiautomatische Synthese eines Strukturprototyps sowie seine Anwendung zur Szenenanalyse wird am Beispiel von Handradiographien vorgestellt. Durch Verwendung der Strukturinformation wird hier eine Steigerung der Erkennungsrate um 17% erzielt.

Einleitung

Aufgrund der bei medizinischem Bildmaterial auftretenden hohen Variabilit¨aten lassen sich nur schwer Normwerte f¨ ur die Attribute zu identiﬁzierender Objekte entwickeln. So sind beispielsweise die Regionenmerkmale von Radiographien der Metakarpalknochen nicht ausreichend trennscharf, um eine korrekte Zuordnung der Regionen zu den einzelnen Fingern zu erlauben. Auch f¨ ur einen menschlichen Betrachter ist es nahezu unm¨ oglich, einen einzelnen Knochen ohne Umgebung korrekt zu klassiﬁzieren. Dies f¨ allt jedoch leicht, wenn die Verh¨altnisse zu anderen Knochen (Position, Gr¨ oße, etc.) bekannt sind. Eine Modellierung einer solchen Szene kann u ¨ ber einen attributierten relationalen Graphen (ARG) erfolgen, der ausschließlich die relevanten Regionen als Knoten und ihre Verh¨ altnisse zueinander als Relationskanten repr¨asentiert. Diese Szene-ARGs kommen daher mit nur wenigen Elementen aus. Bei der Segmentierung von Bildern unbekannten Inhalts m¨ ussen jedoch wegen der unbekannten Aufl¨ osung je Pixel unterschiedliche Regionenzugeh¨origkeiten betrachtet werden. In [1] wird das Bild dazu in einem iterativen Regionenverschmelzungsprozess in einen hierarchischen attributierten Regionenadjazenzgraphen (HARAG) u uhrt. Die resultierenden Graphen sind dadurch jedoch ¨ berf¨ deutlich gr¨ oßer als die Szene-ARGs. Eine Szenenanalyse entspricht dem Vergleich des Szene-ARG mit dem HARAG des zu analysierenden Bildes. Zur allgemeinen Objektsuche anhand von ARGs existiert in der Fachliteratur eine Reihe von L¨osungsans¨atzen [2, 3]. Oft basieren die Referenzen jedoch lediglich auf manuell erzeugten Segmentierungen einer Referenzaufnahme oder auf Kantenbildern, deren Anwendung wegen

B. Fischer et al.

der Delineationsproblematik ausscheidet. H¨ auﬁg k¨onnen mit den vorgestellten L¨osungen auch lediglich kleine Graphen mit weniger als 100 Knoten eﬃzient miteinander verglichen werden oder die hohen Variabilit¨aten in den relevanten Bildstrukturen werden nur unzureichend modelliert. Zudem k¨onnen lediglich Graphen gleichen Typs miteinander verglichen werden.

Methoden

Die Strukturprototypgraphen verallgemeinern die relevanten Bildinhalte f¨ ur einen Anwendungsfall, so dass Variabilit¨ aten nicht nur in Einzelobjekten sondern auch im Zusammenhang ber¨ ucksichtigt werden. Die Repr¨asentation, semiautomatische Synthese und Anwendung des Prototypgraphens wird im Folgenden vorgestellt. Als Anwendungsbeispiel werden Skelettradiographien der linken Hand herangezogen, wie sie bei Kindern zur Bestimmung des Knochenalters verwendet werden. 2.1

Prototyprepr¨ asentation und -training

Der Prototypgraph ist ein vollst¨ andiger Graph, dessen Knoten Einzelobjekte und dessen Kanten Relationen zwischen den Einzelobjekten beschreiben. Im Anwendungsbeispiel entsprechen die Einzelobjekte den einzelnen Handknochen. F¨ ur jedes Einzelobjekt werden die statistischen Verteilungen der regionalen Merkmale als Knotenattribute, die der relationalen Merkmale zwischen den Einzelobjekten als Kantenattribute kodiert. Beispiele f¨ ur ein regionales bzw. relationales Merkmal sind der mittlere Grauwert einer Region bzw. die relative Gr¨oße zweier Regionen. Die Merkmale des Prototypgraphens sind dabei unabh¨angig von den Merkmalen der zur Segmentierung verwendeten HARAGs. Insbesondere werden die relationalen Merkmale erweitert. Die Merkmalsverteilungen werden in einer Trainingsphase unter Annahme von Normalverteilungen gesch¨ atzt: Als Referenz werden zun¨achst die relevanten Knoten f¨ ur jedes Einzelobjekt in automatisch erzeugten HARAGs manuell mit einem entsprechenden Benutzerinterface identiﬁziert. Die bereits w¨ahrend der Segmentierung berechneten regionalen Merkmale werden als Knotenattribute weiterverwendet. Bei Bedarf k¨ onnen auch weitere regionale Merkmale berechnet werden. Danach werden je Datensatz Kanten zwischen den identiﬁzierten Knoten eingef¨ ugt und mit relationalen Merkmalen attributiert. Anschließend kann dann die Verteilung aller Attribute gesch¨ atzt werden. Insgesamt wurden 34 regionale Merkmale zu Form, Intensit¨ at und Textur sowie 12 relationale Merkmale zu Topologie-, Gr¨ oßen- und Intensit¨ atsverh¨ altnissen verwendet. ¨ Mit Hilfe der ermittelten Verteilungen wird dann eine Ahnlichkeitsfunktion ¨ f¨ ur jeden Prototypknoten bzw. jede Prototypkanten deﬁniert, um die Ahnlichkeit unbekannter Regionen bzw. Regionenpaare zu den Einzelobjekten bzw. Einzelobjektpaaren zu bestimmen. Im Anschluss an die Verteilungsberechnung ﬁndet eine Merkmalsselektion statt, um die besten Merkmalsmengen je Region und Relation zu bestimmen und

Strukturprototypen zur Modellierung medizinischer Bildinhalte

Abb. 1. Prozentuales Auftreten der einzelnen Handknochen als eigenst¨ andige Region (links) und Nummerierung der Prototypklassen (rechts). Die Nummerierung besteht aus der Knochenklasse (MK: Metakarpal, UPh: untere, MPh: mittlere, bzw. OPh: obere Phalanx) sowie der Fingernummer (1:Daumen bis 5:kleiner Finger).

83%

84% 85%

85%

82%

88%

58%

OPh-3

80%

OPh-2 MPh-3

72%

67%

86%

76%

44% 67%

70%

73%

MPh-2

61%

UPh-2

OPh-4

MPh-4 OPh-5

UPh-3

OPh-1

UPh-5

UPh-1

36% 42%

MPh-5

UPh-4

MK-2 MK-1

MK-3

MK-4 MK-5

den L¨ osungsraum f¨ ur das Graph-Matching zu verkleinern. Im Anwendungsbeispiel wurde mit Sequential Forward Selection ein ebenso einfaches wie schnelles Verfahren ausgew¨ ahlt [4]. 2.2

Vergleich Prototyp- und Bildgraph

Der Vergleich zwischen Prototyp- und Bildgraph erfolgt auf Basis eines HopﬁeldNetzes [5]. In einem Assoziationsgraphen werden Entsprechungshypothesen zwischen Prototypknoten und HARAG-Knoten als eigener Graphknoten repr¨asen¨ tiert, sofern eine grundlegende Ahnlichkeit der repr¨asentierten Regionen u ¨ ber die regionalen Merkmale vorliegt. Diese dient dann auch gleichzeitig als initiales Potential der Knoten. Das Verfahren wurde modiﬁziert, um die relationalen Merkmale des Prototypen zur St¨ arkung bzw. Hemmung von Hypothesen ebenfalls ber¨ ucksichtigen zu k¨ onnen. F¨ ur die Matching-Hypothesen werden zwischen betrachteten HARAG-Knoten im Bedarfsfall fehlende relationale Merkmale be¨ rechnet. Uber ein Gradientenabstiegsverfahren werden die Potentiale der Matching-Hypothesen solange ver¨ andert, bis die gefundene Zuordnung stabil ist. 2.3

Experimente

Zur Evaluierung wurden 96 Handradiographien verwendet. Sofern in der Segmentierung vorhanden, wurden je Bild bzw. Graph 19 Handknochen (Metakarpalknochen bis obere Phalanx) mit einem Klassenlabel versehen (Abb. 1). Insgesamt wurden 221614 Regionen segmentiert, von denen 1290 den 19 Klassen entsprechen. Abb. 1 illustriert, wie oft die jeweiligen Knochen in der Segmentierung vorhanden waren. Durchschnittlich bestand ein HARAG aus 2300 Knoten.

B. Fischer et al. Tabelle 1. Erkennungsrate der einzelnen Fingerknochen Klasse vorh. lokal strukturell MK-1 70 47 52 64 45 50 MK-2 67 46 53 MK-3 35 24 26 MK-4 40 16 27 MK-5 42 27 29 UPh-1 79 58 59 UPh-2 84 59 63 UPh-3 83 67 71 UPh-4 73 43 60 UPh-5

Klasse vorh. MPh-2 82 MPh-3 82 MPh-4 72 MPh-5 64 OPh-1 56 OPh-2 80 OPh-3 81 OPh-4 77 OPh-5 59

lokal strukturell 45 53 59 68 52 63 38 43 21 32 28 43 44 44 36 47 21 26

Gesamt 1290 776 100% 60.2%

909 70.5%

Zur Evaluierung wurde eine sechsfache Kreuzvalidierung gew¨ahlt, um jeweils einen Prototypen zu berechnen und die Erkennungsleistung zu bestimmen. Eine Regionenzuordnung wird als korrekt gez¨ ahlt, wenn sich die vom Prototypen zugeordnete Region mit der Segmentierung der Ground Truth zu mindestens ¨ 85% u wurde vernachl¨assigt, da ¨ berlappt. Eine hundertprozentige Uberdeckung sich Regionen oft nur in wenigen kleinen Details unterscheiden und auch manuell kaum zu beurteilen ist, welche Region die“ korrekte L¨osung darstellt. So ” werden zum Beispiel bei manchen Regionen die Kapseln in die Segmentierung des Fingers einbezogen, in anderen F¨ allen dagegen nicht.

Ergebnisse

¨ Uber alle Klassen gemittelt wurden ohne relationale Merkmale lediglich 60.2% der Regionen korrekt klassiﬁziert, mit Hilfe des Strukturprototypen konnte dagegen eine Erkennungsrate von 70.5%, d.h. eine Steigerung um u ¨ ber 17%, erzielt werden (Tab. 1). Abbildung 2 gruppiert die Resultate u ¨ ber die jeweilige Knochenklasse. Bei allen Klassen wird eine Verbesserung durch Verwendung der Strukturinformationen erreicht. Das Trainieren der sechs Prototypen beanspruchte durchschnittlich 15 Minuten, das jeweilige Testen der 16 Datens¨atze dagegen lediglich 30 Sekunden.

Diskussion und Ausblick

Die durchgef¨ uhrten Experimente illustrieren die generelle Anwendbarkeit von Strukturprototypen zur Objektsuche. Die Strukturprototypen f¨ uhren in allen Belangen zu einer deutlichen Verbesserung der Erkennungsrate gegen¨ uber Vergleichen ohne strukturelle Information. In Abbildung 2 f¨allt jedoch das schlechte Abschneiden der Klasse OPh auf. Dies ist dadurch zu erkl¨aren, dass die Regionen dieser Klasse sehr klein sind, die Regionengr¨ oße jedoch einen starken Einﬂuss auf

Strukturprototypen zur Modellierung medizinischer Bildinhalte

Abb. 2. Vergleich von lokaler Suche und Struktursuche mit einer geforderten Regionen¨ uberdeckung von 85%. Erkennungsrate nach Knochenklasse

80%

70,47%

60,16%

53,54%

75,67%

78,12%

64,67%

20%

70,36%

40%

75,36%

60% 64,49%

Erkennungsrate

100%

lokal strukturell

0% MK

UPh

MPh

OPh

Gesamt

Knochenklasse

Bildverarbeitung für die Medizin 2006. Algorithmen - Systeme - Anwendungen Proceedings des Workshops vom 19. - 21. März 2006 in Hamburg - PDF Free Download (2024)