Make gaussian blur using separable convolution

Author: Angelos Katharopoulos

blur/blur.go

@@ -35,19 +35,18 @@ func Gaussian(src image.Image, radius float64) *image.RGBA {
 		return clone.AsRGBA(src)
 	}
 
+    // Create the 1-d gaussian kernel
 	length := int(math.Ceil(2*radius + 1))
-	k := convolution.NewKernel(length, length)
-
-	gaussianFn := func(x, y, sigma float64) float64 {
-		return math.Exp(-(x*x/sigma + y*y/sigma))
-	}
-
-	for x := 0; x < length; x++ {
-		for y := 0; y < length; y++ {
-			k.Matrix[y*length+x] = gaussianFn(float64(x)-radius, float64(y)-radius, 4*radius)
-
-		}
-	}
-
-	return convolution.Convolve(src, k.Normalized(), &convolution.Options{Bias: 0, Wrap: false, KeepAlpha: false})
+	k := convolution.NewKernel(length, 1)
+    for i, x := 0, -radius; i < length; i, x = i+1, x+1 {
+        k.Matrix[i] = math.Exp(-(x*x/4/radius))
+    }
+    normK := k.Normalized()
+
+    // Perform separable convolution
+    options := convolution.Options{Bias: 0, Wrap: false, KeepAlpha: false}
+    result := convolution.Convolve(src, normK, &options)
+    result = convolution.Convolve(result, normK.Transpose(), &options)
+
+    return result
 }

blur/blur_test.go

@@ -128,9 +128,9 @@ func TestGaussianBlur(t *testing.T) {
 				Rect:   image.Rect(0, 0, 3, 3),
 				Stride: 3 * 4,
 				Pix: []uint8{
-					0xb1, 0xb1, 0xb1, 0xff, 0xbd, 0xbd, 0xbd, 0xff, 0xcc, 0xcc, 0xcc, 0xff,
-					0x4d, 0x4d, 0x4d, 0xff, 0x68, 0x68, 0x68, 0xff, 0x8b, 0x8b, 0x8b, 0xff,
-					0x0, 0x0, 0x0, 0xff, 0x26, 0x26, 0x26, 0xff, 0x59, 0x59, 0x59, 0xff,
+                    0xb1, 0xb1, 0xb1, 0xff, 0xbc, 0xbc, 0xbc, 0xff, 0xcc, 0xcc, 0xcc, 0xff,
+                    0x4d, 0x4d, 0x4d, 0xff, 0x68, 0x68, 0x68, 0xff, 0x8b, 0x8b, 0x8b, 0xff,
+                    0x0, 0x0, 0x0, 0xff, 0x25, 0x25, 0x25, 0xff, 0x58, 0x58, 0x58, 0xff,
 				},
 			},
 		},
@@ -157,10 +157,10 @@ func TestGaussianBlur(t *testing.T) {
 		},
 	}
 
-	for _, c := range cases {
+	for i, c := range cases {
 		actual := Gaussian(c.value, c.radius)
 		if !util.RGBAImageEqual(actual, c.expected) {
-			t.Errorf("%s: expected: %#v, actual: %#v", "GaussianBlur", c.expected, actual)
+			t.Errorf("%s %d: expected: %#v, actual: %#v", "GaussianBlur", i, c.expected, actual)
 		}
 	}
 }

convolution/kernel.go

@@ -15,6 +15,7 @@ type Matrix interface {
 	Normalized() Matrix
 	MaxX() int
 	MaxY() int
+    Transpose() Matrix
 }
 
 // NewKernel returns a kernel of the provided length.
@@ -63,6 +64,19 @@ func (k *Kernel) At(x, y int) float64 {
 	return k.Matrix[y*k.Width+x]
 }
 
+// Transpose returns a new Kernel that is transposed
+func (k *Kernel) Transpose() Matrix {
+    w := k.Width;
+    h := k.Height;
+    nk := NewKernel(h, w)
+
+    for i := 0; i < w*h; i++ {
+        nk.Matrix[i] = k.Matrix[i];
+    }
+
+    return nk
+}
+
 // String returns the string representation of the matrix.
 func (k *Kernel) String() string {
 	result := ""