Add integral image functionality to enable massive speedup of Sauvola

Note that there are some very small differences to the output compared
to the basic algorithm, but this doesn't make much difference.

This is due to minor differences with the standard deviation
calculation throughout, and with mean calculation at edges, for reasons
I'm unclear about.
WIP integral image speedup. mean is working

Very WIP, but mean is perfect once full window is used

Integral version all working!

Remove debugging info

Organise code better

1 files changed, 25 insertions, 62 deletions

diff --git a/binarize/sauvola.go b/binarize/sauvola.go
index f1d0512..bc311ad 100644
--- a/binarize/sauvola.go
+++ b/binarize/sauvola.go
@@ -3,81 +3,44 @@ package main
 import (
 	"image"
 	"image/color"
-	"math"
 )
 
-func mean(i []int) float64 {
-	sum := 0
-	for _, n := range i {
-		sum += n
-	}
-	return float64(sum) / float64(len(i))
-}
-
-// TODO: is there a prettier way of doing this than float64() all over the place?
-func stddev(i []int) float64 {
-	m := mean(i)
-
-	var sum float64
-	for _, n := range i {
-		sum += (float64(n) - m) * (float64(n) - m)
-	}
-	variance := float64(sum) / float64(len(i) - 1)
-	return math.Sqrt(variance)
-}
-
-func meanstddev(i []int) (float64, float64) {
-	m := mean(i)
-
-	var sum float64
-	for _, n := range i {
-		sum += (float64(n) - m) * (float64(n) - m)
-	}
-	variance := float64(sum) / float64(len(i) - 1)
-	return m, math.Sqrt(variance)
-}
-
-// gets the pixel values surrounding a point in the image
-func surrounding(img *image.Gray, x int, y int, size int) []int {
+// Implements Sauvola's algorithm for text binarization, see paper
+// "Adaptive document image binarization" (2000)
+func Sauvola(img *image.Gray, ksize float64, windowsize int) *image.Gray {
 	b := img.Bounds()
+	new := image.NewGray(b)
 
-	miny := y - size/2
-	if miny < b.Min.Y {
-		miny = b.Min.Y
-	}
-	minx := x - size/2
-	if minx < b.Min.X {
-		minx = b.Min.X
-	}
-	maxy := y + size/2
-	if maxy > b.Max.Y {
-		maxy = b.Max.Y
-	}
-	maxx := x + size/2
-	if maxx > b.Max.X {
-		maxx = b.Max.X
-	}
-
-	var s []int
-	for yi := miny; yi < maxy; yi++ {
-		for xi := minx; xi < maxx; xi++ {
-			s = append(s, int(img.GrayAt(xi, yi).Y))
+	for y := b.Min.Y; y < b.Max.Y; y++ {
+		for x := b.Min.X; x < b.Max.X; x++ {
+			window := surrounding(img, x, y, windowsize)
+			m, dev := meanstddev(window)
+			threshold := m * (1 + ksize * ((dev / 128) - 1))
+			if img.GrayAt(x, y).Y < uint8(threshold) {
+				new.SetGray(x, y, color.Gray{0})
+			} else {
+				new.SetGray(x, y, color.Gray{255})
+			}
 		}
 	}
-	return s
+
+	return new
 }
 
-// TODO: parallelize
-// TODO: switch to using integral images to make faster; see paper
-//       "Efficient Implementation of Local Adaptive Thresholding Techniques Using Integral Images"
-func Sauvola(img *image.Gray, ksize float64, windowsize int) *image.Gray {
+// Implements Sauvola's algorithm using Integral Images, see paper
+// "Effcient Implementation of Local Adaptive Thresholding Techniques Using Integral Images"
+// and
+// https://stackoverflow.com/questions/13110733/computing-image-integral
+func IntegralSauvola(img *image.Gray, ksize float64, windowsize int) *image.Gray {
 	b := img.Bounds()
 	new := image.NewGray(b)
 
+	integral := integralimg(img)
+	integralsq := integralimgsq(img)
+
 	for y := b.Min.Y; y < b.Max.Y; y++ {
 		for x := b.Min.X; x < b.Max.X; x++ {
-			window := surrounding(img, x, y, windowsize)
-			m, dev := meanstddev(window)
+			m, dev := integralmeanstddev(integral, integralsq, x, y, windowsize)
 			threshold := m * (1 + ksize * ((dev / 128) - 1))
 			if img.GrayAt(x, y).Y < uint8(threshold) {
 				new.SetGray(x, y, color.Gray{0})