Algorithmen und Optimierungen

[Tino_Benkenstein]

Hallo Rolf,
erstmal vielen Dank für das tolle Programm und die vorbildliche Dokumentation. Ich habe zwar keine besonderen Programmierkenntnisse in Python, aber ich kenne mich verhältnismäßig gut mit den astronomischen Algorithmen und den technischen Grundlagen aus. Diesbezüglich habe ich bisher haupsäch Makros in ImageJ bzw. AstroImageJ programmiert, mit dem Hintergrund eine vernünftige Datenpipeline für die Auswertung und das Stacking von recht umfänglichen kurzbelichteten Bilddaten einer EM-CCD hauptsächlich bezüglich Deep-Sky zu generieren. Hierbei ist die Qualitätsbestimmung und das Alignment der Einzelbilder problematisch und aufwendig, da das Rauschen sehr dominant und etwas spezieller ist. Das ImageJ hat hierbei den Vorteil, dass man ohne große Programmierkenntnisse Filter und PlugIns per Klick ausprobieren kann und hiermit schon sehr weit kommt.

Es ist aber trotzdem sehr interessant mal in deinem Python-Programm unter die Motorhaube zu schauen, und habe auch gleich mal versucht die Algorithmen und Abläufe nachzuvollziehen und auch schon ein paar Geschwindigkeitsoptimierungen ausprobiert.
Aufgrund der von deinem Programm ausgegebenen Laufzeiten für die Module, ist wohl das eigentliche Qualitäts-Ranking der Bilder der eigentliche Zeitfresser. Als zweites wäre die Berechnung des AP-Shifts zu nennen. Hier ein Beispiel der Ausgabe:

Laplace Stride 2
--------------------------------------------------
Status of time counters:
Execution over all 73.334
Ranking images 33.220
Select optimal alignment patch 0.039
Global frame alignment 4.781
Compute reference frame 0.064
Setting ROI and new reference 0.006
Initialize alignment point object 0.021
Rank frames at alignment points 1.082
Stacking: AP initialization 0.010
Stacking: Initialize background blending 0.032
Stacking: compute AP shifts 16.944
Stacking: remapping and adding 1.528
Stacking: computing background 0.348
Stacking: merging AP buffers 0.010
Stacking: blending APs with background 0.035
Saving the stacked image 0.024
Saving the postprocessed image 0.019
--------------------------------------------------

Soweit ich die Abläufe richtig erfasst habe, passiert beispielweise beim Ranking mittels "bevorzugten" LAPLACE-Operator primär folgendes:
1. Das Video/Bild wird geladen und vorformatiert (FP16 usw.)
2. Im frames.py wird es schonmal vorab mit einem Gauss-Kernel (OpenCV-Funktion im miscellaneous.py) weichgezeichnet
3. Anschließend kommt in rank.frames.py der LAPLACE-Operator (OpenCV-Funktion im miscellaneous.py) zur Anwendung, wobei hier die Gradienten im Bild bestimmt werden (SampleStride skaliert vorher das Bild bspw. auf die Hälfte)
4. Die Standardabeichung des Gradientenbilds ist der Definitionsfaktor für die Bildqualität

Dies ist recht schnell, sortiert aber bei geringer SNR und bestimmten Bildartefakten (HotPixel, Cosmics) diese Problem-Bilder nach vorne. Ansonst funktioniert der Algorithmus ziemlich gut, auch bezüglich Qualität.

Trotzdem habe ich mir nochmal den schon vorgesehenen SOBEL-Filter vorgenommen. Ist ja quasi eine Gradientenbestimmung mit eingebautem Weichzeichner, also alles all-inclusive . Dazu musste ich die Funktion im miscellaneous.py nochmal anpassen, da es ursprünglich nicht richtig funtionierte ("return mag.sum(axis=0)" muss eigentlich durch "return meanStdDev(mag,-1)[1][0][0]" ersetzt werden) und auch im Verhältnis zum Laplace recht langsam war:
Ich habe daher mehrere Anpassungen gemacht:

Original:
frame_int32 = frame[::stride, ::stride].astype('int32')
dx = sobel(frame_int32, 0) # vertical derivative
dy = sobel(frame_int32, 1) # horizontal derivative
mag = hypot(dx, dy) # magnitude
return meanStdDev(mag,-1)[1][0][0]

Optimiert:
dx = Sobel(frame[::stride, ::stride],CV_32F,1,0,ksize=-1) # vertical derivative
# --- Sobel-Funktion nicht aus "from scipy.ndimage import sobel" sondern aus "from cv2 import Sobel" (schneller)
dy = Sobel(frame[::stride, ::stride],CV_32F,0,1,ksize=-1) # horizontal derivative
mag1 = meanStdDev(dx,-1)[1][0][0]
# --- mag = hypot(dx, dy) # hypot-Funktion weggelassen, da sehr zeitaufwendig, meanStdDev der jeweiligen
Richtungsableitungen separat gebildet und abschließend multipliziert, schneller
mag2 = meanStdDev(dy,-1)[1][0][0]
mag = mag1 * mag2
return mag

Zusätzlich kann das künstliche Weichzeichen mittels Gauss im frames.py entfallen. Siehe Laufzeiten:

Ursprünglicher Sobel Stride 2
--------------------------------------------------
Ranking images 40.573
Rank frames at alignment points 18.403
--------------------------------------------------

Optimierter Sobel Stride 2 Ohne GaussianBlur
--------------------------------------------------
Ranking images 26.959
Rank frames at alignment points 7.340
--------------------------------------------------

Zum Vergleich LAPLACE Stride2
--------------------------------------------------
Ranking images 33.220
Rank frames at alignment points 1.082
--------------------------------------------------

Der "Rank frames at alignment points" mit Laplace ist um den Faktor 7 schneller, vermutlich weil das Gauss-weichgezeichnete Bild im
[def frames_mono_blurred_laplacian(self, index)] zwischengespeichert und recycled wird.
Aber jetzt funktioniert der Sobel ertsmal ohne große Abstriche.
Bezüglich Qualität ist übrigens kein Unterschied zu erkennen.

Offene Punkte:
Multiprocessing - es ist meist nur 1 Prozessor in Betrieb

Viele Grüße Tino

 

[Rolf_Hempel]

Hallo Tino,

vielen Dank für den Bericht zu Deinen Experimenten mit dem PSS-Code.

Zunächst möchte ich darauf hinweisen, dass sehr viel Code (insbesondere im Modul "miscellaneous.py") alt ist und in von früheren Tests übriggeblieben ist. Es ist also kein Wunder, wenn diese Funktionen nicht mehr im heutigen Workflow funktionieren.

Die Zeiten, die am Ende ausgegeben werden, sind mit Vorsicht zu genießen. Besonders die scheinbar dominanten Zeiten in "Ranking images" kommen vor allem dadurch zustande, dass an dieser Stelle die Frames erst eingelesen und alle Image-Varianten (grayscale, Gaussian-blurred grayscale) berechnet werden müssen. Je nach Buffering-Level werden diese Daten dann aufgehoben, so dass die späteren Phasen wesentlich schneller ablaufen.

Wenn man die Ausführungszeit optimieren möchte, muss man den ganzen Workflow betrachten. Wie Du richtig bemerkt hast, werden die Daten von "Ranking images" später beim Ranking der AP-Patches wiederverwendet. Ich habe da ziemlich viel Arbeit hineingesteckt, die Gesamtzeit zu optimieren. Insbesondere auch durch Verwendung der zumeist schnelleren OpenCV-Routinen. Ich glaube nicht, dass da noch viel herauszuholen ist.

Übrigens sind die mit dem Gauss-Filter weichgezeichneten Bildversionen nicht nur für das Ranking von Bedeutung, sondern vor allem auch für Verschiebungsberechnungen. Wenn man hier die Weichzeichnung weglässt, landet man unweigerlich in lokalen Optima weit weg von der global optimalen Lösung.

Der einzige Punkt, wo wirklich noch großes Performance-Optimierungpotential liegt, ist das Multiprocessing. Dies ist eine Schwachstelle von Python, das kein brauchbares Parallelisierungskonzept hat. Daher werden Mehrkernprozessoren nicht gut ausgelastet. So ganz stimmt das auch nicht, weil die OpneCV-Bibliotheksroutinen schon parallelisiert sind. Insbesondere bei sehr großen Frames wirkt sich das aus. Aber gerade bei kleinen Framegrößen wäre hier etwas zu holen.

Bisher war mein Entwicklungsziel primär, eine optimale Bildqualität zu erzielen. Eine bis zu einem Faktor 2 längere Ausführungszeit als bei AutoStakkert!3 werden die meisten Benutzer akzeptieren, wenn sie dafür weniger Artefakte in glatten Flächen haben oder etwas mehr Detailzeichnung. Wenn man die Ausführungszeit konsequent optimieren möchte, wird das nur gehen, wenn man zumindest die Berechnungsklassen in C++ umschreibt.

Schöne Grüße
Rolf