Browse Source

Optimize data-locality for a huge increase in processing speed.

This is a complete rewrite! The tight scaling loop needs data locality for optimal performance. The old version used lots of pointer redirections to access image data which was bad for data locality. By providing the complete loop for each image type, this problem is solved. Unfortunately this increases code duplication but the result should be worth it: I could measure a ~6x speed-up for certain test cases!
jst 10 years ago
parent
commit
016a61cd31
6 changed files with 525 additions and 446 deletions
  1. 210
    84
      converter.go
  2. 65
    187
      filters.go
  3. 232
    81
      resize.go
  4. 18
    7
      resize_test.go
  5. 0
    49
      sinc.go
  6. 0
    38
      sinc_test.go

+ 210
- 84
converter.go View File

@@ -21,110 +21,236 @@ import (
21 21
 	"image/color"
22 22
 )
23 23
 
24
-type colorArray [4]float32
25
-
26
-func replicateBorder1d(x, min, max int) int {
27
-	if x < min {
28
-		x = min
29
-	} else if x >= max {
30
-		x = max - 1
24
+// Keep value in [0,255] range.
25
+func clampUint8(in int32) uint8 {
26
+	if in < 0 {
27
+		return 0
31 28
 	}
32
-
33
-	return x
29
+	if in > 255 {
30
+		return 255
31
+	}
32
+	return uint8(in)
34 33
 }
35 34
 
36
-func replicateBorder(x, y int, rect image.Rectangle) (xx, yy int) {
37
-	xx = replicateBorder1d(x, rect.Min.X, rect.Max.X)
38
-	yy = replicateBorder1d(y, rect.Min.Y, rect.Max.Y)
39
-	return
35
+// Keep value in [0,65535] range.
36
+func clampUint16(in int64) uint16 {
37
+	if in < 0 {
38
+		return 0
39
+	}
40
+	if in > 65535 {
41
+		return 65535
42
+	}
43
+	return uint16(in)
40 44
 }
41 45
 
42
-// converter allows to retrieve a colorArray for points of an image.
43
-// the idea is to speed up computation by providing optimized implementations
44
-// for different image types instead of relying on image.Image.At().
45
-type converter interface {
46
-	at(x, y int, color *colorArray)
47
-}
46
+func resizeGeneric(in image.Image, out *image.RGBA64, scale float64, coeffs []int32, filterLength int) {
47
+	oldBounds := in.Bounds()
48
+	newBounds := out.Bounds()
48 49
 
49
-type genericConverter struct {
50
-	src image.Image
51
-}
50
+	for x := newBounds.Min.X; x < newBounds.Max.X; x++ {
51
+		for y := newBounds.Min.Y; y < newBounds.Max.Y; y++ {
52
+			interpX := scale*(float64(y)+0.5) + float64(oldBounds.Min.X)
53
+			start := int(interpX) - filterLength/2 + 1
52 54
 
53
-func (c *genericConverter) at(x, y int, result *colorArray) {
54
-	r, g, b, a := c.src.At(replicateBorder(x, y, c.src.Bounds())).RGBA()
55
-	result[0] = float32(r)
56
-	result[1] = float32(g)
57
-	result[2] = float32(b)
58
-	result[3] = float32(a)
59
-	return
60
-}
55
+			var rgba [4]int64
56
+			var sum int64
57
+			for i := 0; i < filterLength; i++ {
58
+				xx := start + i
59
+				if xx < oldBounds.Min.X {
60
+					xx = oldBounds.Min.X
61
+				} else if xx >= oldBounds.Max.X {
62
+					xx = oldBounds.Max.X - 1
63
+				}
61 64
 
62
-type rgbaConverter struct {
63
-	src *image.RGBA
64
-}
65
+				coeff := coeffs[(y-newBounds.Min.Y)*filterLength+i]
66
+				r, g, b, a := in.At(xx, x).RGBA()
67
+				rgba[0] += int64(coeff) * int64(r)
68
+				rgba[1] += int64(coeff) * int64(g)
69
+				rgba[2] += int64(coeff) * int64(b)
70
+				rgba[3] += int64(coeff) * int64(a)
71
+				sum += int64(coeff)
72
+			}
65 73
 
66
-func (c *rgbaConverter) at(x, y int, result *colorArray) {
67
-	i := c.src.PixOffset(replicateBorder(x, y, c.src.Rect))
68
-	result[0] = float32(uint16(c.src.Pix[i+0])<<8 | uint16(c.src.Pix[i+0]))
69
-	result[1] = float32(uint16(c.src.Pix[i+1])<<8 | uint16(c.src.Pix[i+1]))
70
-	result[2] = float32(uint16(c.src.Pix[i+2])<<8 | uint16(c.src.Pix[i+2]))
71
-	result[3] = float32(uint16(c.src.Pix[i+3])<<8 | uint16(c.src.Pix[i+3]))
72
-	return
74
+			offset := (y-newBounds.Min.Y)*out.Stride + (x-newBounds.Min.X)*8
75
+			value := clampUint16(rgba[0] / sum)
76
+			out.Pix[offset+0] = uint8(value >> 8)
77
+			out.Pix[offset+1] = uint8(value)
78
+			value = clampUint16(rgba[1] / sum)
79
+			out.Pix[offset+2] = uint8(value >> 8)
80
+			out.Pix[offset+3] = uint8(value)
81
+			value = clampUint16(rgba[2] / sum)
82
+			out.Pix[offset+4] = uint8(value >> 8)
83
+			out.Pix[offset+5] = uint8(value)
84
+			value = clampUint16(rgba[3] / sum)
85
+			out.Pix[offset+6] = uint8(value >> 8)
86
+			out.Pix[offset+7] = uint8(value)
87
+		}
88
+	}
73 89
 }
74 90
 
75
-type rgba64Converter struct {
76
-	src *image.RGBA64
77
-}
91
+func resizeRGBA(in *image.RGBA, out *image.RGBA, scale float64, coeffs []int16, filterLength int) {
92
+	oldBounds := in.Bounds()
93
+	newBounds := out.Bounds()
78 94
 
79
-func (c *rgba64Converter) at(x, y int, result *colorArray) {
80
-	i := c.src.PixOffset(replicateBorder(x, y, c.src.Rect))
81
-	result[0] = float32(uint16(c.src.Pix[i+0])<<8 | uint16(c.src.Pix[i+1]))
82
-	result[1] = float32(uint16(c.src.Pix[i+2])<<8 | uint16(c.src.Pix[i+3]))
83
-	result[2] = float32(uint16(c.src.Pix[i+4])<<8 | uint16(c.src.Pix[i+5]))
84
-	result[3] = float32(uint16(c.src.Pix[i+6])<<8 | uint16(c.src.Pix[i+7]))
85
-	return
86
-}
95
+	for x := newBounds.Min.X; x < newBounds.Max.X; x++ {
96
+		row := in.Pix[(x-oldBounds.Min.Y)*in.Stride:]
97
+		for y := newBounds.Min.Y; y < newBounds.Max.Y; y++ {
98
+			interpX := scale*(float64(y)+0.5) + float64(oldBounds.Min.X)
99
+			start := int(interpX) - filterLength/2 + 1
87 100
 
88
-type grayConverter struct {
89
-	src *image.Gray
90
-}
101
+			var rgba [4]int32
102
+			var sum int32
103
+			for i := 0; i < filterLength; i++ {
104
+				xx := start + i
105
+				if xx < oldBounds.Min.X {
106
+					xx = oldBounds.Min.X
107
+				} else if xx >= oldBounds.Max.X {
108
+					xx = oldBounds.Max.X - 1
109
+				}
110
+
111
+				coeff := coeffs[(y-newBounds.Min.Y)*filterLength+i]
112
+				offset := (xx - oldBounds.Min.X) * 4
113
+				rgba[0] += int32(coeff) * int32(row[offset+0])
114
+				rgba[1] += int32(coeff) * int32(row[offset+1])
115
+				rgba[2] += int32(coeff) * int32(row[offset+2])
116
+				rgba[3] += int32(coeff) * int32(row[offset+3])
117
+				sum += int32(coeff)
118
+			}
91 119
 
92
-func (c *grayConverter) at(x, y int, result *colorArray) {
93
-	i := c.src.PixOffset(replicateBorder(x, y, c.src.Rect))
94
-	g := float32(uint16(c.src.Pix[i])<<8 | uint16(c.src.Pix[i]))
95
-	result[0] = g
96
-	result[1] = g
97
-	result[2] = g
98
-	result[3] = float32(0xffff)
99
-	return
120
+			offset := (y-newBounds.Min.Y)*out.Stride + (x-newBounds.Min.X)*4
121
+			out.Pix[offset+0] = clampUint8(rgba[0] / sum)
122
+			out.Pix[offset+1] = clampUint8(rgba[1] / sum)
123
+			out.Pix[offset+2] = clampUint8(rgba[2] / sum)
124
+			out.Pix[offset+3] = clampUint8(rgba[3] / sum)
125
+		}
126
+	}
100 127
 }
101 128
 
102
-type gray16Converter struct {
103
-	src *image.Gray16
129
+func resizeRGBA64(in *image.RGBA64, out *image.RGBA64, scale float64, coeffs []int32, filterLength int) {
130
+	oldBounds := in.Bounds()
131
+	newBounds := out.Bounds()
132
+
133
+	for x := newBounds.Min.X; x < newBounds.Max.X; x++ {
134
+		row := in.Pix[(x-oldBounds.Min.Y)*in.Stride:]
135
+		for y := newBounds.Min.Y; y < newBounds.Max.Y; y++ {
136
+			interpX := scale*(float64(y)+0.5) + float64(oldBounds.Min.X)
137
+			start := int(interpX) - filterLength/2 + 1
138
+
139
+			var rgba [4]int64
140
+			var sum int64
141
+			for i := 0; i < filterLength; i++ {
142
+				xx := start + i
143
+				if xx < oldBounds.Min.X {
144
+					xx = oldBounds.Min.X
145
+				} else if xx >= oldBounds.Max.X {
146
+					xx = oldBounds.Max.X - 1
147
+				}
148
+
149
+				coeff := coeffs[(y-newBounds.Min.Y)*filterLength+i]
150
+				offset := (xx - oldBounds.Min.X) * 8
151
+				rgba[0] += int64(coeff) * int64(uint16(row[offset+0])<<8|uint16(row[offset+1]))
152
+				rgba[1] += int64(coeff) * int64(uint16(row[offset+2])<<8|uint16(row[offset+3]))
153
+				rgba[2] += int64(coeff) * int64(uint16(row[offset+4])<<8|uint16(row[offset+5]))
154
+				rgba[3] += int64(coeff) * int64(uint16(row[offset+6])<<8|uint16(row[offset+7]))
155
+				sum += int64(coeff)
156
+			}
157
+
158
+			offset := (y-newBounds.Min.Y)*out.Stride + (x-newBounds.Min.X)*8
159
+			value := clampUint16(rgba[0] / sum)
160
+			out.Pix[offset+0] = uint8(value >> 8)
161
+			out.Pix[offset+1] = uint8(value)
162
+			value = clampUint16(rgba[1] / sum)
163
+			out.Pix[offset+2] = uint8(value >> 8)
164
+			out.Pix[offset+3] = uint8(value)
165
+			value = clampUint16(rgba[2] / sum)
166
+			out.Pix[offset+4] = uint8(value >> 8)
167
+			out.Pix[offset+5] = uint8(value)
168
+			value = clampUint16(rgba[3] / sum)
169
+			out.Pix[offset+6] = uint8(value >> 8)
170
+			out.Pix[offset+7] = uint8(value)
171
+		}
172
+	}
104 173
 }
105 174
 
106
-func (c *gray16Converter) at(x, y int, result *colorArray) {
107
-	i := c.src.PixOffset(replicateBorder(x, y, c.src.Rect))
108
-	g := float32(uint16(c.src.Pix[i+0])<<8 | uint16(c.src.Pix[i+1]))
109
-	result[0] = g
110
-	result[1] = g
111
-	result[2] = g
112
-	result[3] = float32(0xffff)
113
-	return
175
+func resizeGray(in *image.Gray, out *image.Gray, scale float64, coeffs []int16, filterLength int) {
176
+	oldBounds := in.Bounds()
177
+	newBounds := out.Bounds()
178
+
179
+	for x := newBounds.Min.X; x < newBounds.Max.X; x++ {
180
+		row := in.Pix[(x-oldBounds.Min.Y)*in.Stride:]
181
+		for y := newBounds.Min.Y; y < newBounds.Max.Y; y++ {
182
+			interpX := scale*(float64(y)+0.5) + float64(oldBounds.Min.X)
183
+			start := int(interpX) - filterLength/2 + 1
184
+
185
+			var gray int32
186
+			var sum int32
187
+			for i := 0; i < filterLength; i++ {
188
+				xx := start + i
189
+				if xx < oldBounds.Min.X {
190
+					xx = oldBounds.Min.X
191
+				} else if xx >= oldBounds.Max.X {
192
+					xx = oldBounds.Max.X - 1
193
+				}
194
+
195
+				coeff := coeffs[(y-newBounds.Min.Y)*filterLength+i]
196
+				offset := (xx - oldBounds.Min.X)
197
+				gray += int32(coeff) * int32(row[offset])
198
+				sum += int32(coeff)
199
+			}
200
+
201
+			offset := (y-newBounds.Min.Y)*out.Stride + (x - newBounds.Min.X)
202
+			out.Pix[offset] = clampUint8(gray / sum)
203
+		}
204
+	}
114 205
 }
115 206
 
116
-type ycbcrConverter struct {
117
-	src *image.YCbCr
207
+func resizeGray16(in *image.Gray16, out *image.Gray16, scale float64, coeffs []int32, filterLength int) {
208
+	oldBounds := in.Bounds()
209
+	newBounds := out.Bounds()
210
+
211
+	for x := newBounds.Min.X; x < newBounds.Max.X; x++ {
212
+		row := in.Pix[(x-oldBounds.Min.Y)*in.Stride:]
213
+		for y := newBounds.Min.Y; y < newBounds.Max.Y; y++ {
214
+			interpX := scale*(float64(y)+0.5) + float64(oldBounds.Min.X)
215
+			start := int(interpX) - filterLength/2 + 1
216
+
217
+			var gray int64
218
+			var sum int64
219
+			for i := 0; i < filterLength; i++ {
220
+				xx := start + i
221
+				if xx < oldBounds.Min.X {
222
+					xx = oldBounds.Min.X
223
+				} else if xx >= oldBounds.Max.X {
224
+					xx = oldBounds.Max.X - 1
225
+				}
226
+
227
+				coeff := coeffs[(y-newBounds.Min.Y)*filterLength+i]
228
+				offset := (xx - oldBounds.Min.X) * 2
229
+				gray += int64(coeff) * int64(uint16(row[offset+0])<<8|uint16(row[offset+1]))
230
+				sum += int64(coeff)
231
+			}
232
+
233
+			offset := (y-newBounds.Min.Y)*out.Stride + (x-newBounds.Min.X)*2
234
+			value := clampUint16(gray / sum)
235
+			out.Pix[offset+0] = uint8(value >> 8)
236
+			out.Pix[offset+1] = uint8(value)
237
+		}
238
+	}
118 239
 }
119 240
 
120
-func (c *ycbcrConverter) at(x, y int, result *colorArray) {
121
-	xx, yy := replicateBorder(x, y, c.src.Rect)
122
-	yi := c.src.YOffset(xx, yy)
123
-	ci := c.src.COffset(xx, yy)
124
-	r, g, b := color.YCbCrToRGB(c.src.Y[yi], c.src.Cb[ci], c.src.Cr[ci])
125
-	result[0] = float32(uint16(r) * 0x101)
126
-	result[1] = float32(uint16(g) * 0x101)
127
-	result[2] = float32(uint16(b) * 0x101)
128
-	result[3] = float32(0xffff)
129
-	return
241
+func convertYCbCrToRGBA(in *image.YCbCr) *image.RGBA {
242
+	out := image.NewRGBA(in.Bounds())
243
+	for y := 0; y < out.Bounds().Dy(); y++ {
244
+		for x := 0; x < out.Bounds().Dx(); x++ {
245
+			p := out.Pix[y*out.Stride+4*x:]
246
+			yi := in.YOffset(x, y)
247
+			ci := in.COffset(x, y)
248
+			r, g, b := color.YCbCrToRGB(in.Y[yi], in.Cb[ci], in.Cr[ci])
249
+			p[0] = r
250
+			p[1] = g
251
+			p[2] = b
252
+			p[3] = 0xff
253
+		}
254
+	}
255
+	return out
130 256
 }

+ 65
- 187
filters.go View File

@@ -17,222 +17,100 @@ THIS SOFTWARE.
17 17
 package resize
18 18
 
19 19
 import (
20
-	"image"
21
-	"image/color"
22 20
 	"math"
23 21
 )
24 22
 
25
-// restrict an input float32 to the range of uint16 values
26
-func clampToUint16(x float32) (y uint16) {
27
-	y = uint16(x)
28
-	if x < 0 {
29
-		y = 0
30
-	} else if x > float32(0xfffe) {
31
-		// "else if x > float32(0xffff)" will cause overflows!
32
-		y = 0xffff
23
+func nearest(in float64) float64 {
24
+	if in >= -0.5 && in < 0.5 {
25
+		return 1
33 26
 	}
34
-
35
-	return
27
+	return 0
36 28
 }
37 29
 
38
-// describe a resampling filter
39
-type filterModel struct {
40
-	// resampling is done by convolution with a (scaled) kernel
41
-	kernel func(float32) float32
42
-
43
-	// instead of blurring an image before downscaling to avoid aliasing,
44
-	// the filter is scaled by a factor which leads to a similar effect
45
-	factorInv float32
46
-
47
-	// for optimized access to image points
48
-	converter
49
-
50
-	// temporary used by Interpolate
51
-	tempRow []colorArray
52
-
53
-	kernelWeight []float32
54
-	weightSum    float32
55
-}
56
-
57
-func (f *filterModel) SetKernelWeights(u float32) {
58
-	uf := int(u) - len(f.tempRow)/2 + 1
59
-	u -= float32(uf)
60
-	f.weightSum = 0
61
-
62
-	for j := range f.tempRow {
63
-		f.kernelWeight[j] = f.kernel((u - float32(j)) * f.factorInv)
64
-		f.weightSum += f.kernelWeight[j]
30
+func linear(in float64) float64 {
31
+	in = math.Abs(in)
32
+	if in <= 1 {
33
+		return 1 - in
65 34
 	}
35
+	return 0
66 36
 }
67 37
 
68
-func (f *filterModel) convolution1d() (c colorArray) {
69
-	for j := range f.tempRow {
70
-		for i := range c {
71
-			c[i] += f.tempRow[j][i] * f.kernelWeight[j]
72
-		}
38
+func cubic(in float64) float64 {
39
+	in = math.Abs(in)
40
+	if in <= 1 {
41
+		return in*in*(1.5*in-2.5) + 1.0
73 42
 	}
74
-
75
-	// normalize values
76
-	for i := range c {
77
-		c[i] = c[i] / f.weightSum
43
+	if in <= 2 {
44
+		return in*(in*(2.5-0.5*in)-4.0) + 2.0
78 45
 	}
79
-
80
-	return
46
+	return 0
81 47
 }
82 48
 
83
-func (f *filterModel) Interpolate(u float32, y int) color.RGBA64 {
84
-	uf := int(u) - len(f.tempRow)/2 + 1
85
-	u -= float32(uf)
86
-
87
-	for i := range f.tempRow {
88
-		f.at(uf+i, y, &f.tempRow[i])
49
+func mitchellnetravali(in float64) float64 {
50
+	in = math.Abs(in)
51
+	if in <= 1 {
52
+		return (7.0*in*in*in - 12.0*in*in + 5.33333333333) * 0.16666666666
89 53
 	}
90
-
91
-	c := f.convolution1d()
92
-	return color.RGBA64{
93
-		clampToUint16(c[0]),
94
-		clampToUint16(c[1]),
95
-		clampToUint16(c[2]),
96
-		clampToUint16(c[3]),
54
+	if in <= 2 {
55
+		return (-2.33333333333*in*in*in + 12.0*in*in - 20.0*in + 10.6666666667) * 0.16666666666
97 56
 	}
57
+	return 0
98 58
 }
99 59
 
100
-// createFilter tries to find an optimized converter for the given input image
101
-// and initializes all filterModel members to their defaults
102
-func createFilter(img image.Image, factor float32, size int, kernel func(float32) float32) (f Filter) {
103
-	sizeX := size * (int(math.Ceil(float64(factor))))
104
-
105
-	switch img.(type) {
106
-	default:
107
-		f = &filterModel{
108
-			kernel, 1. / factor,
109
-			&genericConverter{img},
110
-			make([]colorArray, sizeX),
111
-			make([]float32, sizeX),
112
-			0,
113
-		}
114
-	case *image.RGBA:
115
-		f = &filterModel{
116
-			kernel, 1. / factor,
117
-			&rgbaConverter{img.(*image.RGBA)},
118
-			make([]colorArray, sizeX),
119
-			make([]float32, sizeX),
120
-			0,
121
-		}
122
-	case *image.RGBA64:
123
-		f = &filterModel{
124
-			kernel, 1. / factor,
125
-			&rgba64Converter{img.(*image.RGBA64)},
126
-			make([]colorArray, sizeX),
127
-			make([]float32, sizeX),
128
-			0,
129
-		}
130
-	case *image.Gray:
131
-		f = &filterModel{
132
-			kernel, 1. / factor,
133
-			&grayConverter{img.(*image.Gray)},
134
-			make([]colorArray, sizeX),
135
-			make([]float32, sizeX),
136
-			0,
137
-		}
138
-	case *image.Gray16:
139
-		f = &filterModel{
140
-			kernel, 1. / factor,
141
-			&gray16Converter{img.(*image.Gray16)},
142
-			make([]colorArray, sizeX),
143
-			make([]float32, sizeX),
144
-			0,
145
-		}
146
-	case *image.YCbCr:
147
-		f = &filterModel{
148
-			kernel, 1. / factor,
149
-			&ycbcrConverter{img.(*image.YCbCr)},
150
-			make([]colorArray, sizeX),
151
-			make([]float32, sizeX),
152
-			0,
153
-		}
60
+func sinc(x float64) float64 {
61
+	x = math.Abs(x) * math.Pi
62
+	if x >= 1.220703e-4 {
63
+		return math.Sin(x) / x
154 64
 	}
155
-
156
-	return
157
-}
158
-
159
-// Nearest-neighbor interpolation
160
-func NearestNeighbor(img image.Image, factor float32) Filter {
161
-	return createFilter(img, factor, 2, func(x float32) (y float32) {
162
-		if x >= -0.5 && x < 0.5 {
163
-			y = 1
164
-		} else {
165
-			y = 0
166
-		}
167
-
168
-		return
169
-	})
170
-}
171
-
172
-// Bilinear interpolation
173
-func Bilinear(img image.Image, factor float32) Filter {
174
-	return createFilter(img, factor, 2, func(x float32) (y float32) {
175
-		absX := float32(math.Abs(float64(x)))
176
-		if absX <= 1 {
177
-			y = 1 - absX
178
-		} else {
179
-			y = 0
180
-		}
181
-
182
-		return
183
-	})
65
+	return 1
184 66
 }
185 67
 
186
-// Bicubic interpolation (with cubic hermite spline)
187
-func Bicubic(img image.Image, factor float32) Filter {
188
-	return createFilter(img, factor, 4, splineKernel(0, 0.5))
68
+func lanczos2(in float64) float64 {
69
+	if in > -2 && in < 2 {
70
+		return sinc(in) * sinc(in*0.5)
71
+	}
72
+	return 0
189 73
 }
190 74
 
191
-// Mitchell-Netravali interpolation
192
-func MitchellNetravali(img image.Image, factor float32) Filter {
193
-	return createFilter(img, factor, 4, splineKernel(1.0/3.0, 1.0/3.0))
75
+func lanczos3(in float64) float64 {
76
+	if in > -3 && in < 3 {
77
+		return sinc(in) * sinc(in*0.3333333333333333)
78
+	}
79
+	return 0
194 80
 }
195 81
 
196
-func splineKernel(B, C float32) func(float32) float32 {
197
-	factorA := 2.0 - 1.5*B - C
198
-	factorB := -3.0 + 2.0*B + C
199
-	factorC := 1.0 - 1.0/3.0*B
200
-	factorD := -B/6.0 - C
201
-	factorE := B + 5.0*C
202
-	factorF := -2.0*B - 8.0*C
203
-	factorG := 4.0/3.0*B + 4.0*C
204
-	return func(x float32) (y float32) {
205
-		absX := float32(math.Abs(float64(x)))
206
-		if absX <= 1 {
207
-			y = absX*absX*(factorA*absX+factorB) + factorC
208
-		} else if absX <= 2 {
209
-			y = absX*(absX*(absX*factorD+factorE)+factorF) + factorG
210
-		} else {
211
-			y = 0
82
+// range [-256,256]
83
+func createWeights8(dy, minx, filterLength int, blur, scale float64, kernel func(float64) float64) ([]int16, int) {
84
+	filterLength = filterLength * int(math.Max(math.Ceil(blur*scale), 1))
85
+	filterFactor := math.Min(1./(blur*scale), 1)
86
+
87
+	coeffs := make([]int16, dy*filterLength)
88
+	for y := 0; y < dy; y++ {
89
+		interpX := scale*(float64(y)+0.5) + float64(minx)
90
+		start := int(interpX) - filterLength/2 + 1
91
+		for i := 0; i < filterLength; i++ {
92
+			in := (interpX - float64(start) - float64(i)) * filterFactor
93
+			coeffs[y*filterLength+i] = int16(kernel(in) * 256)
212 94
 		}
213
-
214
-		return
215 95
 	}
96
+
97
+	return coeffs, filterLength
216 98
 }
217 99
 
218
-func lanczosKernel(a uint) func(float32) float32 {
219
-	return func(x float32) (y float32) {
220
-		if x > -float32(a) && x < float32(a) {
221
-			y = float32(Sinc(float64(x))) * float32(Sinc(float64(x/float32(a))))
222
-		} else {
223
-			y = 0
100
+// range [-65536,65536]
101
+func createWeights16(dy, minx, filterLength int, blur, scale float64, kernel func(float64) float64) ([]int32, int) {
102
+	filterLength = filterLength * int(math.Max(math.Ceil(blur*scale), 1))
103
+	filterFactor := math.Min(1./(blur*scale), 1)
104
+
105
+	coeffs := make([]int32, dy*filterLength)
106
+	for y := 0; y < dy; y++ {
107
+		interpX := scale*(float64(y)+0.5) + float64(minx)
108
+		start := int(interpX) - filterLength/2 + 1
109
+		for i := 0; i < filterLength; i++ {
110
+			in := (interpX - float64(start) - float64(i)) * filterFactor
111
+			coeffs[y*filterLength+i] = int32(kernel(in) * 65536)
224 112
 		}
225
-
226
-		return
227 113
 	}
228
-}
229
-
230
-// Lanczos interpolation (a=2)
231
-func Lanczos2(img image.Image, factor float32) Filter {
232
-	return createFilter(img, factor, 4, lanczosKernel(2))
233
-}
234 114
 
235
-// Lanczos interpolation (a=3)
236
-func Lanczos3(img image.Image, factor float32) Filter {
237
-	return createFilter(img, factor, 6, lanczosKernel(3))
115
+	return coeffs, filterLength
238 116
 }

+ 232
- 81
resize.go View File

@@ -26,124 +26,275 @@ package resize
26 26
 
27 27
 import (
28 28
 	"image"
29
-	"image/color"
30 29
 	"runtime"
30
+	"sync"
31 31
 )
32 32
 
33
-// Filter can interpolate at points (x,y)
34
-type Filter interface {
35
-	SetKernelWeights(u float32)
36
-	Interpolate(u float32, y int) color.RGBA64
33
+// An InterpolationFunction provides the parameters that describe an
34
+// interpolation kernel. It returns the number of samples to take
35
+// and the kernel function to use for sampling.
36
+type InterpolationFunction func() (int, func(float64) float64)
37
+
38
+// Nearest-neighbor interpolation
39
+func NearestNeighbor() (int, func(float64) float64) {
40
+	return 2, nearest
41
+}
42
+
43
+// Bilinear interpolation
44
+func Bilinear() (int, func(float64) float64) {
45
+	return 2, linear
46
+}
47
+
48
+// Bicubic interpolation (with cubic hermite spline)
49
+func Bicubic() (int, func(float64) float64) {
50
+	return 4, cubic
51
+}
52
+
53
+// Mitchell-Netravali interpolation
54
+func MitchellNetravali() (int, func(float64) float64) {
55
+	return 4, mitchellnetravali
56
+}
57
+
58
+// Lanczos interpolation (a=2)
59
+func Lanczos2() (int, func(float64) float64) {
60
+	return 4, lanczos2
61
+}
62
+
63
+// Lanczos interpolation (a=3)
64
+func Lanczos3() (int, func(float64) float64) {
65
+	return 6, lanczos3
37 66
 }
38 67
 
39
-// InterpolationFunction return a Filter implementation
40
-// that operates on an image. Two factors
41
-// allow to scale the filter kernels in x- and y-direction
42
-// to prevent moire patterns.
43
-type InterpolationFunction func(image.Image, float32) Filter
68
+// values <1 will sharpen the image
69
+var blur = 1.0
44 70
 
45
-// Resize an image to new width and height using the interpolation function interp.
71
+// Resize scales an image to new width and height using the interpolation function interp.
46 72
 // A new image with the given dimensions will be returned.
47 73
 // If one of the parameters width or height is set to 0, its size will be calculated so that
48 74
 // the aspect ratio is that of the originating image.
49 75
 // The resizing algorithm uses channels for parallel computation.
50 76
 func Resize(width, height uint, img image.Image, interp InterpolationFunction) image.Image {
51
-	oldBounds := img.Bounds()
52
-	oldWidth := float32(oldBounds.Dx())
53
-	oldHeight := float32(oldBounds.Dy())
54
-	scaleX, scaleY := calcFactors(width, height, oldWidth, oldHeight)
55
-
56
-	tempImg := image.NewRGBA64(image.Rect(0, 0, oldBounds.Dy(), int(0.7+oldWidth/scaleX)))
57
-	b := tempImg.Bounds()
58
-	adjust := 0.5 * ((oldWidth-1.0)/scaleX - float32(b.Dy()-1))
59
-
60
-	n := numJobs(b.Dy())
61
-	c := make(chan int, n)
62
-	for i := 0; i < n; i++ {
63
-		slice := image.Rect(b.Min.X, b.Min.Y+i*(b.Dy())/n, b.Max.X, b.Min.Y+(i+1)*(b.Dy())/n)
64
-		go resizeSlice(img, tempImg, interp, scaleX, adjust, float32(oldBounds.Min.X), oldBounds.Min.Y, slice, c)
77
+	scaleX, scaleY := calcFactors(width, height, float64(img.Bounds().Dx()), float64(img.Bounds().Dy()))
78
+	if width == 0 {
79
+		width = uint(0.7 + float64(img.Bounds().Dx())/scaleX)
65 80
 	}
66
-	for i := 0; i < n; i++ {
67
-		<-c
81
+	if height == 0 {
82
+		height = uint(0.7 + float64(img.Bounds().Dy())/scaleY)
68 83
 	}
69 84
 
70
-	resultImg := image.NewRGBA64(image.Rect(0, 0, int(0.7+oldWidth/scaleX), int(0.7+oldHeight/scaleY)))
71
-	b = resultImg.Bounds()
72
-	adjust = 0.5 * ((oldHeight-1.0)/scaleY - float32(b.Dy()-1))
85
+	taps, kernel := interp()
86
+	cpus := runtime.NumCPU()
87
+	wg := sync.WaitGroup{}
73 88
 
74
-	for i := 0; i < n; i++ {
75
-		slice := image.Rect(b.Min.X, b.Min.Y+i*(b.Dy())/n, b.Max.X, b.Min.Y+(i+1)*(b.Dy())/n)
76
-		go resizeSlice(tempImg, resultImg, interp, scaleY, adjust, 0, 0, slice, c)
77
-	}
78
-	for i := 0; i < n; i++ {
79
-		<-c
80
-	}
89
+	// Generic access to image.Image is slow in tight loops.
90
+	// The optimal access has to be determined from the concrete image type.
91
+	switch input := img.(type) {
92
+	case *image.RGBA:
93
+		// 8-bit precision
94
+		temp := image.NewRGBA(image.Rect(0, 0, input.Bounds().Dy(), int(width)))
95
+		result := image.NewRGBA(image.Rect(0, 0, int(width), int(height)))
81 96
 
82
-	return resultImg
83
-}
97
+		// horizontal filter, results in transposed temporary image
98
+		coeffs, filterLength := createWeights8(temp.Bounds().Dy(), input.Bounds().Min.X, taps, blur, scaleX, kernel)
99
+		wg.Add(cpus)
100
+		for i := 0; i < cpus; i++ {
101
+			slice := makeSlice(temp, i, cpus).(*image.RGBA)
102
+			go func() {
103
+				defer wg.Done()
104
+				resizeRGBA(input, slice, scaleX, coeffs, filterLength)
105
+			}()
106
+		}
107
+		wg.Wait()
84 108
 
85
-// Resize a rectangle image slice
86
-func resizeSlice(input image.Image, output *image.RGBA64, interp InterpolationFunction, scale, adjust, xoffset float32, yoffset int, slice image.Rectangle, c chan int) {
87
-	filter := interp(input, float32(clampFactor(scale)))
88
-	var u float32
89
-	var color color.RGBA64
90
-	for y := slice.Min.Y; y < slice.Max.Y; y++ {
91
-		u = scale*(float32(y)+adjust) + xoffset
92
-		filter.SetKernelWeights(u)
93
-		for x := slice.Min.X; x < slice.Max.X; x++ {
94
-			color = filter.Interpolate(u, x+yoffset)
95
-			i := output.PixOffset(x, y)
96
-			output.Pix[i+0] = uint8(color.R >> 8)
97
-			output.Pix[i+1] = uint8(color.R)
98
-			output.Pix[i+2] = uint8(color.G >> 8)
99
-			output.Pix[i+3] = uint8(color.G)
100
-			output.Pix[i+4] = uint8(color.B >> 8)
101
-			output.Pix[i+5] = uint8(color.B)
102
-			output.Pix[i+6] = uint8(color.A >> 8)
103
-			output.Pix[i+7] = uint8(color.A)
109
+		// horizontal filter on transposed image, result is not transposed
110
+		coeffs, filterLength = createWeights8(result.Bounds().Dy(), temp.Bounds().Min.X, taps, blur, scaleY, kernel)
111
+		wg.Add(cpus)
112
+		for i := 0; i < cpus; i++ {
113
+			slice := makeSlice(result, i, cpus).(*image.RGBA)
114
+			go func() {
115
+				defer wg.Done()
116
+				resizeRGBA(temp, slice, scaleY, coeffs, filterLength)
117
+			}()
104 118
 		}
105
-	}
119
+		wg.Wait()
120
+		return result
121
+	case *image.YCbCr:
122
+		// 8-bit precision
123
+		// accessing the YCbCr arrays in a tight loop is slow.
124
+		// converting the image before filtering will improve performance.
125
+		inputAsRGBA := convertYCbCrToRGBA(input)
126
+		temp := image.NewRGBA(image.Rect(0, 0, input.Bounds().Dy(), int(width)))
127
+		result := image.NewRGBA(image.Rect(0, 0, int(width), int(height)))
128
+
129
+		// horizontal filter, results in transposed temporary image
130
+		coeffs, filterLength := createWeights8(temp.Bounds().Dy(), input.Bounds().Min.X, taps, blur, scaleX, kernel)
131
+		wg.Add(cpus)
132
+		for i := 0; i < cpus; i++ {
133
+			slice := makeSlice(temp, i, cpus).(*image.RGBA)
134
+			go func() {
135
+				defer wg.Done()
136
+				resizeRGBA(inputAsRGBA, slice, scaleX, coeffs, filterLength)
137
+			}()
138
+		}
139
+		wg.Wait()
140
+
141
+		// horizontal filter on transposed image, result is not transposed
142
+		coeffs, filterLength = createWeights8(result.Bounds().Dy(), temp.Bounds().Min.X, taps, blur, scaleY, kernel)
143
+		wg.Add(cpus)
144
+		for i := 0; i < cpus; i++ {
145
+			slice := makeSlice(result, i, cpus).(*image.RGBA)
146
+			go func() {
147
+				defer wg.Done()
148
+				resizeRGBA(temp, slice, scaleY, coeffs, filterLength)
149
+			}()
150
+		}
151
+		wg.Wait()
152
+		return result
153
+	case *image.RGBA64:
154
+		// 16-bit precision
155
+		temp := image.NewRGBA64(image.Rect(0, 0, input.Bounds().Dy(), int(width)))
156
+		result := image.NewRGBA64(image.Rect(0, 0, int(width), int(height)))
157
+
158
+		// horizontal filter, results in transposed temporary image
159
+		coeffs, filterLength := createWeights16(temp.Bounds().Dy(), input.Bounds().Min.X, taps, blur, scaleX, kernel)
160
+		wg.Add(cpus)
161
+		for i := 0; i < cpus; i++ {
162
+			slice := makeSlice(temp, i, cpus).(*image.RGBA64)
163
+			go func() {
164
+				defer wg.Done()
165
+				resizeRGBA64(input, slice, scaleX, coeffs, filterLength)
166
+			}()
167
+		}
168
+		wg.Wait()
169
+
170
+		// horizontal filter on transposed image, result is not transposed
171
+		coeffs, filterLength = createWeights16(result.Bounds().Dy(), temp.Bounds().Min.X, taps, blur, scaleY, kernel)
172
+		wg.Add(cpus)
173
+		for i := 0; i < cpus; i++ {
174
+			slice := makeSlice(result, i, cpus).(*image.RGBA64)
175
+			go func() {
176
+				defer wg.Done()
177
+				resizeGeneric(temp, slice, scaleY, coeffs, filterLength)
178
+			}()
179
+		}
180
+		wg.Wait()
181
+		return result
182
+	case *image.Gray:
183
+		// 8-bit precision
184
+		temp := image.NewGray(image.Rect(0, 0, input.Bounds().Dy(), int(width)))
185
+		result := image.NewGray(image.Rect(0, 0, int(width), int(height)))
186
+
187
+		// horizontal filter, results in transposed temporary image
188
+		coeffs, filterLength := createWeights8(temp.Bounds().Dy(), input.Bounds().Min.X, taps, blur, scaleX, kernel)
189
+		wg.Add(cpus)
190
+		for i := 0; i < cpus; i++ {
191
+			slice := makeSlice(temp, i, cpus).(*image.Gray)
192
+			go func() {
193
+				defer wg.Done()
194
+				resizeGray(input, slice, scaleX, coeffs, filterLength)
195
+			}()
196
+		}
197
+		wg.Wait()
198
+
199
+		// horizontal filter on transposed image, result is not transposed
200
+		coeffs, filterLength = createWeights8(result.Bounds().Dy(), temp.Bounds().Min.X, taps, blur, scaleY, kernel)
201
+		wg.Add(cpus)
202
+		for i := 0; i < cpus; i++ {
203
+			slice := makeSlice(result, i, cpus).(*image.Gray)
204
+			go func() {
205
+				defer wg.Done()
206
+				resizeGray(temp, slice, scaleY, coeffs, filterLength)
207
+			}()
208
+		}
209
+		wg.Wait()
210
+		return result
211
+	case *image.Gray16:
212
+		// 16-bit precision
213
+		temp := image.NewGray16(image.Rect(0, 0, input.Bounds().Dy(), int(width)))
214
+		result := image.NewGray16(image.Rect(0, 0, int(width), int(height)))
215
+
216
+		// horizontal filter, results in transposed temporary image
217
+		coeffs, filterLength := createWeights16(temp.Bounds().Dy(), input.Bounds().Min.X, taps, blur, scaleX, kernel)
218
+		wg.Add(cpus)
219
+		for i := 0; i < cpus; i++ {
220
+			slice := makeSlice(temp, i, cpus).(*image.Gray16)
221
+			go func() {
222
+				defer wg.Done()
223
+				resizeGray16(input, slice, scaleX, coeffs, filterLength)
224
+			}()
225
+		}
226
+		wg.Wait()
227
+
228
+		// horizontal filter on transposed image, result is not transposed
229
+		coeffs, filterLength = createWeights16(result.Bounds().Dy(), temp.Bounds().Min.X, taps, blur, scaleY, kernel)
230
+		wg.Add(cpus)
231
+		for i := 0; i < cpus; i++ {
232
+			slice := makeSlice(result, i, cpus).(*image.Gray16)
233
+			go func() {
234
+				defer wg.Done()
235
+				resizeGray16(temp, slice, scaleY, coeffs, filterLength)
236
+			}()
237
+		}
238
+		wg.Wait()
239
+		return result
240
+	default:
241
+		// 16-bit precision
242
+		temp := image.NewRGBA64(image.Rect(0, 0, img.Bounds().Dy(), int(width)))
243
+		result := image.NewRGBA64(image.Rect(0, 0, int(width), int(height)))
244
+
245
+		// horizontal filter, results in transposed temporary image
246
+		coeffs, filterLength := createWeights16(temp.Bounds().Dy(), img.Bounds().Min.X, taps, blur, scaleX, kernel)
247
+		wg.Add(cpus)
248
+		for i := 0; i < cpus; i++ {
249
+			slice := makeSlice(temp, i, cpus).(*image.RGBA64)
250
+			go func() {
251
+				defer wg.Done()
252
+				resizeGeneric(img, slice, scaleX, coeffs, filterLength)
253
+			}()
254
+		}
255
+		wg.Wait()
106 256
 
107
-	c <- 1
257
+		// horizontal filter on transposed image, result is not transposed
258
+		coeffs, filterLength = createWeights16(result.Bounds().Dy(), temp.Bounds().Min.X, taps, blur, scaleY, kernel)
259
+		wg.Add(cpus)
260
+		for i := 0; i < cpus; i++ {
261
+			slice := makeSlice(result, i, cpus).(*image.RGBA64)
262
+			go func() {
263
+				defer wg.Done()
264
+				resizeRGBA64(temp, slice, scaleY, coeffs, filterLength)
265
+			}()
266
+		}
267
+		wg.Wait()
268
+		return result
269
+	}
108 270
 }
109 271
 
110
-// Calculate scaling factors using old and new image dimensions.
111
-func calcFactors(width, height uint, oldWidth, oldHeight float32) (scaleX, scaleY float32) {
272
+// Calculates scaling factors using old and new image dimensions.
273
+func calcFactors(width, height uint, oldWidth, oldHeight float64) (scaleX, scaleY float64) {
112 274
 	if width == 0 {
113 275
 		if height == 0 {
114 276
 			scaleX = 1.0
115 277
 			scaleY = 1.0
116 278
 		} else {
117
-			scaleY = oldHeight / float32(height)
279
+			scaleY = oldHeight / float64(height)
118 280
 			scaleX = scaleY
119 281
 		}
120 282
 	} else {
121
-		scaleX = oldWidth / float32(width)
283
+		scaleX = oldWidth / float64(width)
122 284
 		if height == 0 {
123 285
 			scaleY = scaleX
124 286
 		} else {
125
-			scaleY = oldHeight / float32(height)
287
+			scaleY = oldHeight / float64(height)
126 288
 		}
127 289
 	}
128 290
 	return
129 291
 }
130 292
 
131
-// Set filter scaling factor to avoid moire patterns.
132
-// This is only useful in case of downscaling (factor>1).
133
-func clampFactor(factor float32) float32 {
134
-	if factor < 1 {
135
-		factor = 1
136
-	}
137
-	return factor
293
+type imageWithSubImage interface {
294
+	image.Image
295
+	SubImage(image.Rectangle) image.Image
138 296
 }
139 297
 
140
-// Set number of parallel jobs
141
-// but prevent resize from doing too much work
142
-// if #CPUs > width
143
-func numJobs(d int) (n int) {
144
-	n = runtime.NumCPU()
145
-	if n > d {
146
-		n = d
147
-	}
148
-	return
298
+func makeSlice(img imageWithSubImage, i, n int) image.Image {
299
+	return img.SubImage(image.Rect(img.Bounds().Min.X, img.Bounds().Min.Y+i*img.Bounds().Dy()/n, img.Bounds().Max.X, img.Bounds().Min.Y+(i+1)*img.Bounds().Dy()/n))
149 300
 }

+ 18
- 7
resize_test.go View File

@@ -14,13 +14,6 @@ func init() {
14 14
 	img.Set(1, 1, color.White)
15 15
 }
16 16
 
17
-func Test_Nearest(t *testing.T) {
18
-	m := Resize(6, 0, img, NearestNeighbor)
19
-	if m.At(1, 1) == m.At(2, 2) {
20
-		t.Fail()
21
-	}
22
-}
23
-
24 17
 func Test_Param1(t *testing.T) {
25 18
 	m := Resize(0, 0, img, NearestNeighbor)
26 19
 	if m.Bounds() != img.Bounds() {
@@ -53,6 +46,24 @@ func Test_CorrectResize(t *testing.T) {
53 46
 	}
54 47
 }
55 48
 
49
+func Test_SameColor(t *testing.T) {
50
+	img := image.NewRGBA(image.Rect(0, 0, 20, 20))
51
+	for y := img.Bounds().Min.Y; y < img.Bounds().Max.Y; y++ {
52
+		for x := img.Bounds().Min.X; x < img.Bounds().Max.X; x++ {
53
+			img.SetRGBA(x, y, color.RGBA{0x80, 0x80, 0x80, 0xFF})
54
+		}
55
+	}
56
+	out := Resize(10, 10, img, Lanczos3)
57
+	for y := out.Bounds().Min.Y; y < out.Bounds().Max.Y; y++ {
58
+		for x := out.Bounds().Min.X; x < out.Bounds().Max.X; x++ {
59
+			color := img.At(x, y).(color.RGBA)
60
+			if color.R != 0x80 || color.G != 0x80 || color.B != 0x80 || color.A != 0xFF {
61
+				t.Fail()
62
+			}
63
+		}
64
+	}
65
+}
66
+
56 67
 func Benchmark_BigResizeLanczos3(b *testing.B) {
57 68
 	var m image.Image
58 69
 	for i := 0; i < b.N; i++ {

+ 0
- 49
sinc.go View File

@@ -1,49 +0,0 @@
1
-/*
2
-Copyright (c) 2012, Jan Schlicht <jan.schlicht@gmail.com>
3
-
4
-Permission to use, copy, modify, and/or distribute this software for any purpose
5
-with or without fee is hereby granted, provided that the above copyright notice
6
-and this permission notice appear in all copies.
7
-
8
-THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
9
-REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
10
-FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
11
-INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
12
-OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
13
-TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
14
-THIS SOFTWARE.
15
-*/
16
-
17
-package resize
18
-
19
-import (
20
-	"math"
21
-)
22
-
23
-var (
24
-	epsilon      = math.Nextafter(1.0, 2.0) - 1.0 // machine epsilon
25
-	taylor2bound = math.Sqrt(epsilon)
26
-	taylorNbound = math.Sqrt(taylor2bound)
27
-)
28
-
29
-// unnormalized sinc function
30
-func Sinc1(x float64) (y float64) {
31
-	if math.Abs(x) >= taylorNbound {
32
-		y = math.Sin(x) / x
33
-	} else {
34
-		y = 1.0
35
-		if math.Abs(x) >= epsilon {
36
-			x2 := x * x
37
-			y -= x2 / 6.0
38
-			if math.Abs(x) >= taylor2bound {
39
-				y += (x2 * x2) / 120.0
40
-			}
41
-		}
42
-	}
43
-	return
44
-}
45
-
46
-// normalized sinc function
47
-func Sinc(x float64) float64 {
48
-	return Sinc1(x * math.Pi)
49
-}

+ 0
- 38
sinc_test.go View File

@@ -1,38 +0,0 @@
1
-package resize
2
-
3
-import (
4
-	"fmt"
5
-	"math"
6
-	"testing"
7
-)
8
-
9
-const limit = 1e-12
10
-
11
-func Test_SincOne(t *testing.T) {
12
-	zero := Sinc(1)
13
-	if zero >= limit {
14
-		t.Error("Sinc(1) != 0")
15
-	}
16
-}
17
-
18
-func Test_SincZero(t *testing.T) {
19
-	one := Sinc(0)
20
-	if math.Abs(one-1) >= limit {
21
-		t.Error("Sinc(0) != 1")
22
-	}
23
-}
24
-
25
-func Test_SincDotOne(t *testing.T) {
26
-	res := Sinc(0.1)
27
-	if math.Abs(res-0.983631643083466) >= limit {
28
-		t.Error("Sinc(0.1) wrong")
29
-	}
30
-}
31
-
32
-func Test_SincNearZero(t *testing.T) {
33
-	res := Sinc(0.000001)
34
-	if math.Abs(res-0.9999999999983551) >= limit {
35
-		fmt.Println(res)
36
-		t.Error("Sinc near zero not stable")
37
-	}
38
-}