© Alexander Kaverin / 123RF.com

Bildbearbeitungen wie Gimp verändern Bilder primär nach künstlerischen Gesichtspunkten. Dagegen bearbeitet die wissenschaftliche Bildverarbeitung mit entsprechenden Algorithmen vor allem in Matrizen gespeicherte Bildpunkte.

Sieht man sich den Aufbau eines Bilds an, fällt auf, dass jede Bildebene aus in Form einer Matrix angeordneter Punkte besteht. Um solche Matrizen zu bearbeiten und umzuformen, eignet sich die Python-Bibliothek NumPy besonders gut. Ein Graustufenbild enthält nur eine Ebene, ein Farbbild drei für die Farbauszüge Rot, Grün und Blau. Eine optionale vierte Ebene, der Alpha-Kanal, steuert die Transparenz der Punkte. Die gängigsten Bildformate sind das JPEG-Format mit drei Farbebenen und das PNG-Format mit einer (Graustufenbild), drei (Farbbild) oder vier Ebenen (Farbe mit Alpha-Kanal).

Dieser Artikel beschränkt sich auf Bilder mit einer Farbtiefe von 8 Bit. Das entspricht 255 Farbwerten je Kanal oder Ebene, die sich entweder als 8-Bit-Integer-Zahlen oder als Dezimalzahlen im Intervall von 0 bis 1 darstellen lassen. Übrigens: Sollte ein Bild nach dem Speichern gänzlich schwarz erscheinen, dann liegt das häufig an einem Konversionsfehler: Weißanteile mit dem dezimalen Maximalwert 1.0 werden dabei als Integer-Zahl 1 interpretiert, was auf der Skala von 0 bis 255 fast reinem Schwarz entspricht.

Den Import und Export von Bildern delegieren die meisten Python-Skripte an die Bibliothek Imageio. Neben den oben genannten Standardformaten kennt sie »Tiff« für Bilder mit beliebiger Farbtiefe, »Hdf5« für komplexe Bildserien, »Dicom« für medizinische Fotodokumentationen, »Fits« aus dem Bereich der Astronomie und diverse Raw-Formate. Sie schreibt und liest außerdem Metadaten aus dem Bild-Header. Die Anweisung »imageio.immeta(Datei.jpg)« sucht beispielsweise nach EXIF-Datenfeldern und schreibt sie in ein Dictionary.

Pythons ebenfalls bekannte Bildbibliothek Pil (kurz für Pillow) stammt aus einer Zeit, als NumPy noch in den Kinderschuhen steckte. Deshalb definiert sie einen Objekttyp für Bilder. Diesen Umweg braucht es heute nicht mehr: Die Bibliothek Imageio steuert die Ein- und Ausgabe von Bildern und erzeugt direkt NumPy-Arrays. Sie liegt in der Version 3 vor und wird als »imageio.v3« importiert. Auf diese Bibliothek greift zudem »matplotlib.image« zurück.

Ein Minimalskript (Listing 1) zur Bilddarstellung mit Speicherung lädt zunächst die Bibliotheken. Die Variable »img« (Zeile 5) speichert die Bilddaten als NumPy-Array, hier in der Dimension (601,800,3). Das Bild mit drei Farbkanälen ist 601 Pixel breit, 800 Pixel hoch und kodiert Farbwerte mit 8-Bit-Integern zwischen 0 und 255. Das Kommando »imshow()« in Zeile 7 gibt das Bild aus (Abbildung 1). Die zwei Speicherbefehle erstellen eine Kopie des Bilds.

Abbildung 1: Der Python-Befehl »imshow()« gibt das Bild schließlich aus.

import numpy as np
import matplotlib.pyplot as plt
import imageio.v3 as iio
fn1 = 'images/A2ost.jpg'
img = iio.imread(fn1)
print(img.shape, img.min().min(), img.max().max())
plt.imshow(img)
plt.savefig('imgs/test1.png')
plt.imsave('imgs/test2.jpg', img)
plt.show()

Der Befehl »img[:,:,0]« wählt die Ebene für den Rot-Kanal aus. Ein realistischeres Graustufenbild erzeugt der Aufruf »img.mean(axis=2)«, indem er über die drei Farbebenen mittelt. Damit das Bild zudem in der Ausgabe grau erscheint, erwartet Matplotlib die Angabe einer passenden Farbtabelle, konkret »cmap = ‘gray’«.

Einfache Array-Operationen wie die aus Listing 2 drehen und spiegeln Bilder, überblenden sie oder setzen sie zu neuen Bildern zusammen (Abbildung 2). Einige Rechenschritte würden bei einer Beschränkung des Wertebereichs auf 8 Bit zu einem Überlauf führen. Mit geschickt eingesetzten Modulo-Rechnungen ließe sich der Grauwertbereich auf das Intervall [0..1] transformieren. Einfacher erledigt das der Aufruf der Funktion »img_norm()« in Zeile 16).

Abbildung 2: Die verschiedenen Array-Operationen aus Listing 2 drehen, spiegeln oder verschieben das Bild.

# 1. rotate by 90 deg
imgp = np.moveaxis(img, 0, 1)
imgp = imgp[:, ::-1,:]
# 2. rotation 180 deg
imgp = img[::-1, ::-1,:]
# 3. rotate by -90 deg
imgp = np.moveaxis(img, 0, 1)
imgp = imgp[::-1, :,:]
# 4. mirror along y-axis
imgp = img[:, ::-1,:]
# 5. reshape
imgp = img[:,::4,:]
imgp = np.repeat(imgp, 2, axis=0)
# 6. fading
imgr = img[::-1, ::-1,:]
imgp = img_norm(img_norm(img)+ img_norm(imgr))
3
# 7. invert
imgp = img_norm(-img)
# 8. cut
imgp = img[::-1, ::1,:] #mirror
imgp = imgp[imgp.shape[0]//2:,:,:]
# 9. join
img9 = img[:img.shape[0]//2,:,:]
img10 = img[::-1, ::1,:]
img10 = img10[img10.shape[0]//2:,:,:]
img11 = np.r_[img9, img10]

Die Matrix-Operationen in Listing 2 basieren auf der Art und Weise, wie NumPy Arrays anspricht. Der Ausdruck »[Anfang:Ende]« legt die Größe eines Arrays fest. Beispielsweise erfasst »[0:17]« die Elemente von 0 bis 16, die Kurzschreibweise »>[:]« dagegen das gesamte Array. Dimensionen trennen Sie durch Kommas. Der Ausdruck »[:,:,0]« liefert die erste Bildebene zurück (bei Farbbildern der Rot-Kanal), »[:,:,2]« den Blau-Auszug.

Neben einem Anfangs- und Endwert kennt NumPy eine Schrittweite. Fehlt die Angabe, nimmt die Bibliothek automatisch eine Schrittweite von 1 an. In der Zeile 12 von Listing 2 vergrößert die Angabe »[::4,::4,:]« die Schrittweite auf 4, nur noch jedes vierte Pixel wird angesprochen. Bei den Elementardrehungen in Listing 2 taucht als Schrittweite »-1« auf. Die Angabe einer negativen Schrittweite kehrt die Reihenfolge um: Das letzte Pixel erscheint zuerst bis zum ursprünglich ersten Pixel, das nun an letzter Stelle steht.

Histogrammoperationen

Der Abbildung 3 liegt der Befehl »plt.hist(img.flatten(), bins=256«) zugrunde. Er sammelt die Helligkeitswerte der Farbkanäle in einem Histogramm mit 256 Abstufungen. Dessen Verlauf legt nahe, die Helligkeit anzupassen und Werte oberhalb von 210 abzuschneiden.

Abbildung 3: Das Histogramm der Farbwerte aller Kanäle des Beispielbilds.

Der Befehl »np.clip()« in der letzten Zeile von Listing 3 beschränkt die Maximalwerte aller drei Farbebenen auf 150. In der Ausgabe bemerkt man davon zunächst nichts, da der Helligkeitsraum weiter im Intervall [0..255]+ beziehungsweise [0..1]+ rangiert. Die Funktion »img_norm()« in der ersten Zeile spreizt die Werte in das Intervall [0..1].

def img_norm(img):
"""
normalize image 0..1
in: img
out: img
"""
# normalize to 0...1
imin = img.min().min().min()
imax = img.max().max().min()
imgo = (img - imin)/(imax-imin)
return imgo
plt.imshow(img_norm(np.clip(img, 0,150))

Eine weitere Option, den Helligkeitsraum besser auszunutzen, bietet eine Gamma-Korrektur. Der Name Gamma steht für den Exponenten, der die Helligkeitswerte gewichtet (Listing 4). Ist er kleiner als eins, hellt er das Bild auf, sonst verdunkelt er es. Die Helligkeitseinstellung eines Monitors basiert auf diesem Prinzip. Setzen Sie sich über den künstlerischen Anspruch der Nachtaufnahme aus Abbildung 4 hinweg, dann fördert eine Gamma-Korrektur mit dem Wert »0.4« (Zeile 4) den zuvor nahezu unsichtbaren Nadelbaum im Vordergrund zutage (Abbildung 5).

Abbildung 4: Eine unbearbeitete Nachtaufnahme – Details verschwinden im Dunkel.

Abbildung 5: Die via Gammakorrektur aufgehellte Nachtaufnahme macht einen Nadelbaum im Vordergrund sichtbar.

gamma = 0.4
img2 = img_norm(img)
img2 = img2 ** gamma
img2 = img_norm(np.clip(img2, 0.2,0.7))
plt.imshow(img2)
plt.imsave('imgs/nightimage_gamma04.jpg', img2)
plt.show()

Ein Histogrammausgleich (Histogram Equalization) passt das Bild selbstständig an, ohne Vorgabe eines Gammawerts (Listing 5). Er versucht, alle Helligkeitswerte gleichmäßig im Bild zu verteilen. Aufgrund des begrenzten Helligkeitsraums gelingt das allerdings nur näherungsweise. Die Funktion »np.histogram()« ordnet die Grauwerte den 256 verfügbaren Helligkeitsstufen zu (Zeile 7). Die kumulierte Summe »cumsum()« wird gestreckt und liefert interpoliert die neuen Grauwerte (Zeile 9).

def img_histequ(img):
""" histogram equalization
in: image
out: equalized image
"""
img = img_norm(img)
histo, bins = np.histogram(img.flatten(),
bins= 256, density=True)
cdf = histo.cumsum()
image_eq = np.interp(img.flatten(), bins[:-1], cdf)
imgeq=image_eq.reshape(img.shape)
imgeq = img_norm(imgeq)
return imgeq

Abbildung 6 vergleicht die kumulierten Grauwerte. In der Ursprungskurve erschöpfen sich die Helligkeitsstufen bereits bei einem Level von 50, die Anpassung nutzt das Gesamtintervall. Listing 5 berechnet den Ausgleich über alle Ebenen. Das Jupyter Notebook aus dem Download-Bereich zu diesem Artikel erweitert das Skript und gleicht die Ebenen unabhängig voneinander aus. Der Blaustich tritt dann stärker hervor als in Abbildung 6.

Abbildung 6: Ein Vergleich der kumulierten Grauwertverteilung vor und nach dem Histogrammausgleich.

Abtastprobleme

Die Bildbibliotheken Scikit-image und Opencv-python erlauben zahlreiche Faltungsoperatoren. Hier geht es nicht um das schlichte Nachprogrammieren: Der Fokus dieses Artikels liegt auf NumPy und Bildoperatoren, die sich mit wenigen Zeilen Code begnügen. NumPy beherrscht lineare Interpolationen, Fourier-Transformationen und Faltungsoperatoren, allerdings ausschließlich für eindimensionale Fälle. Das Nachprogrammieren lohnt nicht, allein schon wegen des Dokumentationsaufwands. Deswegen sei an dieser Stelle auf die beiden erwähnten Bibliotheken verwiesen.

Beim Beurteilen der Wirkungsweise von Bildfiltern helfen Testbilder mit leicht erkennbaren Raumfrequenzen. Listing 6 erzeugt eine passende Funktion (Abbildung 7), aus der eine 2D-Grafik entsteht (Abbildung 8). Der Filter »np.where()« (Zeile 3) wandelt Grauwerte in ein Binärbild um, das nur noch zwei Helligkeitsstufen kennt.

x = np.linspace(0.12,0.45,512)
y = np.sin(30/x)
yr = np.where(y>0.4, 1, 0)
img = np.tile(yr, (yr.shape[0],1))

Abbildung 7: Eine eindimensionale Funktion mit großem Frequenzspektrum.

Abbildung 8: Die Erweiterung der 1D-Funktion zu einer 2D-Grafik.

Die Funktion »np.meshgrid()« aus der dritten Zeile von Listing 7 generiert ein Gitter, hier mit den Abständen zum Bildzentrum. Wenden Sie darauf die 1D-Funktion aus Listing 7 an, gelangen Sie zu einer rotationssymmetrischen Darstellung des Linienbilds (Abbildung 9). Die Grafiken mit ihren variablen Raumfrequenzen unterstützen darin, die Auflösung zu beurteilen, beispielsweise bei einer Druckwiedergabe. Ebenso offenbaren sie die Probleme, wenn Fotos ungeschickt verkleinert werden. Der Aufruf »img2 = imgr[::4,::4]« wählt nur jeden vierten Wert des Originalbilds aus, das Bild schrumpft von 512 auf 128 Pixel Breite.

imsize=512
def mesh_distances(n):
x,y = np.meshgrid(range(n),range(n))
return np.sqrt((x-(n/2)+1)**2+(y-(n/2)+1)**2)
def radialfu(x):
x = x/255*(0.45-0.12) +0.12
return np.sin(30/x) #30/x # moiree: 90/x
d = mesh_distances(imsize)
imgr = radialfu(d)

Abbildung 9: Eine eindimensionale Funktion, rotiert um den Bildmittelpunkt.

Faltungsoperatoren

Ein gleitender Mittelwert (Moving Average) ist ein einfacher Vertreter eines Faltungsoperators. Die Elemente in einem Fenster (auch Kernel genannt) verrechnen die Punkte eines Signals innerhalb dieses Fensters. Beim gleitenden Mittelwert steckt dahinter der Mittelwert der Punkte, der dem Signalpunkt in der Fenstermitte zugeordnet wird. Im nächsten Schritt wandert das Fenster um einen Punkt weiter und erzeugt den nächsten Wert aus der Mittelung innerhalb des Fensters. Zum selben Ergebnis gelangen Sie nach einer Fourier-Transformation und einer Multiplikation mit dem transformierten Kernel.

NumPy kennt weder die 2D-Fourier-Transformation noch 2D-Faltungsoperatoren. Solche Operationen stellt bei Bedarf die Bibliothek Scipy.ndimage zur Verfügung. Zumindest stellt NumPy Methoden für eine eindimensionale Faltung bereit. Damit kann es Kernel berechnen, die sich in zwei 1D-Anteile unterteilen lassen. Dazu gehört beispielsweise der Gauß-Filter.

Listing 8 demonstriert, wie ein 5×5-Kernel aussehen könnte. Der Parameter »sigma« (erste Zeile) steuert die Breite der Verteilung. Der Kernel lässt sich in zwei Komponenten zerlegen, die bei einem Gauß-Filter identisch ausfallen. Das Matrixprodukt aus »[1,6,10,6,1].T @ [1,6,10,6,1]« liefert exakt das Array aus Listing 8. Der Skalierungsfaktor spielt dabei keine Rolle.

imgauss_size, sigma=5, 3
x = np.linspace(-6,6,imgauss_size)
y = np.linspace(-6,6,imgauss_size)
x,y = np.meshgrid(x,y)
imgauss = np.exp(-(x*x+y*y)/(2*sigma**2)) * 10
print(np.linalg.matrix_rank(imgauss))
print(imgauss)
>>> 1
[[ 0 0 1 0 0]
[ 0 3 6 3 0]
[ 1 6 10 6 1]
[ 0 3 6 3 0]
[ 0 0 1 0 0]]

Als doppelte Ableitung einer Funktion verstärkt der Laplace-Filter Bildkanten. Die hier vorgestellten Werkzeuge können ihn nicht nachbilden. Sein Rang beträgt »2«, und er lässt sich nicht separieren (Listing 9, Zeile 12).

sobelx = np.array([1,0,-1])
sobely = np.array([1,2,1])
sobelxy = sobely.reshape(3,-1) @ sobelx.reshape(-1,3)
prewitx = np.array([1, 0,-1])
prewity = np.array([1, 1, 1])
preweitxy = prewity.reshape(3,-1) @ prewitx.reshape(-1,3)
>>>
Laplace
[[ 0 1 0]
[ 1 -4 1]
[ 0 1 0]]
rank: 2
Sobel
[[ 1 0 -1]
[ 2 0 -2]
[ 1 0 -1]]
rank: 1
Prewit
[[ 1 0 -1]
[ 1 0 -1]
[ 1 0 -1]]
rank: 1

Alternativ übernehmen der Sobel- oder der Prewitt-Filter eine ähnliche Aufgabe. Sie zeigen sich rechentechnisch anspruchslos und vom Ergebnis her einem Laplace-Filter ebenbürtig. Der Faltungsfilter aus Listing 10 zusammen mit den x- und y-Sobel-Kerneln aus Listing 9 erzeugen mit dem Aufruf »convol_1d(myimage, sobelx, sobely)« das Ergebnis in Abbildung 10. Im Bild markieren dünne Linien die Straßen, Autos setzen sich aus feinen Strichen zusammen, Grünflächen verwandeln sich in Punktwolken.

Abbildung 10: Das Ergebnis einer Kantendetektion mit dem Sobel-Filter aus Listing 9.

def gaussian1d(x, sigma=1):
return np.exp(-(x*x)/(2*sigma**2))
def convol_1d(img, kernelx, kernely): # , size=9, sigma=0.5):
"""
Convolution filter
in: image, kernelx, kernely
out: img
"""
imgo = np.apply_along_axis(lambda x: np.convolve(x,
kernelx, mode='same'), 0, img)
imgo = np.apply_along_axis(lambda x: np.convolve(x,
kernely, mode='same'), 1, imgo)
# normalize to 0..1
imin = imgo.min().min()
imax = imgo.max().max()
imgo = (imgo

Für die Feineinstellung sorgen ein Clipping und ein anschließendes Normalisieren mit dem Kommando »img_norm(np.clip(imgsob,0.49, 0.51))«. Dabei reduzieren die Parameter »0.49« und »0.51« den Ausgabebereich auf ein enges Intervall um den mittleren Wert von »0.5«. Das Normalisieren spreizt ihn auf das Ausgabeintervall von [0..1]+ (Listing 10, Zeile 13). Ein Speichern im JPEG-Format konvertiert die Werte automatisch auf das Integer-Intervall [0..255].

Filter kombinieren

Der Laplace-Filter operiert selten direkt auf dem Ausgangsbild, da er das Bildrauschen überbetonen würde. Stattdessen glättet zuerst ein Gauß-Filter die hochfrequenten Anteile. Die Kombination beider Filter nennt sich Laplacian-of-Gaussian-Filter (LoG). Er steht hier in NumPy aus den bereits genannten Gründen nicht zur Verfügung.

Ein fast gleichwertiger Ersatz liegt im Difference-of-Gaussian-Filter (DoG). In einem zweistufigen Prozess glätten zwei Gauß-Filter das Ausgangsbild mit unterschiedlichem Glättungsfaktor, beispielsweise »sigma« und »sigma/2«. Der zweite Schritt berechnet die Differenz beider Bilder. Das Ergebnis nach Kontrastanhebung wie beim Sobel-Filter zeigt Abbildung 11. Auch hier treten die Kanten deutlich hervor. Anders als beim Sobel-Filter in Abbildung 10 erscheinen die Linien gröber, der Filter wandelt das Foto in eine Skizze um.

Abbildung 11: Kombiniert ergeben der Laplace- und der Gauß-Filter den Difference-of-Gaussian-Filter.

Nach einem ähnlichen Prinzip arbeitet das unscharfe Maskieren, einer der gebräuchlichsten Schärfungsfilter, der vom Ausgangsbild eine weichgezeichnete Kopie subtrahiert. Das Schärfen eines Bilds erweist sich jedoch als ein zweischneidiges Schwert. Es hebt Kanten hervor, was den Gesamteindruck des Bilds verbessert, da das Auge darauf trainiert ist, Linien zu erkennen. Die Kehrseite der Medaille liegt in der mit dem Schärfen verbundenen Degradation des Bilds. Die entwickelt sich insbesondere dann zum Problem, wenn dem Schärfen weitere Bearbeitungsschritte folgen. Bereits eine Veränderung der Bildgröße, etwa durch Anpassung bei der Druckausgabe, verschlechtert den Bildeindruck.

Der Filter besitzt eine besondere Eigenschaft, sobald er nicht wie beim Bildschärfen moderat auf ein Bild einwirkt, sondern mit hoher Intensität arbeitet. Er verstärkt den Kontrast adaptiv und blendet lokal den Hintergrund aus. Die folgenden Bilder zeichnen die entsprechenden Arbeitsschritte nach. In Abbildung 12 trübt ein Gaußscher Weichzeichner mit dem oben erläuterten Aufruf »convol_1d(imgi, kernel_gauss, kernel_gauss)« das Eingangsbild.

Abbildung 12: Eine weißgezeichnete Version des schlecht leserlichen Demobilds.

Der Beitrag des Differenzbilds aus Abbildung 12 und Abbildung 13 bestimmt die Stärke des Filters. Der Schärfungsgrad legt fest, wie stark das Differenzbild das Ausgangsbild überlagert. In Listing 11 liegt dieser Anteil bei »sharpen=3«. Die zweite Größe »clipping=0.13« definiert, ab welchem Grauwert Pixel weiß erscheinen. Das Ergebnis des übersteuerten unscharfen Maskierens steigert den Kontrast auf Schwarz-Weiß-Niveau und verbessert die Lesbarkeit des Texts erheblich.

Abbildung 13: Ein unscharfes Maskieren fungiert hier letztlich als Kontrastverstärker und macht das Notenblatt erst lesbar.

sharpen=3
clipping=0.13
img_norm(np.clip( imgi + sharping * (imgi-imgb), 0, clipping))

Fazit

In Python steuert die Bibliothek Imageio die Ein- und Ausgabe von Bildern. Sie überführt die Bildpunkte unmittelbar in NumPy-Arrays. In diesem Beitrag haben Sie das Spektrum elementarer Bildbearbeitungsfunktionen kennengelernt, die ausschließlich auf NumPy basieren. Beispielsweise erfolgt ein unscharfes Maskieren direkt aus der Differenzbildung entsprechend vorbearbeiteter Bilder. Dadurch offenbart sich die Logik der Bearbeitungsschritte, die Sie beim Verwenden dezidierter Bildbearbeitungsbibliotheken nicht ohne Weiteres erkennen. (jcb/csi)