Implementieren von technischen Papieralgorithmen in C ++ oder MATLAB

14

Ich bin ein Student der Elektrotechnik. Ich habe viele technische Artikel über Signal- und Bildverarbeitungsalgorithmen (Rekonstruktion, Segmentierung, Filterung usw.) gelesen. Die meisten der in diesen Artikeln gezeigten Algorithmen sind über eine kontinuierliche Zeit und eine kontinuierliche Frequenz definiert und geben häufig die Lösungen in Form von komplizierten Gleichungen an. Wie würden Sie ein technisches Papier von Grund auf in C ++ oder MATLAB implementieren, um die in diesem Papier erzielten Ergebnisse zu replizieren?

Insbesondere habe ich mir die Arbeit "Ein allgemeiner Algorithmus zur Rekonstruktion von Kegelstrahlen" von Wang et al. ( IEEE Trans Med Imaging. 1993; 12 (3): 486-96 ) angesehen und mich gefragt, wie ich überhaupt anfangen soll Implementierung ihres Algorithmus? Gleichung 10 gibt Ihnen die Formel des rekonstruierten Bildes bei. Wie würden Sie das kodieren? Hätten Sie eine for-Schleife, die durch jedes Voxel geht und die entsprechende Formel berechnet? Wie würden Sie Funktionen von Funktionen in dieser Formel codieren? Wie bewerten Sie Funktionen an beliebigen Stellen?

Ich habe das Buch "Digital Image Processing" von Gonzalez und Woods gelesen, bin aber immer noch ratlos. Ich habe auch über die Buchreihe Numerical Recipes gelesen. Wäre das der richtige Weg?

Welche Erfahrungen haben Sie mit der Programmierung von Algorithmen aus Forschungsarbeiten gemacht? Irgendwelche Tipps oder Vorschläge?

Damian
quelle
1
Ich werde mir die Zeitung ansehen, wenn ich die Gelegenheit dazu habe. Aber ich glaube, hier geht es um XYZ-Punkte in einer bestimmten Grafik. Sie definieren einen Scheitelpunkt und arbeiten dann von dort aus.
2
Typischerweise diskretisiert man die Signale durch Abtasten und wandelt dann Integrale in Summen um.
Nibot
Ich habe also etwas über das Abtasten und Umwandeln der Integrale in Summen gelesen. Aber wie bewertet man den Integranden an jedem Abtastpunkt, wenn die Funktionen im Integranden als Matrizen gespeichert sind?
1
Damian, hast du gesehen, wie die Radontransformation durch Rückprojektion invertiert wird? Dies ist ein etwas einfacheres Beispiel, das ich erklären könnte, wenn es Sie interessieren würde. Es wird eher für die Tomographie mit ebenen Wellen als für die in der von Ihnen veröffentlichten Veröffentlichung beschriebene konische Probenahme verwendet. en.wikipedia.org/wiki/Radon_transform
nibot
1
@ mr-crt, kann man stattdessen auf dsp.SE migrieren?
Nibot

Antworten:

15

Signalverarbeitungsalgorithmen, die in kontinuierlicher Zeit / Raum / Frequenz definiert sind, werden typischerweise implementiert, indem das Signal in einem diskreten Gitter abgetastet und Integrale in Summen (und Ableitungen in Differenzen) umgewandelt werden. Raumfilter werden durch Faltung mit einem Faltungskern (dh gewichtete Summe von Nachbarn) implementiert.

Es gibt eine Fülle von Kenntnissen über das Filtern von abgetasteten Zeitdomänensignalen. Zeitbereichsfilter werden entweder als Filter mit endlicher Impulsantwort implementiert , wobei der aktuelle Ausgangsabtastwert als eine gewichtete Summe der vorherigen N Eingangsabtastwerte berechnet wird; oder Filter mit unendlicher Impulsantwort, wobei die aktuelle Ausgabe eine gewichtete Summe der vorherigen Eingaben und vorherigen Ausgaben ist . Formal werden die zeitdiskreten Filter mit der z-Transformation beschrieben , die die zeitdiskrete Analogie zur Laplace-Transformation darstellt . Das bilineare Transformation bildet eine auf die andere ( c2dund d2cin Matlab) ab.

Wie bewerten Sie Funktionen an beliebigen Stellen?

Wenn Sie den Wert eines Signals an einem Punkt benötigen, der nicht direkt auf Ihrem Abtastraster liegt, interpolieren Sie den Wert von nahegelegenen Punkten. Die Interpolation kann so einfach sein, wie die Auswahl der nächstgelegenen Probe, die Berechnung eines gewichteten Durchschnitts der nächstgelegenen Proben oder die Anpassung einer beliebig komplizierten Analysefunktion an die abgetasteten Daten und die Auswertung dieser Funktion an den erforderlichen Koordinaten. Das Interpolieren auf ein gleichmäßig feineres Gitter ist ein Upsampling . Wenn Ihr ursprüngliches (kontinuierliches) Signal keine Details (dh Frequenzen) enthält, die kleiner als die Hälfte des Abtastgitters sind, kann die kontinuierliche Funktion aus der abgetasteten Version (der perfekt rekonstruiert werden Nyquist-Shannon-Abtasttheorem ) . Ein Beispiel, wie Sie in 2D interpolieren können, finden Sie unterbilineare Interpolation .

In Matlab können Sie verwenden interp1oder interp21D oder regelmäßig abgetastet 2D - Daten (jeweils) oder zu interpolieren griddatazu interpolieren aus unregelmäßig abgetasteten 2D - Daten.

Würden Sie eine for-Schleife haben, die durch jedes Voxel geht und die entsprechende Formel berechnet?

Ja genau.

Matlab erspart Ihnen dies durch explizite for-Schleifen, da es für Matrizen und Vektoren (dh mehrdimensionale Arrays) ausgelegt ist. In Matlab wird dies als "Vektorisierung" bezeichnet. Definite Integrale mit angenähert werden sum, cumsum, trapz, cumtrapzetc.

Ich habe das Buch "Digital Image Processing" von Gonzalez und Woods gelesen, bin aber immer noch ratlos. Ich habe auch über die Buchreihe Numerical Recipes gelesen. Wäre das der richtige Weg?

Ja, numerische Rezepte wären ein guter Anfang. Es ist sehr praktisch und deckt die meisten numerischen Methoden ab, die Sie am Ende benötigen. (Sie werden feststellen, dass Matlab bereits alles implementiert, was Sie benötigen, aber numerische Rezepte liefern einen hervorragenden Hintergrund.)

Ich habe eine Klasse "Algorithmen und Datenstrukturen" besucht, sehe aber keinen Zusammenhang zwischen dem dort vorgestellten Material und der Implementierung wissenschaftlicher Algorithmen.

Das in den Kursen "Algorithmen und Datenstrukturen" behandelte Material konzentriert sich in der Regel auf Strukturen wie Listen, Arrays, Bäume und Diagramme, die Ganzzahlen oder Zeichenfolgen enthalten, sowie auf Operationen wie Sortieren und Auswählen: Probleme, für die es normalerweise nur ein einziges korrektes Ergebnis gibt. Wenn es um wissenschaftliche Algorithmen geht, ist dies nur die Hälfte der Geschichte. Die andere Hälfte betrifft Methoden zur Schätzung von reellen Zahlen und analytischen Funktionen. Sie finden dies in einem Kurs über "Numerische Methoden" (oder "Numerische Analyse") wie diesem- Scrollen Sie für die Folien nach unten): wie Sie Sonderfunktionen schätzen, wie Sie Integrale und Ableitungen schätzen usw. Hier besteht eine der Hauptaufgaben darin, die Genauigkeit Ihres Ergebnisses zu schätzen, und ein gängiges Muster besteht darin, eine Routine zu durchlaufen, die die Funktion verbessert schätzen, bis es ausreichend genau ist. (Sie könnten sich fragen, wie Matlab etwas so Einfaches tun kann, als einen Wert von sin(x)für einige zu schätzen x.)


Als einfaches Beispiel sehen Sie hier ein kurzes Skript, das eine Radontransformation eines Bildes in Matlab berechnet. Die Radontransformation nimmt Projektionen eines Bildes über einen Satz von Projektionswinkeln auf. Anstatt zu versuchen, die Projektion entlang eines beliebigen Winkels zu berechnen, drehe ich stattdessen das gesamte Bild mit imrotate, sodass die Projektionsaufnahme immer vertikal ist. Dann können wir die Projektion einfach mit nehmen sum, da diesum einer Matrix einen Vektor zurückgibt, der die Summe über jede Spalte enthält.

Sie können Ihre eigene schreiben, imrotatewenn Sie möchten, mit interp2.

%%# Home-made Radon Tranform

%# load a density map (image).  
A = phantom;

n_pixels = size(A, 1);  %# image width (assume square)

%# At what rotation angles do we want to take projections?
n_thetas = 101;
thetas = linspace(0, 180, n_thetas);

result = zeros(n_thetas, n_pixels);

%# Loop over angles
for ii=1:length(thetas)
    theta = thetas(ii);
    rotated_image = imrotate(A, theta, 'crop');
    result(ii, :) = sum(rotated_image);
end

%# display the result
imagesc(thetas, 1:n_pixels, result.');
xlabel('projection angle [degrees]');

Was früher ein Integral der Dichte entlang eines Strahls war, ist jetzt eine Summe über einer Spalte eines diskret abgetasteten Bildes, was wiederum durch Interpolation des Originalbildes über ein transformiertes Koordinatensystem gefunden wurde.

Nibot
quelle
Wow @nibot, danke für die so ausführliche Antwort. Ich habe eine Klasse "Algorithmen und Datenstrukturen" besucht, sehe aber keinen Zusammenhang zwischen dem dort vorgestellten Material und der Implementierung wissenschaftlicher Algorithmen. Ich lese die Links, die Sie mir gegeben haben, und beginne mit einfacheren Algorithmen (aus Büchern statt Papieren) zu üben.
Damian
Hallo Damian, ich habe meine Antwort bearbeitet, um deinen Kommentar zu adressieren. Ich denke, Sie finden, was Sie in einem Kurs oder Buch über numerische Methoden / numerische Analyse suchen.
Nibot
Während der gesamten Antwort!
Victor Sorokin
@nibot: danke für die bearbeitung. Ich mag den Kurs zur numerischen Analyse, den Sie verlinkt haben. Warum sind die "Filter mit endlicher Impulsantwort" mit der Interpolation verbunden? Ich frage mich, warum dies nicht Teil des Lehrplans als EE-Student ist. Naja. Vielen Dank!
Damian
@Damian: Abtasttheorie, Interpolation / Dezimation, Z-Transformationen, die bilineare Transformation und FIR / IIR-Filter werden in EE-Klassen / Labors unterrichtet, z. B. in Signalen und Systemen, Kommunikationssystemen, linearen Steuerungssystemen und Einführung in DSP. Ich nahm numerische Methoden im Rahmen eines dualen Studiengangs in Computertechnik; Ich denke nicht, dass es von EEs im Allgemeinen verlangt werden sollte.
Eryk Sun
3

Zur hervorragenden Erklärung von nibot noch ein paar weitere Punkte.

  • Numerische Computerumgebungen wie MATLAB, Octave oder SciPy / NumPy ersparen Ihnen eine Menge Mühe, verglichen mit einer generischen Programmiersprache wie C ++. Das Jonglieren mit doubleArrays und Loops ist einfach nicht mit Datentypen wie komplexen Zahlen und Operationen wie Integralen vergleichbar. (Es ist mit Sicherheit machbar und guter C ++ - Code kann eine Größenordnung schneller sein. Mit guten Bibliotheksabstraktionen und -vorlagen kann es sogar einigermaßen sauber und klar sein, aber es ist definitiv einfacher, mit z. B. MATLAB zu beginnen.)

  • MATLAB bietet auch "Toolkits" für z. B. Bildverarbeitung und digitale Signalverarbeitung an , die abhängig von Ihren Aufgaben sehr hilfreich sein können.

  • Die digitale Signalverarbeitung von Mitra ist ein gutes Buch, um (in MATLAB!) Die Grundlagen der diskreten Zeit, Filter, Transformationen usw. zu erlernen.
Joonas Pulakka
quelle
Ja, ich habe die Dokumentation der Image Processing Toolboox gelesen. Ich scheine sehr nützlich zu sein, aber meine Frage war darauf ausgerichtet, so etwas umzusetzen. Grundsätzlich wollte ich wissen, wie man einen mathematischen Algorithmus / eine mathematische Formel verwendet und implementiert (wie es Mathworks mit dem IPT getan hat). Ich wollte etwas über das Gedankenmuster oder einige Richtlinien wissen. Ich werde Mitras Buch anschauen. Vielen Dank!
Damian
1
Um die obige Antwort zu ergänzen, können C ++ - Toolkits wie Armadillo die Konvertierung von Matlab-Code in schnellen C ++ - Code erheblich vereinfachen. Die Syntax von Armadillo ähnelt der von Matlab. Sie können Matlab- und C ++ - Code auch über die mex-Oberfläche von Armadillo mischen.
15.
2

Numerische Methoden. Es ist in der Regel ein Hochschullehrgang und ein Lehrbuch.

DSP befindet sich normalerweise in der Nähe der Schnittstelle zwischen numerischen Methoden und effizienter Implementierung. Wenn Sie die Effizienz ignorieren, suchen Sie möglicherweise nach einer numerischen Approximationsmethode, die ein ausreichend genaues Ergebnis für die interessierenden Gleichungen des technischen Papiers liefert. Manchmal kann es sich um Stichprobendaten handeln, bei denen die Stichprobensätze sowohl die Datenerfassungsmethode (Vorfilterung) als auch den Bereich oder die Qualität der Ergebnisse einschränken, die mit diesen Daten erzielt werden können.

Manchmal verfügen Matlab, numerische Rezepte oder verschiedene Bild- / Signalverarbeitungsbibliotheken über effiziente Algorithmen oder Codes für die gewünschte numerische Lösung. Aber manchmal müssen Sie vielleicht Ihre eigenen rollen, daher ist es hilfreich, die Mathematik hinter den verschiedenen numerischen Lösungsmethoden zu kennen. Und das ist schon ein großes Thema.

hotpaw2
quelle