summaryrefslogtreecommitdiff
path: root/integral.go
diff options
context:
space:
mode:
Diffstat (limited to 'integral.go')
-rw-r--r--integral.go213
1 files changed, 213 insertions, 0 deletions
diff --git a/integral.go b/integral.go
new file mode 100644
index 0000000..e3d09bb
--- /dev/null
+++ b/integral.go
@@ -0,0 +1,213 @@
+// Copyright 2020 Nick White.
+// Use of this source code is governed by the GPLv3
+// license that can be found in the LICENSE file.
+
+// integral is a package for processing integral images, aka
+// summed area tables. These are structures which precompute the
+// sum of pixels to the left and above each pixel, which can make
+// several common image processing operations much faster.
+//
+// A lot of image processing operations rely on many calculations
+// of the sum or mean of a set of pixels. As these have been
+// precalculated for an integral image, these calculations are
+// much faster. Image.Sum() and Image.Mean() functions are provided
+// by this package to take advantage of this.
+//
+// Another common requirement is standard deviation over an area
+// of an image. This can be calculated by creating an integral
+// image and squared integral image (SqImage) for a base image, and
+// passing them to the MeanStdDev() function provided.
+package integral
+
+import (
+ "image"
+ "image/color"
+ "math"
+)
+
+// Image is an integral image
+type Image [][]uint64
+
+// SqImage is a Square integral image.
+// A squared integral image is an integral image for which the square of
+// each pixel is saved; this is useful for efficiently calculating
+// Standard Deviation.
+type SqImage [][]uint64
+
+func (i Image) ColorModel() color.Model { return color.Gray16Model }
+
+func (i Image) Bounds() image.Rectangle {
+ return image.Rect(0, 0, len(i[0]), len(i))
+}
+
+// at64 is used to return the raw uint64 for a given pixel. Accessing
+// this separately to a (potentially lossy) conversion to a Gray16 is
+// necessary for SqImage to function accurately.
+func (i Image) at64(x, y int) uint64 {
+ if !(image.Point{x, y}.In(i.Bounds())) {
+ return 0
+ }
+
+ var prevx, prevy, prevxy uint64
+ prevx, prevy, prevxy = 0, 0, 0
+ if x > 0 {
+ prevx = i[y][x-1]
+ }
+ if y > 0 {
+ prevy = i[y-1][x]
+ }
+ if x > 0 && y > 0 {
+ prevxy = i[y-1][x-1]
+ }
+ orig := i[y][x] + prevxy - prevx - prevy
+ return orig
+}
+
+func (i Image) At(x, y int) color.Color {
+ c := i.at64(x, y)
+ return color.Gray16{uint16(c)}
+}
+
+func (i Image) set64(x, y int, c uint64) {
+ var prevx, prevy, prevxy uint64
+ prevx, prevy, prevxy = 0, 0, 0
+ if x > 0 {
+ prevx = i[y][x-1]
+ }
+ if y > 0 {
+ prevy = i[y-1][x]
+ }
+ if x > 0 && y > 0 {
+ prevxy = i[y-1][x-1]
+ }
+ final := c + prevx + prevy - prevxy
+ i[y][x] = final
+}
+
+func (i Image) Set(x, y int, c color.Color) {
+ gray := color.Gray16Model.Convert(c).(color.Gray16).Y
+ i.set64(x, y, uint64(gray))
+}
+
+// NewImage returns a new integral image with the given bounds.
+func NewImage(r image.Rectangle) *Image {
+ w, h := r.Dx(), r.Dy()
+ var rows Image
+ for i := 0; i < h; i++ {
+ col := make([]uint64, w)
+ rows = append(rows, col)
+ }
+ return &rows
+}
+
+func (i SqImage) ColorModel() color.Model { return Image(i).ColorModel() }
+
+func (i SqImage) Bounds() image.Rectangle {
+ return Image(i).Bounds()
+}
+
+func (i SqImage) At(x, y int) color.Color {
+ c := Image(i).at64(x, y)
+ rt := math.Sqrt(float64(c))
+ return color.Gray16{uint16(rt)}
+}
+
+func (i SqImage) Set(x, y int, c color.Color) {
+ gray := uint64(color.Gray16Model.Convert(c).(color.Gray16).Y)
+ Image(i).set64(x, y, gray*gray)
+}
+
+// NewSqImage returns a new squared integral image with the given bounds.
+func NewSqImage(r image.Rectangle) *SqImage {
+ i := NewImage(r)
+ s := SqImage(*i)
+ return &s
+}
+
+func lowest(a, b int) int {
+ if a < b {
+ return a
+ }
+ return b
+}
+
+func highest(a, b int) int {
+ if a > b {
+ return a
+ }
+ return b
+}
+
+func (i Image) topLeft(r image.Rectangle) uint64 {
+ b := i.Bounds()
+ x := r.Min.X - 1
+ y := r.Min.Y - 1
+ x = lowest(x, b.Max.X-1)
+ y = lowest(y, b.Max.Y-1)
+ if x < 0 || y < 0 {
+ return 0
+ }
+ return i[y][x]
+}
+
+func (i Image) topRight(r image.Rectangle) uint64 {
+ b := i.Bounds()
+ x := lowest(r.Max.X-1, b.Max.X-1)
+ y := r.Min.Y - 1
+ y = lowest(y, b.Max.Y-1)
+ if x < 0 || y < 0 {
+ return 0
+ }
+ return i[y][x]
+}
+
+func (i Image) bottomLeft(r image.Rectangle) uint64 {
+ b := i.Bounds()
+ x := r.Min.X - 1
+ x = lowest(x, b.Max.X-1)
+ y := lowest(r.Max.Y-1, b.Max.Y-1)
+ if x < 0 || y < 0 {
+ return 0
+ }
+ return i[y][x]
+}
+
+func (i Image) bottomRight(r image.Rectangle) uint64 {
+ b := i.Bounds()
+ x := lowest(r.Max.X-1, b.Max.X-1)
+ y := lowest(r.Max.Y-1, b.Max.Y-1)
+ return i[y][x]
+}
+
+// Sum returns the sum of all pixels in a section of an image
+func (i Image) Sum(r image.Rectangle) uint64 {
+ return i.bottomRight(r) + i.topLeft(r) - i.topRight(r) - i.bottomLeft(r)
+}
+
+// Mean returns the average value of pixels in a section of an image
+func (i Image) Mean(r image.Rectangle) float64 {
+ in := r.Intersect(i.Bounds())
+ return float64(i.Sum(r)) / float64(in.Dx()*in.Dy())
+}
+
+// Sum returns the sum of all pixels in a section of an image
+func (i SqImage) Sum(r image.Rectangle) uint64 {
+ return Image(i).Sum(r)
+}
+
+// Mean returns the average value of pixels in a section of an image
+func (i SqImage) Mean(r image.Rectangle) float64 {
+ return Image(i).Mean(r)
+}
+
+// MeanStdDev calculates the mean and standard deviation of a
+// section of an image, using the corresponding regular and square
+// integral images.
+func MeanStdDev(i Image, sq SqImage, r image.Rectangle) (float64, float64) {
+ imean := i.Mean(r)
+ smean := sq.Mean(r)
+
+ variance := smean - (imean * imean)
+
+ return imean, math.Sqrt(variance)
+}