Calculate 1D Gaussian Kernel?

BlitzMax Forums/BlitzMax Programming/Calculate 1D Gaussian Kernel?

Gabriel(Posted 2009) [#1]
Hi,

I'm trying to tweak my blur shaders a bit, and I decided to try a gaussian blur. I separate the blur into two passes which means I can get a 9x9 kernel with 18 samples instead of 81, and it also means I need a 1d kernel. The only algorithm I managed to find was for a 2d kernel, and had a couple of symbols I didn't recognize anyway. (Always my problem.)

Could anyone point me to some code to calculate a 1D gaussian kernel of a variable size? As I say, I'm using a size of 9 right now, but that might change so I had better find the code rather than just a kernel someone created for me.


Gabriel(Posted 2009) [#2]
I found some source code in Java which I quickly ported over. It produces a kernel whose total is not 1, which I guess is ok in some applications, but probably not what I want here. So I added a tweak function which scales the 1d output array so that the total is as close to 1 as possible.

If someone has better code, then by all means share.

Function MakeKernel:Float[] (Radius:Float)
	Local R:Int = Ceil(Radius)
	Local Rows:Int = R * 2 + 1
	Local Matrix:Float[] = New Float[Rows]
	Local Sigma:Float = Radius / 3.0
	Local Sigma22:Float = 2 * Sigma * Sigma
	Local SigmaPi2:Float = 2 * Pi * Sigma
	Local SqrtSigmaPi2:Float = Sqr(SigmaPi2)
	Local Radius2:Float = Radius * Radius
	Local Total:Float = 0
	Local Index:Int = 0
	For Local Row:Int = -R To R
		Local Distance:Float = Row * Row
		If (Distance > Radius2)
			Matrix[Index] = 0
		Else
			Matrix[Index] = Exp(- (Distance) / Sigma22) / SqrtSigmaPi2
		End If
		Total:+Matrix[Index]
		Index:+1
	Next
	Return Matrix
End Function

Function TweakKernel:Float[](In:Float[])
	Local Out:Float[]=New Float[In.Length]
	Local Total:Float=0
	For Local I:Int=0 To In.Length-1
		Total:+In[I]
	Next
	Print "Total:"+Total
	For Local J:Int=0 To In.Length-1
		Out[J]=In[J]/Total
		Print In[J]+" Becomes "+ Out[J]
	Next
	Return Out
End Function


Local Mat:Float[]=MakeKernel(4)
Local Mat2:Float[]=TweakKernel(Mat)
CheckKernel(Mat2)



Warpy(Posted 2009) [#3]
The main reason you had to tweak it was that the normalisation coefficient wasn't quite right. The code you found included sigma in the square root, when it should in fact be outside. Anyway, there's still a slight error in the

From a quick google for "gaussian kernel", here's the wikipedia page which has the 1d and 2d equations on it, and here's a pdf which is very interesting to me and explains the maths behind it.

So, it's just a matter of writing out the equation. The only bit of trouble I had was with picking the right sigma value, but the wikipedia page says
In practice, when computing a discrete approximation of the Gaussian function, pixels at a distance of more than 3*sigma are small enough to be considered effectively zero.

so that's that solved.

Function GaussianKernel#[](radius#)
	width=Ceil(radius)	'kernel will have a middle cell, and width on either side
	Local matrix#[width*2+1]
	sigma# = radius/3		'apparently this is all you need to get a good approximation
	norm# = 1.0 / (Sqr(2*Pi) * sigma)		'normalisation constant makes sure total of matrix is 1
	coeff# = 2*sigma*sigma	'the bit you divide x^2 by in the exponential
	total#=0
	For x = -width To width	'fill in matrix!
		g# = norm * Exp( -x*x/coeff )
		matrix[x+width] = g
		total:+g
	Next
	For x=0 To 2*width	'rescale things to get a total of 1, because of discretisation error
		matrix[x]:/total
	Next
	Return matrix
End Function


Enjoy!


Gabriel(Posted 2009) [#4]
Thanks, that looks great. I did convert that java code at about 5am, so it's possible that it was my mistake when converting. I don't want the original author to feel slighted if it was my error.