initial commit

filters.go

@@ -0,0 +1,142 @@
+/*
+Copyright (c) 2012, Jan Schlicht <jan.schlicht@gmail.com>
+
+Permission to use, copy, modify, and/or distribute this software for any purpose
+with or without fee is hereby granted, provided that the above copyright notice
+and this permission notice appear in all copies.
+
+THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
+REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
+FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
+INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
+OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
+TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+THIS SOFTWARE.
+*/
+
+package resize
+
+import (
+	"image"
+	"image/color"
+	"math"
+)
+
+// color.RGBA64 as array
+type RGBA [4]uint16
+
+// build RGBA from an arbitrary color
+func toRGBA(c color.Color) RGBA {
+	n := color.RGBA64Model.Convert(c).(color.RGBA64)
+	return RGBA{n.R, n.G, n.B, n.A}
+}
+
+func clampToUint16(x float32) (y uint16) {
+	y = uint16(x)
+	if x < 0 {
+		y = 0
+	} else if x > float32(0xffff) {
+		y = 0xffff
+	}
+	return
+}
+
+// Nearest-neighbor interpolation.
+// Approximates a value by returning the value of the nearest point.
+func NearestNeighbor(x, y float32, img image.Image) color.RGBA64 {
+	xn, yn := int(x), int(y)
+	c := toRGBA(img.At(xn, yn))
+	return color.RGBA64{c[0], c[1], c[2], c[3]}
+}
+
+// Linear interpolation.
+func linearInterp(x float32, p *[2]RGBA) (c RGBA) {
+	x -= float32(math.Floor(float64(x)))
+	for i := range c {
+		c[i] = clampToUint16(float32(p[0][i])*(1.0-x) + x*float32(p[1][i]))
+	}
+	return
+}
+
+// Bilinear interpolation.
+func Bilinear(x, y float32, img image.Image) color.RGBA64 {
+	xf, yf := int(math.Floor(float64(x))), int(math.Floor(float64(y)))
+
+	var row [2]RGBA
+	var col [2]RGBA
+	row = [2]RGBA{toRGBA(img.At(xf, yf)), toRGBA(img.At(xf+1, yf))}
+	col[0] = linearInterp(x, &row)
+	row = [2]RGBA{toRGBA(img.At(xf, yf+1)), toRGBA(img.At(xf+1, yf+1))}
+	col[1] = linearInterp(x, &row)
+
+	c := linearInterp(y, &col)
+	return color.RGBA64{c[0], c[1], c[2], c[3]}
+}
+
+// cubic interpolation
+func cubicInterp(x float32, p *[4]RGBA) (c RGBA) {
+	x -= float32(math.Floor(float64(x)))
+	for i := range c {
+		c[i] = clampToUint16(float32(p[1][i]) + 0.5*x*(float32(p[2][i])-float32(p[0][i])+x*(2.0*float32(p[0][i])-5.0*float32(p[1][i])+4.0*float32(p[2][i])-float32(p[3][i])+x*(3.0*(float32(p[1][i])-float32(p[2][i]))+float32(p[3][i])-float32(p[0][i])))))
+	}
+	return
+}
+
+// Bicubic interpolation.
+func Bicubic(x, y float32, img image.Image) color.RGBA64 {
+	xf, yf := int(math.Floor(float64(x))), int(math.Floor(float64(y)))
+
+	var row [4]RGBA
+	var col [4]RGBA
+	row = [4]RGBA{toRGBA(img.At(xf-1, yf-1)), toRGBA(img.At(xf, yf-1)), toRGBA(img.At(xf+1, yf-1)), toRGBA(img.At(xf+2, yf-1))}
+	col[0] = cubicInterp(x, &row)
+	row = [4]RGBA{toRGBA(img.At(xf-1, yf)), toRGBA(img.At(xf, yf)), toRGBA(img.At(xf+1, yf)), toRGBA(img.At(xf+2, yf))}
+	col[1] = cubicInterp(x, &row)
+	row = [4]RGBA{toRGBA(img.At(xf-1, yf+1)), toRGBA(img.At(xf, yf+1)), toRGBA(img.At(xf+1, yf+1)), toRGBA(img.At(xf+2, yf+1))}
+	col[2] = cubicInterp(x, &row)
+	row = [4]RGBA{toRGBA(img.At(xf-1, yf+2)), toRGBA(img.At(xf, yf+2)), toRGBA(img.At(xf+1, yf+2)), toRGBA(img.At(xf+2, yf+2))}
+	col[3] = cubicInterp(x, &row)
+
+	c := cubicInterp(y, &col)
+	return color.RGBA64{c[0], c[1], c[2], c[3]}
+}
+
+// 1-d convolution with windowed sinc for a=3.
+func lanczos_x(x float32, p *[6]RGBA) (c RGBA) {
+	x -= float32(math.Floor(float64(x)))
+	var v float32
+	l := [4]float32{0.0, 0.0, 0.0, 0.0}
+	for j := range p {
+		v = float32(Sinc(float64(x-float32(j-2)))) * float32(Sinc(float64((x-float32(j-2))/3.0)))
+		for i := range c {
+			l[i] += float32(p[j][i]) * v
+		}
+	}
+	for i := range c {
+		c[i] = clampToUint16(l[i])
+	}
+	return
+}
+
+// Lanczos interpolation (a=3).
+func Lanczos3(x, y float32, img image.Image) color.RGBA64 {
+	xf, yf := int(math.Floor(float64(x))), int(math.Floor(float64(y)))
+
+	var row [6]RGBA
+	var col [6]RGBA
+	row = [6]RGBA{toRGBA(img.At(xf-2, yf-2)), toRGBA(img.At(xf-1, yf-2)), toRGBA(img.At(xf, yf-2)), toRGBA(img.At(xf+1, yf-2)), toRGBA(img.At(xf+2, yf-2)), toRGBA(img.At(xf+3, yf-2))}
+	col[0] = lanczos_x(x, &row)
+	row = [6]RGBA{toRGBA(img.At(xf-2, yf-1)), toRGBA(img.At(xf-1, yf-1)), toRGBA(img.At(xf, yf-1)), toRGBA(img.At(xf+1, yf-1)), toRGBA(img.At(xf+2, yf-1)), toRGBA(img.At(xf+3, yf-1))}
+	col[1] = lanczos_x(x, &row)
+	row = [6]RGBA{toRGBA(img.At(xf-2, yf)), toRGBA(img.At(xf-1, yf)), toRGBA(img.At(xf, yf)), toRGBA(img.At(xf+1, yf)), toRGBA(img.At(xf+2, yf)), toRGBA(img.At(xf+3, yf))}
+	col[2] = lanczos_x(x, &row)
+	row = [6]RGBA{toRGBA(img.At(xf-2, yf+1)), toRGBA(img.At(xf-1, yf+1)), toRGBA(img.At(xf, yf+1)), toRGBA(img.At(xf+1, yf+1)), toRGBA(img.At(xf+2, yf+1)), toRGBA(img.At(xf+3, yf+1))}
+	col[3] = lanczos_x(x, &row)
+	row = [6]RGBA{toRGBA(img.At(xf-2, yf+2)), toRGBA(img.At(xf-1, yf+2)), toRGBA(img.At(xf, yf+2)), toRGBA(img.At(xf+1, yf+2)), toRGBA(img.At(xf+2, yf+2)), toRGBA(img.At(xf+3, yf+2))}
+	col[4] = lanczos_x(x, &row)
+	row = [6]RGBA{toRGBA(img.At(xf-2, yf+3)), toRGBA(img.At(xf-1, yf+3)), toRGBA(img.At(xf, yf+3)), toRGBA(img.At(xf+1, yf+3)), toRGBA(img.At(xf+2, yf+3)), toRGBA(img.At(xf+3, yf+3))}
+	col[5] = lanczos_x(x, &row)
+
+	c := lanczos_x(y, &col)
+	return color.RGBA64{c[0], c[1], c[2], c[3]}
+}

resize.go

@@ -0,0 +1,108 @@
+/*
+Copyright (c) 2012, Jan Schlicht <jan.schlicht@gmail.com>
+
+Permission to use, copy, modify, and/or distribute this software for any purpose
+with or without fee is hereby granted, provided that the above copyright notice
+and this permission notice appear in all copies.
+
+THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
+REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
+FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
+INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
+OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
+TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+THIS SOFTWARE.
+*/
+
+// Package resize implements various image resizing methods.
+//
+// The package works with the Image interface described in the image package.
+// Various interpolation methods are provided and multiple processors may be
+// utilized in the computations.
+//
+// Example:
+//     imgResized := resize.Resize(1000, -1, imgOld, Lanczos3)
+package resize
+
+import (
+	"image"
+	"image/color"
+	"runtime"
+)
+
+var (
+	// NCPU holds the number of available CPUs at runtime.
+	NCPU = runtime.NumCPU()
+)
+
+// Trans2 is a 2-dimensional linear transformation.
+type Trans2 [6]float32
+
+// Apply the transformation to a point (x,y).
+func (t *Trans2) Eval(x, y float32) (u, v float32) {
+	u = t[0]*x + t[1]*y + t[2]
+	v = t[3]*x + t[4]*y + t[5]
+	return
+}
+
+// Calculate scaling factors using old and new image dimensions.
+func calcFactors(w, h int, wo, ho float32) (sx, sy float32) {
+	if w == -1 {
+		if h == -1 {
+			sx = 1.0
+			sy = 1.0
+		} else {
+			sy = ho / float32(h)
+			sx = sy
+		}
+	} else {
+		sx = wo / float32(w)
+		if h == -1 {
+			sy = sx
+		} else {
+			sy = ho / float32(h)
+		}
+	}
+	return
+}
+
+// InterpolationFunction return a color for an arbitrary point inside
+// an image
+type InterpolationFunction func(float32, float32, image.Image) color.RGBA64
+
+// Resize an image to new width w and height h using the interpolation function interp.
+// A new image with the given dimensions will be returned.
+// If one of the parameters w or h is set to -1, its size will be calculated so that
+// the aspect ratio is that of the originating image.
+// The resizing algorithm uses slices for parallel computation.
+func Resize(w int, h int, img image.Image, interp InterpolationFunction) image.Image {
+	b_old := img.Bounds()
+	w_old := float32(b_old.Dx())
+	h_old := float32(b_old.Dy())
+
+	scaleX, scaleY := calcFactors(w, h, w_old, h_old)
+	t := Trans2{scaleX, 0, float32(b_old.Min.X), 0, scaleY, float32(b_old.Min.Y)}
+
+	m := image.NewRGBA64(image.Rect(0, 0, int(w_old/scaleX), int(h_old/scaleY)))
+	b := m.Bounds()
+
+	c := make(chan int, NCPU)
+	for i := 0; i < NCPU; i++ {
+		go func(b image.Rectangle, c chan int) {
+			var u, v float32
+			for y := b.Min.Y; y < b.Max.Y; y++ {
+				for x := b.Min.X; x < b.Max.X; x++ {
+					u, v = t.Eval(float32(x), float32(y))
+					m.SetRGBA64(x, y, interp(u, v, img))
+				}
+			}
+			c <- 1
+		}(image.Rect(b.Min.X, b.Min.Y+i*(b.Dy())/4, b.Max.X, b.Min.Y+(i+1)*(b.Dy())/4), c)
+	}
+
+	for i := 0; i < NCPU; i++ {
+		<-c
+	}
+
+	return m
+}

resize_test.go

@@ -0,0 +1,18 @@
+package resize
+
+import (
+	"image"
+	"image/color"
+	"testing"
+)
+
+func Test_Nearest(t *testing.T) {
+	img := image.NewGray16(image.Rect(0,0, 3,3))
+	img.Set(1,1, color.White)
+	
+	m := Resize(6,-1, img, NearestNeighbor)
+	
+	if m.At(2,2) != m.At(3,3) {
+		t.Fail()
+	}
+}

sinc.go

@@ -0,0 +1,49 @@
+/*
+Copyright (c) 2012, Jan Schlicht <jan.schlicht@gmail.com>
+
+Permission to use, copy, modify, and/or distribute this software for any purpose
+with or without fee is hereby granted, provided that the above copyright notice
+and this permission notice appear in all copies.
+
+THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
+REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
+FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
+INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
+OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
+TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+THIS SOFTWARE.
+*/
+
+package resize
+
+import (
+	"math"
+)
+
+var (
+	epsilon      = math.Nextafter(1.0, 2.0) - 1.0 // machine epsilon
+	taylor2bound = math.Sqrt(epsilon)
+	taylorNbound = math.Sqrt(taylor2bound)
+)
+
+// unnormalized sinc function
+func Sinc1(x float64) (y float64) {
+	if math.Abs(x) >= taylorNbound {
+		y = math.Sin(x) / x
+	} else {
+		y = 1.0
+		if math.Abs(x) >= epsilon {
+			x2 := x * x
+			y -= x2 / 6.0
+			if math.Abs(x) >= taylor2bound {
+				y += (x2 * x2) / 120.0
+			}
+		}
+	}
+	return
+}
+
+// normalized sinc function
+func Sinc(x float64) float64 {
+	return Sinc1(x * math.Pi)
+}

sinc_test.go

@@ -0,0 +1,38 @@
+package resize
+
+import (
+	"fmt"
+	"math"
+	"testing"
+)
+
+const limit = 1e-12
+
+func Test_SincOne(t *testing.T) {
+	zero := Sinc(1)
+	if zero >= limit {
+		t.Error("Sinc(1) != 0")
+	}
+}
+
+func Test_SincZero(t *testing.T) {
+	one := Sinc(0)
+	if math.Abs(one-1) >= limit {
+		t.Error("Sinc(0) != 1")
+	}
+}
+
+func Test_SincDotOne(t *testing.T) {
+	res := Sinc(0.1)
+	if math.Abs(res-0.983631643083466) >= limit {
+		t.Error("Sinc(0.1) wrong")
+	}
+}
+
+func Test_SincNearZero(t *testing.T) {
+	res := Sinc(0.000001)
+	if math.Abs(res-0.9999999999983551) >= limit {
+		fmt.Println(res)
+		t.Error("Sinc near zero not stable")
+	}
+}