Verfasst am: 10.11.2010, 16:09
Titel: Erstellen des Sinogrammes für Iradon
Hallo zusammen
Ich habe gleich 2 Fragen...
Erstens will ich eine Bildrekonstruktion von Röntgenaufnahmen (Tomogram) machen, dafür will ich Iradon verwenden. Jedoch hab ich die Bilder in Form einzeilner Slides, sprich ein Bild von jedem Winkel statt in Form eines Sinogrames -> Wie mach ich aus den einzelnen Projektionen ein sinogram? Oder gibt es hier eine einfachere Möglichkeit ohne über das Sinogram zu gehen?
Zweitens handelt es sich hierbei um 1024x1024 Projektionsbilder und 376 bzw 652 an der Zahl, hierbei hat Matlab nicht so Freude wenn ich sie alle Laden will und in ein 3D Array packe. Ihrgendwelche guten Alternativvorschläge?
Ich verstehe die erste Frage nicht ganz. Was genau meinst du mit "slides"? In einem 2D Datensatz sollte unter einem Winkel das Schwächungsprofil des Objektes vorliegen. Meinst du das mit "slide"? Dann wäre das Sinogramm ja einfach nur eine Matrix, in der in den Spalten die sortierten Schwächungsprofile angeordnet sind.
Wenn du die iradon() Funktion verwenden willst, dann musst du glaube ich schon die Sinogramme als Eingabeparameter übergeben. Alternativ kannst du natürlich die gefilterte Rückprojektion selber schreiben. So kannst du die Implementierung relativ einfach auf deine Eingabewerte anpassen - wobei für die Filterung alle Daten unter einem Winkel vorhanden sein müssen.
Zur zweiten Frage:
Du hast schon recht, das ganze wird ziemlich schnell ziemlich groß. Eine Alternative besteht darin nach der Rekonstruktion eines Bildes die Daten seriell in eine Datei zu schreiben. Falls die Bilder im Anschluss nicht mehr groß bearbeitet, sondern hauptsächlich angeschaut werden sollen, eignet sich ein Programm wie z.B. amide (--> http://amide.sourceforge.net/installation_windows.html)
für das zweite Problem (bilder laden und auslesen)kann ich dir helfen aber wie man aus einzelnen Projektionen ein Sinogram Kriegt weiss ich nicht wenn du die lösung gefunden hast melde dich bitte weil ich habe gleich problem ich mache auch cone beam bild rekonstruktion
Verfasst am: 31.05.2011, 12:58
Titel: zum Beispiel so...
Hallo,
es ist zwar schon länger her, aber vielleicht hilft es doch noch dem ein oder anderen. Gesetzt den Fall du hast 180 Bilder aufgenommen mit einer Kameraauflösung von zum Beispiel 1200 x 1600 pixeln und zwar jeweils eins im Winkel x = 0,1,2,3,4...,179. Um nun die Tomografie machen zu können, schneidest du dein Bild jeweils in Slices. Im einfachsten Fall (aber mit dem größten Speicherplatzbedarf) definierst du jede Pixelspalte deines Bildes als Slice. Nun gehst du zu jedem Bild und schneidest genau diese Spalte heraus und ordnest sie in einer Matrix nacheinander an. Also zB.: du willst den 50 Slice zurück transformieren, dann baust du die Matrix folgender Maßen auf:
erste Spalte RadonMatrix = 50 Spalte Bild 1 aus dem Winkel 0 Grad aufgenommen
zweite Spalte RadonMatrix = 50 Spalte Bild 2 aus dem Winkel 1 Grad aufgenommen
dritte Spalte RadonMatrix = 50 Spalte Bild 3 aus dem Winkel 2 Grad aufgenommen
usw.
so bekommst du nach und nach eine Matrix welche die Größe (AnzahlZeilenBild,180) hat.
Nun brauchst du noch einen Winkelvektor phi. In ihn schreibst du für jede Spalte deines Radonplots, aus welchem Winkel die Projektion aufgenommen wurde. In unserem Beispiel ist das einfach phi = 0:179;
nun kannst du mit der so entstandenen Matrix und dem phi-vector für jede Scheibe iradon aufrufen:
meinSlice = iradon(RadonMatrix,phi);
Falls es sich bei deinen Daten um reale Meßdaten handelt musst du dir unbedingt noch Gedanken darüber machen welchen Filter du verwenden willst.
Zum Speicherplatzproblem:
In vielen Fällen wirst du nicht das ganze Volumen (also alle Slices) im Arbeitsspeicher halten können, ohne das Matlab dir ein 'out of memory' zurück gibt. Dann hast du mehrere Alternativen:
1. Berechne die Slices einzeln und schreibe sie nach jeder Berechnung in einen dafür vorgesehenen Ordner (entweder binär oder einfach nur als Text)
2. Fasse mehrere Slices zusammen und mache dadurch das Gesamtvolumen kleiner
3. Überlege dir genau, was du mit dem Volumen anfangen willst, manchmal reicht es auch aus einzelne Slices oder Teilvolumen zu betrachten, du sparst dir dann Zeit und Raum, wenn du auch nur diese zurücktransformierst
4. Falls du einen ganzen Volumen-Durchlauf machst, speichere die Slices weg und heb sie gut auf, dann brauchst du es nicht noch einmal machen und kannst fortan mit den Slices anstatt auf den Rohdaten arbeiten, außerdem kannst du die berechneten slices dann auch in externe Tomografieanzeigeprogramme die die slices benötigen laden und dort weiter anschauen und analysieren.
Zum Schluss noch ein Beispiel, wie so ein Radonplot bzw. Slice Reader aussehn könnte:
% funktion generiert den entsprechenden Radonplot % pfad := pfad zu dem Ordner in dem sich die aufgenommenen Bilder befinden % musterv := Buchstabenmuster des Dateinamens vor der Winkelnummerierung % musterh := Buchstabenmuster des Dateinamens hinter der Winkelnummerierung % Beispiel: der Dateiname ist : Winkel_1.txt, dann ist musterv = 'Winkel_', und musterh = '.txt' % Generell sollten die Bilder bis auf die Winkelnummerierung alle gleich heißen, damit % man das auslesen automatisieren kann % wSchritt := der Winkelschritt, zum Beispiel 1, wenn du jeweils 1° von 0° bis 179° gegangen bist, du kannst auch in 2,4,6 Grad Schritten gehen und gucken ob das Ergebnis noch befriedigend ist (-> Speicherplatz einsparen), irgendwann wirst du aber auf Artefakte stoßen, du kannst nur Winkelschritte eingeben, die du auch gemessen hast, 0.5 macht also nur Sinn wenn du auch die entsprechenden Bilder dazu hast % wStart := du kannst entscheiden von welchem Winkel dein Radonplot starten soll % wStop := du kannst entscheiden bei welchem Winkel der Radonplot enden soll % sliceNum := welcher Slice gelesen werden soll. Hast du Beispielsweise 1600 Spalten zur Verfügung und gehst jeweils einen Schritt kann der Wert zwischen 1 und 1600 liegen % SliceGr := bestimme, wieviele Spalten zu einem Slice zusammengefasst werden sollen
position = 1;
for iWinkel = wStart:wSchritt:wStopp
datei = strcat(pfad,musterv,num2str(iWinkel),musterh); %Name und Pfad des aktuellen Bildes
rohDaten = load(datei,'-ascii'); % hole das Bild, wenn du es als ascii hast, ansonsten musst du es anders lesen
einSlice = sliceExtract(rohDaten,sliceGr,sliceNum);
radonplot(:,position) = einSlice;
position = position+1;
end
% diese Funktion generiert dir einen Slice, im einfachsten Fall nimmt er einfach die Spalte aus deinem Bild die deinem Slice entspricht, wenn du aber Spalten zu einem Slice zusammen fassen willst, dann mittelt dir diese Funktion diese um einen einzigen Slice-Vektor zu erhalten
% rohDaten := die Daten des aufgenommenen Bildes als Matrix % sliceGr := gib an, wieviele Spalten du zu einem Slice zusammenfassen willst % sliceNum := gib an, welchen Slice du extrahieren möchtest
datenGr = size(rohDaten); % Größe des Bildes bestimmen
ifmod(datenGr(2),sliceGr) == 0 % teste ob deine gewünschte Slice-Größe ohne Rest in die Daten-Größe passt % also wenn dein Bild 1200 lang wäre und du eine Slice Größe von 7 wählen würdest % würden 3 Slices übrig bleiben da 7*171 = 1197 % wenn es passt ist alles ok und hier muss nichts weiter getan werden else % falls etwas übrig bleibt müss das Bild zugeschnitten werden, hier habe ich mich dafür entschieden, die überstehenden Spalten abzuschneiden
schneidWeg = mod(datenGr(2),sliceGr); % Modulo, gibt den Rest an der nicht mehr in die Datengröße passt
rohDaten = rohDaten(1:datenGr(1),1:(datenGr(2)-schneidWeg); % und abschneiden end
% mit den Zeilen, machen wir das selbe, weil wir im letztendlich im 3-dimensionalen sind haben wir keine 2-dimensionalen pixel mehr sondern 3-dimensionale voxel, also noch mal
von der anderen Seite:
% jetzt müsste alles passen und wir können an die vorher gesehene Stelle springen (dieser ganze Aufwand ist wie gesagt nur nötig, weil wir eventuell Slices zusammenfassen wollen, wenn man einfach nur Einzelspalten nimmt, kann man sich das alles auch schenken). Also als nächstes an die entsprechende Stelle springen:
if sliceNum ==1% wir sind vorne und müssen nicht springen spring = 0;
else spring = sliceGr*(sliceNum-1); %springe zu dem Block der Zusammengefasst werden soll
end
% jetzt endlich wird gehüpft und der entsprechende Block herausgenommen:
sliceBlock = rohDaten(:,spring+1:spring+sliceGr);
position = 1;
% nun mitteln über alle voxel im Block um einen einzigen Slice zu bekommen for iSpalte = 1:sliceGr:datenGr(1)% über alle Zeilen des Blocks
voxel = sliceBlock(i:i+sliceGr-1,1:sliceGr);
voxelWert = sum(sum(voxel))/sliceGr^2;
slice(position,1) = voxelWert;
position = position+1;
end
function volumen = getVolumen(pfad,musterv,musterh,wSchritt, wStart, wStop, eSlice,lSlice,sliceSchritt,sliceGr)
% nun kann man mit obigen beiden Funktionen das ganze Volumen berechnen. Die meisten Parameter sind von oben schon bekannt, es sind hinzu gekommen:
% eSlice : = lege fest mit welchem Slice dein Volumen beginnen kann (falls du nur Teilvolumen berechnen willst, ansonsten 1) % lSlice := lege fest mit welchem Slice dein Volumen enden soll % sliceSchritt := lege fest in welchen Abständen die Slices gelesen werden sollen (2 = jeder zweite, 3 = jeder 3. etc.)
phi = wStart:wSchritt:wStop; % der oben angesprochene Winkelvektor
position = 1;
for iSlice = eSlice:sliceSchritt: lSlice % über alle ausgewählten slices disp(['slice' num2str(iSlice)]); % wenn man ungeduldig ist wie ich, will man unbedingt wissen wo matlab gerade rechnet, weil es je nach Größe des Volumens länger dauern kann, so weiß ich immer, ob noch Zeit ist mir einen Kaffee zu machen
slice = getRadonPlot(pfad,musterv,musterh,wSchritt,wStart,wStop,iSlice,sliceGr);
% der slice wird geholt
sliceFertig = iradon(slice, phi,'hamming'); % für meine Zwecke benutze ich am liebsten den Hamming-Filter, aber es kommt auf die Anwendung an
v(:,:,position) = sliceFertig; % slice in das Volumen einfügen end
% falls du das ganze Volumen noch auf 0...1 normieren willst:
% hier kannst du in die eckingen Klammern eintragen welche slices du in x, y und z Richtung sehen willst um die verschiedensten Schnitte durch dein Volumen anzuzeigen, ich habe jetzt eine sehr ungünstige Sichtweise für viele Slices eingestellt, nämlich dass [b]alle[/b] slices in z-Richtung angezeigt werden sollen (eventuell packt dies auch das Matlab nicht mehr), probier ein wenig damit rum um die für dich gerade sinnvollste Anzeige zu finden
So, hoffe ich habe jetzt nicht zu viele Tippfehler gemacht. Viel Spaß beim Tomografieren!
Viele Grüße
Hermine
beamtomo7f_cutlevelxyz.jpg
Beschreibung:
Hier hab ich mal mit obigem Code das Volumen eines Teilchenstrahls in einem Beschleuniger brechnet und mit 6 Schnitten in verschiedenen Richtungen angezeigt. Alle 512 Slices anzuzeigen war da zwar möglich hat aber nicht viel Sinn gemacht. Für die Medizi
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
RSS