commit 428642c9f1e452287f8cb4d9789282ae6a1afca9 Author: nfnt Date: Thu Aug 2 21:59:40 2012 +0200 initial commit diff --git a/README b/README new file mode 100644 index 0000000..e69de29 diff --git a/filters.go b/filters.go new file mode 100644 index 0000000..27c321f --- /dev/null +++ b/filters.go @@ -0,0 +1,142 @@ +/* +Copyright (c) 2012, Jan Schlicht + +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]} +} diff --git a/resize.go b/resize.go new file mode 100644 index 0000000..36fd2f6 --- /dev/null +++ b/resize.go @@ -0,0 +1,108 @@ +/* +Copyright (c) 2012, Jan Schlicht + +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 +} diff --git a/resize_test.go b/resize_test.go new file mode 100644 index 0000000..c40aed2 --- /dev/null +++ b/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() + } +} diff --git a/sinc.go b/sinc.go new file mode 100644 index 0000000..723e6af --- /dev/null +++ b/sinc.go @@ -0,0 +1,49 @@ +/* +Copyright (c) 2012, Jan Schlicht + +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) +} diff --git a/sinc_test.go b/sinc_test.go new file mode 100644 index 0000000..9372853 --- /dev/null +++ b/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") + } +}