Gaussian Blur

[zawran]

Do anyone have a working example, basic language variant if possible, of a gaussian blur function?

[stormbringer]

if you know how to make a simple box blur (simple average of pixels in a region), then apply the algorthim 3 times, each pass reuse the previous result => fast approximation to gaussian blur. The more passes you use, the closer you get to a perfect gaussian blur...

[va!n]

@zawran:
possible one of this nice sites will help you: 

http://www.catenary.com/howto/gaussian.html
http://www.jhlabs.com/ip/blurring.html
http://canopus.physik.uni-potsdam.de/~axm/bgsmooth.html
http://en.wikipedia.org/wiki/Gaussian_blur
http://www.gamasutra.com/features/20010209/evans_01.htm
http://www.pages.drexel.edu/~nk752/cannyTut2.html
http://www.gamedev.net/community/forums/topic.asp?topic_id=427696
http://ciri.be/blog/?p=14

[zawran]

I found some blitzbasic code, which I have converted to BlitzMax. It seems to be working like it is supposed to, so I will use this for now.

Code Select Expand


' pixmap gaussian blur function
' based on blitzbasic code by Elias_T
' maskWidth & maskHeight should always be uneven numbers [3,4,7,9]
Function gaussianBlur:TPixmap(pMap:TPixmap,maskWidth:Int,maskHeight:Int)

Local tmp:TPixmap = CopyPixmap(pMap)
Local width:Int = pMap.width
Local height:Int = pMap.height
Local texel:Float[width,height,3]
Local result:Float[width,height,3]
Local maskData:Float[width*height]

Local x:Int,y:Int,ym:Int,xm:Int
Local cy:Float,cx:Float,rt:Float
Local r1:Float,g1:Float,b1:Float
Local rr:Float,gg:Float,bb:Float
Local mult:Float = 0.0

For x = 0 To width-1
For y = 0 To height-1
rgb = ReadPixel(pMap,x,y)
texel[x,y,0] = (rgb Shr 16) & 255
texel[x,y,1] = (rgb Shr 8) & 255
texel[x,y,2] = rgb & 255
Next
Next

For ym = 0 To maskHeight-1
For xm = 0 To maskWidth-1
cx = Float(xm - (maskWidth - 1) / 2.0)
cy = Float(ym - (maskHeight - 1) / 2.0)
rt = cx*cx + cy*cy
mult :+ Exp(-0.35 * rt)
Next
Next

mult = 1.0 / mult

For ym = 0 To maskHeight-1
For xm = 0 To maskWidth-1
cx = xm - (maskWidth-1) / 2.0
cy = ym - (maskHeight-1) / 2.0
rt = cx*cx + cy*cy
maskData[ym * maskWidth + xm] = mult * Exp(-0.35 * rt)
Next
Next

For ym = 0 To height-1
For xm = 0 To width-1
rr = 0.0
gg = 0.0
bb = 0.0
For yy = 0 To maskHeight-1
For xx = 0 To maskWidth-1
If (xm+xx-Floor(maskWidth/2.0)<0) Or (ym+yy-Floor(maskHeight/2.0)<0) Or (xm+xx-Floor(maskWidth/2.0)>width-1) Or (ym+yy-Floor(maskHeight/2.0)>height-1) Then
r1 = 0.0
g1 = 0.0
b1 = 0.0
Else
r1 = texel[xm + xx - Floor(maskWidth/2.0), ym + yy - Floor(maskHeight/2.0),0]
g1 = texel[xm + xx - Floor(maskWidth/2.0), ym + yy - Floor(maskHeight/2.0),1]
b1 = texel[xm + xx - Floor(maskWidth/2.0), ym + yy - Floor(maskHeight/2.0),2]
End If
rr :+ r1 * maskData[xx + yy * maskWidth]
gg :+ g1 * maskData[xx + yy * maskWidth]
bb :+ b1 * maskData[xx + yy * maskWidth]
Next
Next
result[xm,ym,0] = rr
result[xm,ym,1] = gg
result[xm,ym,2] = bb
Next
Next

For x = 0 To width-1
For y = 0 To height-1
WritePixel(tmp,x,y,Int(result[x,y,0]) Shl 16 + Int(result[x,y,1]) Shl 8 + Int(result[x,y,2]))
Next
Next
Return tmp
End Function


Its not the fastest thing around, and I should probably try and convert it into using pointers instead of read-/write-pixels. And it just might be faster to do a regular box blur a couple of times with almost the same end result. But for now, this will do fine for what I needed it for, where it isn't too critical if its fast or not, at least not at this point.

[stormbringer]

first thing to optimize in your code is the reading of pixels. Take advantage of the fact that a gaussian blur is separable, which means you can process lines and columns separately. I won't go in the math details here why it is separable but it is, so do the following:

1) blur the lines first (blur on x axis)
2) take the result and feed it into the next blur for the columns (blur on y axis)

this will reduce the pixel access by a great factor and definitely improve the speed of the algorithm.

[p01]

Which reduces the bluring to one add, one sub and one div/lookup per pixel and per pass.

[relsoft]

http://blackpawn.com/texts/blur/default.html

fastest one I've ever used. :*)

[Shockwave]

Hell yes!

Nice one Rel, that's very useful. I came home from work to find this forum full of lovely useful replies today

[p01]

relsoft: this is exactly what I meant by one add + one sub + one div/lookUp per pixel per pass.

[zawran]

@relsoft, thanks for the link, I will take a look.

[zawran]

I did one small change to the blur code, which knocked a lot of the processing time of it. Don't know, why I didn't do that in the first place:

Code Select Expand


Graphics 1024,768
' test image = 256x256 pixels
' first version = 4956ms
' second version = 108ms
Cls
Local pixm:TPixmap = LoadPixmap("testimage.png")
Local time:Int = MilliSecs()
Local pixm2:TPixmap = gaussianBlur(pixm,9,9)
Print MilliSecs()-time
DrawPixmap(pixm,0,0)
DrawPixmap(pixm2,256,0)
Flip
WaitKey()
End

' pixmap gaussian blur function
' based on blitzbasic code by Elias_T
' maskWidth & maskHeight should always be uneven numbers [3,5,7,9]
Function gaussianBlur:TPixmap(pMap:TPixmap,maskWidth:Int,maskHeight:Int)

Local tmp:TPixmap = CopyPixmap(pMap)
Local width:Int = pMap.width
Local height:Int = pMap.height
Local texel:Float[width,height,3]
Local result:Float[width,height,3]
Local maskData:Float[width*height]

Local x:Int,y:Int,ym:Int,xm:Int
Local cy:Float,cx:Float,rt:Float
Local r1:Float,g1:Float,b1:Float
Local rr:Float,gg:Float,bb:Float
Local mult:Float = 0.0

For x = 0 To width-1
For y = 0 To height-1
rgb = ReadPixel(pMap,x,y)
texel[x,y,0] = (rgb Shr 16) & 255
texel[x,y,1] = (rgb Shr 8) & 255
texel[x,y,2] = rgb & 255
Next
Next

For ym = 0 To maskHeight-1
For xm = 0 To maskWidth-1
cx = Float(xm - (maskWidth - 1) / 2.0)
cy = Float(ym - (maskHeight - 1) / 2.0)
rt = cx*cx + cy*cy
mult :+ Exp(-0.35 * rt)
Next
Next

mult = 1.0 / mult

For ym = 0 To maskHeight-1
For xm = 0 To maskWidth-1
cx = xm - (maskWidth-1) / 2.0
cy = ym - (maskHeight-1) / 2.0
rt = cx*cx + cy*cy
maskData[ym * maskWidth + xm] = mult * Exp(-0.35 * rt)
Next
Next

Local mhhalf:Int = Floor(maskHeight/2.0)
Local mwhalf:Int = Floor(maskWidth/2.0)
For ym = 0 To height-1
For xm = 0 To width-1
rr = 0.0
gg = 0.0
bb = 0.0
For yy = 0 To maskHeight-1
For xx = 0 To maskWidth-1
If (xm+xx-mwhalf<0) Or (ym+yy-mhhalf<0) Or (xm+xx-mwhalf>width-1) Or (ym+yy-mhhalf>height-1) Then
r1 = 0.0
g1 = 0.0
b1 = 0.0
Else
r1 = texel[xm + xx - mwhalf, ym + yy - mhhalf,0]
g1 = texel[xm + xx - mwhalf, ym + yy - mhhalf,1]
b1 = texel[xm + xx - mwhalf, ym + yy - mhhalf,2]
End If
rr :+ r1 * maskData[xx + yy * maskWidth]
gg :+ g1 * maskData[xx + yy * maskWidth]
bb :+ b1 * maskData[xx + yy * maskWidth]
Next
Next
result[xm,ym,0] = rr
result[xm,ym,1] = gg
result[xm,ym,2] = bb
Next
Next

For x = 0 To width-1
For y = 0 To height-1
WritePixel(tmp,x,y,Int(result[x,y,0]) Shl 16 + Int(result[x,y,1]) Shl 8 + Int(result[x,y,2]))
Next
Next
Return tmp
End Function


Difference is that I calculate that / 2.0 value before the loop instead of each time within the loop. And it went from almost 5000ms to about 110ms, so thats a big saving just for that.

[Stonemonkey]

since maskwidth and maskheight are ints you could use shift right instead of divide, also since xm and ym are ints you could make cx and cy ints and convert to float in the line:

rt = float(cx*cx + cy*cy)

so you only have one int to float in each of the loops.

Another thing I would do is something like this:

Code Select Expand


y_index=ym-mhhalf

For yy = 0 To maskHeight-1

If (y_index>=0) and (y_index<height) Then

mask_index=yy * maskwidth
x_index=xm-mwhalf

For xx = 0 To maskWidth-1

If (x_index>=0) and (x_index<width)  Then

rr :+ texel[x_index, y_index,0] * maskData[mask_index]
gg :+ texel[x_index, y_index,1] * maskData[mask_index]
bb :+ texel[x_index, y_index,2] * maskData[mask_index]

End If

mask_index :+ 1
x_index :+ 1

Next

End If

y_index :+ 1

Next


to take as much of the calculations as possible out of the loops (although no idea if i have that all right as i just edited it in notepad), it should even be possible to remove the If from the inner loop completely.

[zawran]

By making the changes you suggested, it came down another 20ms.

Code Select Expand


Graphics 1024,768
' test image = 256x256 pixels
' first version = 4956ms
' second version = 108ms
' third version = 82ms

Cls
Local pixm:TPixmap = LoadPixmap("testimage.png")
Local time:Int = MilliSecs()
Local pixm2:TPixmap = gaussianBlur(pixm,9,9)
Print "processtime was: "+(MilliSecs()-time)+"ms"
DrawPixmap(pixm,0,0)
DrawPixmap(pixm2,256,0)
Flip
WaitKey()
End

' pixmap gaussian blur function
' based on blitzbasic code by Elias_T
' maskWidth & maskHeight should always be uneven numbers [3,5,7,9]
Function gaussianBlur:TPixmap(pMap:TPixmap,maskWidth:Int,maskHeight:Int)

Local tmp:TPixmap = CopyPixmap(pMap)
Local width:Int = pMap.width
Local height:Int = pMap.height
Local texel:Float[width,height,3]
Local result:Float[width,height,3]
Local maskData:Float[width*height]

Local x:Int,y:Int,ym:Int,xm:Int
Local cy:Int,cx:Int,rt:Float
Local r1:Float,g1:Float,b1:Float
Local rr:Float,gg:Float,bb:Float
Local mult:Float = 0.0
Local mwhalf:Int,mhhalf:Int
Local y_index:Int,x_index:Int,mask_index:Int

mwhalf = Floor((maskWidth-1) / 2)
mhhalf = Floor((maskHeight-1) / 2)

For x = 0 To width-1
For y = 0 To height-1
rgb = ReadPixel(pMap,x,y)
texel[x,y,0] = (rgb Shr 16) & 255
texel[x,y,1] = (rgb Shr 8) & 255
texel[x,y,2] = rgb & 255
Next
Next

For ym = 0 To maskHeight-1
For xm = 0 To maskWidth-1
cx = xm - mwhalf
cy = ym - mhhalf
rt = Float(cx*cx + cy*cy)
mult :+ Exp(-0.35 * rt)
Next
Next

mult = 1.0 / mult

For ym = 0 To maskHeight-1
For xm = 0 To maskWidth-1
cx = xm - mwhalf
cy = ym - mhhalf
rt = Float(cx*cx + cy*cy)
maskData[ym * maskWidth + xm] = mult * Exp(-0.35 * rt)
Next
Next

mhhalf = Floor(maskHeight/2.0)
mwhalf = Floor(maskWidth/2.0)
For ym = 0 To height-1
For xm = 0 To width-1
rr = 0.0
gg = 0.0
bb = 0.0
y_index = ym-mhhalf
For yy = 0 To maskHeight-1
If y_index => 0 And y_index < height Then
mask_index = yy * maskwidth
x_index = xm-mwhalf
For xx = 0 To maskWidth-1
If x_index => 0 And x_index < width Then
rr :+ texel[x_index, y_index,0] * maskData[mask_index]
gg :+ texel[x_index, y_index,1] * maskData[mask_index]
bb :+ texel[x_index, y_index,2] * maskData[mask_index]
End If
mask_index :+ 1
x_index :+ 1
Next
End If
y_index :+ 1
Next
result[xm,ym,0] = rr
result[xm,ym,1] = gg
result[xm,ym,2] = bb
Next
Next

For x = 0 To width-1
For y = 0 To height-1
WritePixel(tmp,x,y,Int(result[x,y,0]) Shl 16 + Int(result[x,y,1]) Shl 8 + Int(result[x,y,2]))
Next
Next
Return tmp
End Function


I tried removing the if/then check in the inner loop, but it crashes, so its still needed. But at least it took off a good chunk of the process time, about 24%, which isn't half bad. Thanks for the tips.

[Stonemonkey]

Sorry, yeah it'd need a few other changes before you can remove the If but good to hear the other changes made some difference.

[zawran]

I went ahead and changed it from read/write pixels to pointers, and shave another 8ms of the process time, not a huge difference, but still

Code Select Expand


Graphics 1024,768
' test image = 256x256 pixels
' first version = 4956ms
' second version = 108ms
' third version = 82ms
' forth version = 74ms

Cls
Local pixm:TPixmap = LoadPixmap("testimage.png")
Local time:Int = MilliSecs()
Local pixm2:TPixmap = gaussianBlur(pixm,9,9)
Print "processtime was: "+(MilliSecs()-time)+"ms"
DrawPixmap(pixm,0,0)
DrawPixmap(pixm2,256,0)
Flip
WaitKey()
End

' pixmap gaussian blur function
' based on blitzbasic code by Elias_T
' maskWidth & maskHeight should always be uneven numbers [3,5,7,9]
Function gaussianBlur:TPixmap(pMap:TPixmap,maskWidth:Int,maskHeight:Int)

Local tmp:TPixmap = CopyPixmap(pMap)
Local pMapPtr:Byte Ptr = PixmapPixelPtr(pMap,0,0)
Local tmpPtr:Byte Ptr = PixmapPixelPtr(tmp,0,0)

Local width:Int = pMap.width
Local height:Int = pMap.height
Local texel:Float[width,height,3]
Local result:Float[width,height,3]
Local maskData:Float[width*height]

Local x:Int,y:Int,ym:Int,xm:Int
Local cy:Int,cx:Int,rt:Float
Local r1:Float,g1:Float,b1:Float
Local rr:Float,gg:Float,bb:Float
Local mult:Float = 0.0
Local mwhalf:Int,mhhalf:Int
Local y_index:Int,x_index:Int,mask_index:Int

mwhalf = Floor((maskWidth-1) / 2)
mhhalf = Floor((maskHeight-1) / 2)

For x = 0 To width-1
For y = 0 To height-1
texel[x,y,0] = pMapPtr[x*4+y*4*width]
texel[x,y,1] = pMapPtr[x*4+y*4*width+1]
texel[x,y,2] = pMapPtr[x*4+y*4*width+2]
Next
Next

For ym = 0 To maskHeight-1
For xm = 0 To maskWidth-1
cx = xm - mwhalf
cy = ym - mhhalf
rt = Float(cx*cx + cy*cy)
mult :+ Exp(-0.35 * rt)
Next
Next

mult = 1.0 / mult

For ym = 0 To maskHeight-1
For xm = 0 To maskWidth-1
cx = xm - mwhalf
cy = ym - mhhalf
rt = Float(cx*cx + cy*cy)
maskData[ym * maskWidth + xm] = mult * Exp(-0.35 * rt)
Next
Next

mhhalf = Floor(maskHeight/2.0)
mwhalf = Floor(maskWidth/2.0)
For ym = 0 To height-1
For xm = 0 To width-1
rr = 0.0
gg = 0.0
bb = 0.0
y_index = ym-mhhalf
For yy = 0 To maskHeight-1
If y_index => 0 And y_index < height Then
mask_index = yy * maskwidth
x_index = xm-mwhalf
For xx = 0 To maskWidth-1
If x_index => 0 And x_index < width Then
rr :+ texel[x_index, y_index,0] * maskData[mask_index]
gg :+ texel[x_index, y_index,1] * maskData[mask_index]
bb :+ texel[x_index, y_index,2] * maskData[mask_index]
End If
mask_index :+ 1
x_index :+ 1
Next
End If
y_index :+ 1
Next
result[xm,ym,0] = rr
result[xm,ym,1] = gg
result[xm,ym,2] = bb
Next
Next

For x = 0 To width-1
For y = 0 To height-1
tmpPtr[x*4+y*4*width] = Int(result[x,y,0])
tmpPtr[x*4+y*4*width+1] = Int(result[x,y,1])
tmpPtr[x*4+y*4*width+2] = Int(result[x,y,2])
Next
Next

Return tmp
End Function


I am going to move on to other filters now. I am slowly working my way through a lot of the different filters used by a lot of art packages. Eventually they are to be used in my own paint program, but its a long way ahead and its taking its time. But I will get there eventually. Thanks for the help so far. If I run into a wall on some of the other filters, I will be back posing what I got to see if it can be fixed/improved.