Selektive Gauss'sche Unschaerfe - Umsetzung beschleunigen

[etuli]

moin,

ich versuche seit Tagen - hin und wieder - meine Umsetzung fuer den Bildfilter zu beschleunigen. Leider mit nur sehr maessigem Erfolg. Bei einem Bild der Dimension 1600x1200 (24 bit), Radius 10 komme ich auf 16 Sekunden mit meiner besten Umsetzung, was aber voellig inakzeptabel ist! :)

Wer den Algorithmus nicht kennt, hier eine Einfuehrung:

m Gegensatz zu den anderen Weichzeichnungsfiltern wirkt der Selektive Gaußsche Weichzeichner nicht auf alle Pixel des Bildes, der Auswahl oder der aktuellen Ebene. Das Filter wirkt nur auf den Pixel, deren Farbe höchstens um einen definierten Deltawert von der Farbe der Nachbarpixel abweicht. Daher werden Kanten im Bild erhalten.

[...] Radius Hier können Sie den Radius in Pixeln angeben, der zum Berechnen des Filters verwendet wird. Dieser Wert beeinflusst massgeblich die Intensität der Wirkung.
Max. Delta Mit dieser Eigenschaft können Sie die maximale Farbdifferenz im Bereich von 0 bis 255 einstellen. Diese Farbdifferenz bestimmt welche Pixel aus der Umgebung weichgezeichnet werden. Dieser Wert beeinflusst massgeblich wie gut Kanten gegen das Weichzeichnen geschützt werden.





Sei 'src' das Eingabebild, 'dst' das Ausgabebild mit
'bytesPerPixel; bit je Pixel jeweils. Sowohl 'src',
als auch 'dst' haben eine Weite von 'width', wie eine
Hoehe von 'height' Pixeln.

maxDelta ist eine Ganzzahl aus {0,1,...,255}.

numrad ist eine Ganzzahl. Ist r der Radius, so gilt:
numrad = floor (|r| + 2)

pixel (i,x,y,b) liefere das b-te Byte des Pixels
an Stelle (x,y) aus Bild i. i ist gleich 'src', oder
'dst'. Bei einem 24 Bit Bild koennte b=0 rot sein,
b=1 gruen, und b=2 blau. Zum Beispiel. Dies spielt
allerdings keine grosse Rolle.

M sei die Koeffizientenmatrix mit 2*numrad-1 Zeilen
und Spalten. M(i,j) sei der Wert in der i-ten Zeile
und der j-ten Spalte.

Siehe unten fuer ein Beispiel der Matrix und den
zugehoerigen Algorithmus zum Erstellen einer
solchen.

Eine Implementierung des Algorithmus findet sich im
Quellcode von Gimp (plug-ings/common/selgauss.c),
welche mir als Ausgangsbasis diente.

foreach 0 <= y < height
foreach 0 <= x < width
foreach 0 <= b < bytesPerPixel
sum = 0;
fact = 0;

foreach -numrad < i < numrad
if (y + i < 0) or (y + i) >= height then
continue;

foreach -numrad < j < numrad
if (x + j < 0) or (x + j) >= width then
continue;

dist = pixel (src, x,y, b) - pixel (src, x+j,y+i, b);
if (dist > maxDelta) or (dist < -maxDelta) then
continue;

coeff = M(i,j);

sum += coeff * pixel (src, x+j,y+i, b);
fact += coeff;

if (0 == fact) then
pixel (dst, x,y, b) = pixel (src, x,y, b);
else
pixel (dst, x,y, b) = floor (sum / fact);

Der fuehrt in seiner innersten Schleife h*w*b*(2*numrad-1)^2 Schritte aus. Dies entspricht bei meiner obigen Eingabe grob 2.079e9 Iterationen. :D Zuviele!

Leider sehe ich auch keinen Ausweg fuer dieses Problem. Vielleicht ja einer von euch.

Der Algo. fuer die Initialisierung der Matrix enspr. ungefaehr jenem von oben angesprochenem Plugin (insb. bzgl. orig. Kommentaren).

int dx,dy;
double sd, c1, c2;

double ** mat;

mat = new double* [num * 2 - 1];
for (int i = 0 ; i < num * 2 - 1 ; i ++) {
mat [i] = new double [num * 2 - 1];
mat [i] += num - 1;
}

mat += num - 1;

/* This formula isn't really correct, but it'll do */
sd = radius / 3.329042969;
c1 = 1.0 / sqrt (2.0 * 3.1415927 * sd);
c2 = -2.0 * (sd * sd);

for (dy = 0 ; dy < num ; dy++) {
for (dx = dy ; dx < num ; dx++) {
double a = c1 * exp ((dx * dx + dy * dy)/ c2);

mat [ dy][ dx] = a;
mat [ dx][ dy] = a;
mat [-dy][ dx] = a;
mat [-dx][ dy] = a;
mat [-dy][-dx] = a;
mat [-dx][-dy] = a;
mat [ dy][-dx] = a;
mat [ dx][-dy] = a;
}
}

Eine Beipielauspraegung fuer numrad = 7 waere:

0.000 0.000 0.000 0.000 0.001 0.001 0.001 0.001 0.001 0.000 0.000 0.000 0.000
0.000 0.000 0.001 0.002 0.003 0.005 0.006 0.005 0.003 0.002 0.001 0.000 0.000
0.000 0.001 0.002 0.006 0.014 0.022 0.025 0.022 0.014 0.006 0.002 0.001 0.000
0.000 0.002 0.006 0.019 0.040 0.064 0.074 0.064 0.040 0.019 0.006 0.002 0.000
0.001 0.003 0.014 0.040 0.087 0.138 0.161 0.138 0.087 0.040 0.014 0.003 0.001
0.001 0.005 0.022 0.064 0.138 0.218 0.255 0.218 0.138 0.064 0.022 0.005 0.001
0.001 0.006 0.025 0.074 0.161 0.255 0.297 0.255 0.161 0.074 0.025 0.006 0.001
0.001 0.005 0.022 0.064 0.138 0.218 0.255 0.218 0.138 0.064 0.022 0.005 0.001
0.001 0.003 0.014 0.040 0.087 0.138 0.161 0.138 0.087 0.040 0.014 0.003 0.001
0.000 0.002 0.006 0.019 0.040 0.064 0.074 0.064 0.040 0.019 0.006 0.002 0.000
0.000 0.001 0.002 0.006 0.014 0.022 0.025 0.022 0.014 0.006 0.002 0.001 0.000
0.000 0.000 0.001 0.002 0.003 0.005 0.006 0.005 0.003 0.002 0.001 0.000 0.000
0.000 0.000 0.000 0.000 0.001 0.001 0.001 0.001 0.001 0.000 0.000 0.000 0.000


Gruss, Stefan

[badphantom]

Ich bin jetz nich mehr ganz dabei.

heißt das, wenn b = 2.2 habe ich ein satteres blau?

Oder wo liest du denn die Differenz zur Referenzsättigung aus?

[etuli]

Das b ist gleichbedeutend mit 'Kanal'. Der Algo arbeitet aso fuer jeden die Unscharfe durch.

[DarkTom]

Ich habe mal ein wenig darüber gegrübelt. Das Problem sehe ich in der Matrix, deren Werte sich leider nicht für eine Optimierung anbieten. Wenn du da z.B. von der Mitte ausgehend Richtung Rand immer eine Halbierung des Wertes hättest, könnte man da mehr machen. Aber das würde sich halt schon stark auf das Ergebnis auswirken.

Anbieten könnte ich dir eine Verbesserung für den Fall, wenn nur ein Teil der umgebenden Punkte relevant ist. Dazu benötigst du einen Array der Grösse 256 in dem jedes Element eine Liste von Punkten darstellt, die den zugehörigen Farbwert haben. Wenn du einen Pixel weiter nach rechts gehst, packst du die Pixel aus der Spalte, die jetzt neu zu deiner radius-Umgebung gehört, vorne in die jeweiligen Listen. Dann betrachtest du den Wert w des aktuellen Pixels und guckst dir dann die Pixel in den Listen w-delta bis w+delta an. Erreichst du in einer Liste einen Pixel , der zu weit entfernt ist, kannst du da auch mit der Liste aufhören.

Das war jetzt nur kurz angerissen, aber ich glaube, du verstehst das so.
Leider bringt das halt keinen Vorteil, wenn du viele einfarbige Flächen in deinem Bild hast.

[etuli]

Danke erstmal fuer eure Antworten. Deinen Ansatz, Darktom, werde ich mir nochmal genauer anschauen. Das zeitkritsche Problem wird wohl sein, genau diejenigen Pixel aus den Listen zu schmeissen, welche aus der r-Umgebung fallen.

Vielleicht findet sich ja auch noch ein Naeherungsalgorithmus mit einer gescheiten Laufzeit in der Anwendung. :)

Gruss, Stefan

[fraal]

Traurig spreche ich nicht viel Deutsches. Ich benutze einen on-line-Übersetzer. Ich schrieb diesen Algorithmus unter Verwendung IDL neu. Meine neue Version lässt ungefähr 30% laufen, das schneller als die Vorlage ist. Es ist viel schwieriger zu verstehen, aber den gleichen Ausgang schneller zu erzeugen ein wenig. Merken Sie, dass die Vermehrung der Matrizen NICHT zutreffendes Matrix multiplcation ist - es ist einfache Skalarvermehrung der Reihen der selben Größe. Hoffen Sie, dass dieses hilft, aber es vermutlich wird, nicht weil es in IDL geschrieben wird! Alex.

(Anmerkung - wahlweise freigestellter sobel Randabfragungscode auch eingeschlossen)


;+
; NAME: MySelGauss.pro
;
;
; PURPOSE: Selectively Gaussian-blurs an image, with a user definable "delta" blur threshold, and a user-definable blur radius
; This algorithm will actually blur 'along' an edge, but not 'across' it!
; Optional Sobel edge detection.
;
; AUTHOR: Alex Fraser, ACE CRC / IASOS / University of Tasmania
;
;
; CATEGORY: Image processing FUNCTION
;
;
; CALLING SEQUENCE: newfile = MySelGauss('image.img', x=1024, y=768, radius=5, delta=1, /sobel)
;
;
; INPUTS: file: The filename of the image you want to process.
; Must be an unformatted (uncompressed) binary (like an ENVI .img) of size x*y
; Can be a range of data types - just change the type of the variable you are reading the file into to match (eg, float, double)
; x: the x-dimension of the image
; y: the y-dimension of the image
; radius: the (integer, please!) "radius" of the Gaussian function used for convolution
; delta: the maximum value that the difference between two pixels can take for blurring to take place
;
; KEYWORDS:
; /sobel: will sobel-filter the output image before returning
;
;
; OUTPUTS:
; function returns a new data variable (2D array) of the same time as the input data
;
;
; RESTRICTIONS:
; heaps
;
;
; EXAMPLE:
; see calling sequence
;
;
; MODIFICATION HISTORY:
;
; Written by Alex Fraser, 2/6/2008.
;-

FUNCTION MySelGauss, input, x=xdim, y=ydim, radius=rad, delta=del, sobel=sob

;===========================================================================================================
; Generate the Gaussian convolution array. Will be (square, and) of size (2*radius)-1, so as to have the convolution array 'centered' around the px of interest
xs = FIndGen((2*rad)-1)+1-rad
ys = FIndGen((2*rad)-1)+1-rad

GaussArray = FltArr((2*rad)-1, (2*rad)-1)

; set up the Gaussian forumla constants
sigma = rad/3.32904 ;this line is a rough way of getting "rad" to be around what we think it should be (the forumla uses sigma, not rad. Of
;course, the Gaussian function never drops to exactly zero.)
frontConstant = 1/(2.0 * 3.14159 * sigma * sigma)
expconst = -2.0*sigma*sigma

FOR i=0, (2*rad)-2 DO BEGIN
FOR j=0, (2*rad)-2 DO BEGIN
GaussArray[i,j]=frontConstant*exp(((xs[i]*xs[i])+(ys[j]*ys[j]))/expconst)
ENDFOR
ENDFOR
;===========================================================================================================



;===========================================================================================================
;OK, now read the file of interest
if (n_elements(input) eq 0) then begin
input = dialog_pickfile(TITLE="Choose the image file to open")
endif

image = FltArr(xdim,ydim)

;read it in (assume unformatted 32-bit float of size xdim*ydim)
openr, lun, input, /get_lun
ReadU, lun, image
free_lun, lun
;===========================================================================================================



;===========================================================================================================
; Perform the smoothing.

;generate a new output image
outimg = FltArr(xdim,ydim)

;loop through every pixel
for y=0, ydim-1 DO BEGIN
print, y
for x=0, xdim-1 DO BEGIN

;I've written 2 algorithms. The algorithm with the longer code is about 30% quicker. I've left it uncommented-out.

;generate a blank blur binary matrix
blurmatrix = BytArr((2*rad)-1, (2*rad)-1)

;now loop through the kernel
for i=1-rad, rad-1 DO BEGIN
for j=1-rad, rad-1 DO BEGIN
; if any coordinate is out of bounds
if (((y+i) lt 0) or ((y+i) ge ydim-1) or ((x+j) lt 0) or ((x+j) ge xdim-1)) THEN begin
;set its blur matrix entry to zero
blurmatrix[j+rad-1,i+rad-1] = 0
endif else begin
;otherwise, perform the difference test
diff = abs(image[x,y] - image[x+j, y+i])

;if their difference is less than delta, mark it for smoothing by putting a 1 in the blur matrix
if (diff le del) THEN BEGIN
blurmatrix[j+rad-1, i+rad-1]=1
ENDIF
endelse
endfor
endfor

; now, construct the "image data" array to be the same size as the kernel/blurmatrix. Pad at the edges (with any value)
imagedata = FltArr((2*rad)-1, (2*rad)-1)

;normal case of no over or underflows:
xoldindexlow = x+1-rad
xoldindexhigh = x-1+rad
xnewindexlow = 0
xnewindexhigh = 2*rad - 2
yoldindexlow = y+1-rad
yoldindexhigh = y-1+rad
ynewindexlow = 0
ynewindexhigh = 2*rad - 2


;detect under and overflows in both directions, calculate their magnitudes

xunder = 0
yunder = 0
xover = 0
yover = 0

if x-rad lt -1 then begin
;print, 'x underflow found'
xunder = abs(x+1-rad)
xoldindexlow = 0
xoldindexhigh = 2*rad-(2+xunder)
xnewindexlow = xunder
xnewindexhigh = 2*rad - 2
endif
if y-rad lt 0 then begin
;print, 'y underflow found'
yunder = abs(y+1-rad)
yoldindexlow = 0
yoldindexhigh = 2*rad-(2+yunder)
ynewindexlow = yunder
ynewindexhigh = 2*rad - 2
endif
if x+rad gt xdim-1 then begin
;print, 'x overflow found'
xover = abs(x+rad-xdim)
xoldindexlow = x+1-rad
xoldindexhigh = xdim-1
xnewindexlow = 0
xnewindexhigh = 2*rad - (2 + xover)
endif
if y+rad gt ydim-1 then begin
;print, 'y overflow found'
yover = abs(y+rad-ydim)
yoldindexlow = y+1-rad
yoldindexhigh = ydim-1
ynewindexlow = 0
ynewindexhigh = 2*rad - (2 + yover)
endif

;print, 'xoldindexlow '+string(xoldindexlow)
;print, 'xoldindexhigh '+string(xoldindexhigh)
;print, 'xnewindexlow '+string(xnewindexlow)
;print, 'xnewindexhigh '+string(xnewindexhigh)
;print, 'yoldindexlow '+string(yoldindexlow)
;print, 'yoldindexhigh '+string(yoldindexhigh)
;print, 'ynewindexlow '+string(ynewindexlow)
;print, 'ynewindexhigh '+string(ynewindexhigh)
;print, '========================='

;general formula, here we come. Note that you can't have an overflow and an underflow in the same direction at the same time.
imagedata[xnewindexlow:xnewindexhigh, ynewindexlow:ynewindexhigh] = image[xoldindexlow:xoldindexhigh, yoldindexlow:yoldindexhigh]

;Now, multiply then sum the blur and gaussian arrays to get the factor, and all 3 arrays to get the total
intermed = blurmatrix * GaussArray
factor = total(intermed)
total = total(intermed * imagedata)

outimg[x,y] = total/factor


;THIS IS THE ALTERNATIVE ALGORITHM (it's about 30% slower)
; total=0.0
; factor=0.0
;
; ;now loop through the kernel
; for i=1-rad, rad-1 DO BEGIN
; for j=1-rad, rad-1 DO BEGIN
; ;check that we're within the limits of the img array for this combo of x+j, y+i
; if (~ (((y+i) lt 0) or ((y+i) ge ydim-1) or ((x+j) lt 0) or ((x+j) ge xdim-1) ) ) THEN begin
;
; ;all ok, so continue. Calculate the difference between the pixel values
; diff = abs(image[x,y] - image[x+j, y+i])
;
; ;if their difference is less than delta, smooth.
; if (diff le del) THEN BEGIN
;
;
; coeff = GaussArray[i+rad-1,j+rad-1]
; total = total + (coeff*image[x+j, y+i])
; factor = factor + coeff
; ENDIF
; ENDIF
; ENDFOR
; ENDFOR
;
; if (factor eq 0.0) then begin
; ;don't smooth
; outimg[x,y] = image[x,y]
; ENDIF ELSE begin
; outimg[x,y] = total/factor
; endelse

ENDFOR
ENDFOR
;===========================================================================================================



;===========================================================================================================
; Sobel edge detect if the keyword is set
if (keyword_set(sob)) THEN BEGIN
; generate the sobel convolution kernels (one in x, one in y)
sobelx = [[-1,0,1], [-2,0,2], [-1,0,1]]
sobely = [[1,2,1], [0,0,0], [-1,-2,-1]]

;perform the convolution
outputx = convol(outimg, sobelx, center=1)
outputy = convol(outimg, sobely, center=1)

;combine the convolutions
outimg = (outputx^2 + outputy^2)^0.5

ENDIF

;===========================================================================================================

return, outimg
END