convolution.go

  1package gift
  2
  3import (
  4	"image"
  5	"image/draw"
  6	"math"
  7)
  8
  9type uweight struct {
 10	u      int
 11	weight float32
 12}
 13
 14type uvweight struct {
 15	u      int
 16	v      int
 17	weight float32
 18}
 19
 20func prepareConvolutionWeights(kernel []float32, normalize bool) (int, []uvweight) {
 21	size := int(math.Sqrt(float64(len(kernel))))
 22	if size%2 == 0 {
 23		size--
 24	}
 25	if size < 1 {
 26		return 0, []uvweight{}
 27	}
 28	center := size / 2
 29
 30	weights := []uvweight{}
 31	for i := 0; i < size; i++ {
 32		for j := 0; j < size; j++ {
 33			k := j*size + i
 34			w := float32(0)
 35			if k < len(kernel) {
 36				w = kernel[k]
 37			}
 38			if w != 0 {
 39				weights = append(weights, uvweight{u: i - center, v: j - center, weight: w})
 40			}
 41		}
 42	}
 43
 44	if !normalize {
 45		return size, weights
 46	}
 47
 48	var sum, sumpositive float32
 49	for _, w := range weights {
 50		sum += w.weight
 51		if w.weight > 0 {
 52			sumpositive += w.weight
 53		}
 54	}
 55
 56	var div float32
 57	if sum != 0 {
 58		div = sum
 59	} else if sumpositive != 0 {
 60		div = sumpositive
 61	} else {
 62		return size, weights
 63	}
 64
 65	for i := 0; i < len(weights); i++ {
 66		weights[i].weight /= div
 67	}
 68
 69	return size, weights
 70}
 71
 72type convolutionFilter struct {
 73	kernel    []float32
 74	normalize bool
 75	alpha     bool
 76	abs       bool
 77	delta     float32
 78}
 79
 80func (p *convolutionFilter) Bounds(srcBounds image.Rectangle) (dstBounds image.Rectangle) {
 81	dstBounds = image.Rect(0, 0, srcBounds.Dx(), srcBounds.Dy())
 82	return
 83}
 84
 85func (p *convolutionFilter) Draw(dst draw.Image, src image.Image, options *Options) {
 86	if options == nil {
 87		options = &defaultOptions
 88	}
 89
 90	srcb := src.Bounds()
 91	dstb := dst.Bounds()
 92
 93	if srcb.Dx() <= 0 || srcb.Dy() <= 0 {
 94		return
 95	}
 96
 97	ksize, weights := prepareConvolutionWeights(p.kernel, p.normalize)
 98	kcenter := ksize / 2
 99
100	if ksize < 1 {
101		copyimage(dst, src, options)
102		return
103	}
104
105	pixGetter := newPixelGetter(src)
106	pixSetter := newPixelSetter(dst)
107
108	parallelize(options.Parallelization, srcb.Min.Y, srcb.Max.Y, func(pmin, pmax int) {
109		// init temp rows
110		starty := pmin
111		rows := make([][]pixel, ksize)
112		for i := 0; i < ksize; i++ {
113			rowy := starty + i - kcenter
114			if rowy < srcb.Min.Y {
115				rowy = srcb.Min.Y
116			} else if rowy > srcb.Max.Y-1 {
117				rowy = srcb.Max.Y - 1
118			}
119			row := make([]pixel, srcb.Dx())
120			pixGetter.getPixelRow(rowy, &row)
121			rows[i] = row
122		}
123
124		for y := pmin; y < pmax; y++ {
125			// calculate dst row
126			for x := srcb.Min.X; x < srcb.Max.X; x++ {
127				var r, g, b, a float32
128				for _, w := range weights {
129					wx := x + w.u
130					if wx < srcb.Min.X {
131						wx = srcb.Min.X
132					} else if wx > srcb.Max.X-1 {
133						wx = srcb.Max.X - 1
134					}
135					rowsx := wx - srcb.Min.X
136					rowsy := kcenter + w.v
137
138					px := rows[rowsy][rowsx]
139					r += px.r * w.weight
140					g += px.g * w.weight
141					b += px.b * w.weight
142					if p.alpha {
143						a += px.a * w.weight
144					}
145				}
146				if p.abs {
147					r = absf32(r)
148					g = absf32(g)
149					b = absf32(b)
150					if p.alpha {
151						a = absf32(a)
152					}
153				}
154				if p.delta != 0 {
155					r += p.delta
156					g += p.delta
157					b += p.delta
158					if p.alpha {
159						a += p.delta
160					}
161				}
162				if !p.alpha {
163					a = rows[kcenter][x-srcb.Min.X].a
164				}
165				pixSetter.setPixel(dstb.Min.X+x-srcb.Min.X, dstb.Min.Y+y-srcb.Min.Y, pixel{r, g, b, a})
166			}
167
168			// rotate temp rows
169			if y < pmax-1 {
170				tmprow := rows[0]
171				for i := 0; i < ksize-1; i++ {
172					rows[i] = rows[i+1]
173				}
174				nextrowy := y + ksize/2 + 1
175				if nextrowy > srcb.Max.Y-1 {
176					nextrowy = srcb.Max.Y - 1
177				}
178				pixGetter.getPixelRow(nextrowy, &tmprow)
179				rows[ksize-1] = tmprow
180			}
181		}
182	})
183}
184
185// Convolution creates a filter that applies a square convolution kernel to an image.
186// The length of the kernel slice must be the square of an odd kernel size (e.g. 9 for 3x3 kernel, 25 for 5x5 kernel).
187// Excessive slice members will be ignored.
188// If normalize parameter is true, the kernel will be normalized before applying the filter.
189// If alpha parameter is true, the alpha component of color will be filtered too.
190// If abs parameter is true, absolute values of color components will be taken after doing calculations.
191// If delta parameter is not zero, this value will be added to the filtered pixels.
192//
193// Example:
194//
195//	// Apply the emboss filter to an image.
196//	g := gift.New(
197//		gift.Convolution(
198//			[]float32{
199//				-1, -1, 0,
200//				-1, 1, 1,
201//				0, 1, 1,
202//			},
203//			false, false, false, 0,
204//		),
205//	)
206//	dst := image.NewRGBA(g.Bounds(src.Bounds()))
207//	g.Draw(dst, src)
208//
209func Convolution(kernel []float32, normalize, alpha, abs bool, delta float32) Filter {
210	return &convolutionFilter{
211		kernel:    kernel,
212		normalize: normalize,
213		alpha:     alpha,
214		abs:       abs,
215		delta:     delta,
216	}
217}
218
219// prepare pixel weights using convolution kernel. weights equal to 0 are excluded
220func prepareConvolutionWeights1d(kernel []float32) (int, []uweight) {
221	size := len(kernel)
222	if size%2 == 0 {
223		size--
224	}
225	if size < 1 {
226		return 0, []uweight{}
227	}
228	center := size / 2
229	weights := []uweight{}
230	for i := 0; i < size; i++ {
231		w := float32(0)
232		if i < len(kernel) {
233			w = kernel[i]
234		}
235		if w != 0 {
236			weights = append(weights, uweight{i - center, w})
237		}
238	}
239	return size, weights
240}
241
242// calculate pixels for one line according to weights
243func convolveLine(dstBuf []pixel, srcBuf []pixel, weights []uweight) {
244	max := len(srcBuf) - 1
245	if max < 0 {
246		return
247	}
248	for dstu := 0; dstu < len(srcBuf); dstu++ {
249		var r, g, b, a float32
250		for _, w := range weights {
251			k := dstu + w.u
252			if k < 0 {
253				k = 0
254			} else if k > max {
255				k = max
256			}
257			c := srcBuf[k]
258			wa := c.a * w.weight
259			r += c.r * wa
260			g += c.g * wa
261			b += c.b * wa
262			a += wa
263		}
264		if a != 0 {
265			r /= a
266			g /= a
267			b /= a
268		}
269		dstBuf[dstu] = pixel{r, g, b, a}
270	}
271}
272
273// fast vertical 1d convolution
274func convolve1dv(dst draw.Image, src image.Image, kernel []float32, options *Options) {
275	srcb := src.Bounds()
276	dstb := dst.Bounds()
277	if srcb.Dx() <= 0 || srcb.Dy() <= 0 {
278		return
279	}
280	if kernel == nil || len(kernel) < 1 {
281		copyimage(dst, src, options)
282		return
283	}
284	_, weights := prepareConvolutionWeights1d(kernel)
285	pixGetter := newPixelGetter(src)
286	pixSetter := newPixelSetter(dst)
287	parallelize(options.Parallelization, srcb.Min.X, srcb.Max.X, func(pmin, pmax int) {
288		srcBuf := make([]pixel, srcb.Dy())
289		dstBuf := make([]pixel, srcb.Dy())
290		for x := pmin; x < pmax; x++ {
291			pixGetter.getPixelColumn(x, &srcBuf)
292			convolveLine(dstBuf, srcBuf, weights)
293			pixSetter.setPixelColumn(dstb.Min.X+x-srcb.Min.X, dstBuf)
294		}
295	})
296}
297
298// fast horizontal 1d convolution
299func convolve1dh(dst draw.Image, src image.Image, kernel []float32, options *Options) {
300	srcb := src.Bounds()
301	dstb := dst.Bounds()
302	if srcb.Dx() <= 0 || srcb.Dy() <= 0 {
303		return
304	}
305	if kernel == nil || len(kernel) < 1 {
306		copyimage(dst, src, options)
307		return
308	}
309	_, weights := prepareConvolutionWeights1d(kernel)
310	pixGetter := newPixelGetter(src)
311	pixSetter := newPixelSetter(dst)
312	parallelize(options.Parallelization, srcb.Min.Y, srcb.Max.Y, func(pmin, pmax int) {
313		srcBuf := make([]pixel, srcb.Dx())
314		dstBuf := make([]pixel, srcb.Dx())
315		for y := pmin; y < pmax; y++ {
316			pixGetter.getPixelRow(y, &srcBuf)
317			convolveLine(dstBuf, srcBuf, weights)
318			pixSetter.setPixelRow(dstb.Min.Y+y-srcb.Min.Y, dstBuf)
319		}
320	})
321}
322
323func gaussianBlurKernel(x, sigma float32) float32 {
324	return float32(math.Exp(-float64(x*x)/float64(2*sigma*sigma)) / (float64(sigma) * math.Sqrt(2*math.Pi)))
325}
326
327type gausssianBlurFilter struct {
328	sigma float32
329}
330
331func (p *gausssianBlurFilter) Bounds(srcBounds image.Rectangle) (dstBounds image.Rectangle) {
332	dstBounds = image.Rect(0, 0, srcBounds.Dx(), srcBounds.Dy())
333	return
334}
335
336func (p *gausssianBlurFilter) Draw(dst draw.Image, src image.Image, options *Options) {
337	if options == nil {
338		options = &defaultOptions
339	}
340
341	srcb := src.Bounds()
342	if srcb.Dx() <= 0 || srcb.Dy() <= 0 {
343		return
344	}
345
346	if p.sigma <= 0 {
347		copyimage(dst, src, options)
348		return
349	}
350
351	radius := int(math.Ceil(float64(p.sigma * 3)))
352	size := 2*radius + 1
353	center := radius
354	kernel := make([]float32, size)
355
356	kernel[center] = gaussianBlurKernel(0, p.sigma)
357	sum := kernel[center]
358
359	for i := 1; i <= radius; i++ {
360		f := gaussianBlurKernel(float32(i), p.sigma)
361		kernel[center-i] = f
362		kernel[center+i] = f
363		sum += 2 * f
364	}
365
366	for i := 0; i < len(kernel); i++ {
367		kernel[i] /= sum
368	}
369
370	tmp := createTempImage(srcb)
371	convolve1dh(tmp, src, kernel, options)
372	convolve1dv(dst, tmp, kernel, options)
373}
374
375// GaussianBlur creates a filter that applies a gaussian blur to an image.
376// The sigma parameter must be positive and indicates how much the image will be blurred.
377// Blur affected radius roughly equals 3 * sigma.
378//
379// Example:
380//
381//	g := gift.New(
382//		gift.GaussianBlur(1.5),
383//	)
384//	dst := image.NewRGBA(g.Bounds(src.Bounds()))
385//	g.Draw(dst, src)
386//
387func GaussianBlur(sigma float32) Filter {
388	return &gausssianBlurFilter{
389		sigma: sigma,
390	}
391}
392
393type unsharpMaskFilter struct {
394	sigma     float32
395	amount    float32
396	threshold float32
397}
398
399func (p *unsharpMaskFilter) Bounds(srcBounds image.Rectangle) (dstBounds image.Rectangle) {
400	dstBounds = image.Rect(0, 0, srcBounds.Dx(), srcBounds.Dy())
401	return
402}
403
404func unsharp(orig, blurred, amount, threshold float32) float32 {
405	dif := (orig - blurred) * amount
406	if absf32(dif) > absf32(threshold) {
407		return orig + dif
408	}
409	return orig
410}
411
412func (p *unsharpMaskFilter) Draw(dst draw.Image, src image.Image, options *Options) {
413	if options == nil {
414		options = &defaultOptions
415	}
416
417	srcb := src.Bounds()
418	dstb := dst.Bounds()
419
420	if srcb.Dx() <= 0 || srcb.Dy() <= 0 {
421		return
422	}
423
424	blurred := createTempImage(srcb)
425	blur := GaussianBlur(p.sigma)
426	blur.Draw(blurred, src, options)
427
428	pixGetterOrig := newPixelGetter(src)
429	pixGetterBlur := newPixelGetter(blurred)
430	pixelSetter := newPixelSetter(dst)
431
432	parallelize(options.Parallelization, srcb.Min.Y, srcb.Max.Y, func(pmin, pmax int) {
433		for y := pmin; y < pmax; y++ {
434			for x := srcb.Min.X; x < srcb.Max.X; x++ {
435				pxOrig := pixGetterOrig.getPixel(x, y)
436				pxBlur := pixGetterBlur.getPixel(x, y)
437
438				r := unsharp(pxOrig.r, pxBlur.r, p.amount, p.threshold)
439				g := unsharp(pxOrig.g, pxBlur.g, p.amount, p.threshold)
440				b := unsharp(pxOrig.b, pxBlur.b, p.amount, p.threshold)
441				a := unsharp(pxOrig.a, pxBlur.a, p.amount, p.threshold)
442
443				pixelSetter.setPixel(dstb.Min.X+x-srcb.Min.X, dstb.Min.Y+y-srcb.Min.Y, pixel{r, g, b, a})
444			}
445		}
446	})
447}
448
449// UnsharpMask creates a filter that sharpens an image.
450// The sigma parameter is used in a gaussian function and affects the radius of effect.
451// Sigma must be positive. Sharpen radius roughly equals 3 * sigma.
452// The amount parameter controls how much darker and how much lighter the edge borders become. Typically between 0.5 and 1.5.
453// The threshold parameter controls the minimum brightness change that will be sharpened. Typically between 0 and 0.05.
454//
455// Example:
456//
457//	g := gift.New(
458//		gift.UnsharpMask(1, 1, 0),
459//	)
460//	dst := image.NewRGBA(g.Bounds(src.Bounds()))
461//	g.Draw(dst, src)
462//
463func UnsharpMask(sigma, amount, threshold float32) Filter {
464	return &unsharpMaskFilter{
465		sigma:     sigma,
466		amount:    amount,
467		threshold: threshold,
468	}
469}
470
471type meanFilter struct {
472	ksize int
473	disk  bool
474}
475
476func (p *meanFilter) Bounds(srcBounds image.Rectangle) (dstBounds image.Rectangle) {
477	dstBounds = image.Rect(0, 0, srcBounds.Dx(), srcBounds.Dy())
478	return
479}
480
481func (p *meanFilter) Draw(dst draw.Image, src image.Image, options *Options) {
482	if options == nil {
483		options = &defaultOptions
484	}
485
486	srcb := src.Bounds()
487	if srcb.Dx() <= 0 || srcb.Dy() <= 0 {
488		return
489	}
490
491	ksize := p.ksize
492	if ksize%2 == 0 {
493		ksize--
494	}
495
496	if ksize <= 1 {
497		copyimage(dst, src, options)
498		return
499	}
500
501	if p.disk {
502		diskKernel := genDisk(p.ksize)
503		f := Convolution(diskKernel, true, true, false, 0)
504		f.Draw(dst, src, options)
505	} else {
506		kernel := make([]float32, ksize*ksize)
507		for i := range kernel {
508			kernel[i] = 1
509		}
510		f := Convolution(kernel, true, true, false, 0)
511		f.Draw(dst, src, options)
512	}
513}
514
515// Mean creates a local mean image filter.
516// Takes an average across a neighborhood for each pixel.
517// The ksize parameter is the kernel size. It must be an odd positive integer (for example: 3, 5, 7).
518// If the disk parameter is true, a disk-shaped neighborhood will be used instead of a square neighborhood.
519func Mean(ksize int, disk bool) Filter {
520	return &meanFilter{
521		ksize: ksize,
522		disk:  disk,
523	}
524}
525
526type hvConvolutionFilter struct {
527	hkernel, vkernel []float32
528}
529
530func (p *hvConvolutionFilter) Bounds(srcBounds image.Rectangle) (dstBounds image.Rectangle) {
531	dstBounds = image.Rect(0, 0, srcBounds.Dx(), srcBounds.Dy())
532	return
533}
534
535func (p *hvConvolutionFilter) Draw(dst draw.Image, src image.Image, options *Options) {
536	if options == nil {
537		options = &defaultOptions
538	}
539
540	srcb := src.Bounds()
541	dstb := dst.Bounds()
542
543	if srcb.Dx() <= 0 || srcb.Dy() <= 0 {
544		return
545	}
546
547	tmph := createTempImage(srcb)
548	Convolution(p.hkernel, false, false, true, 0).Draw(tmph, src, options)
549	pixGetterH := newPixelGetter(tmph)
550
551	tmpv := createTempImage(srcb)
552	Convolution(p.vkernel, false, false, true, 0).Draw(tmpv, src, options)
553	pixGetterV := newPixelGetter(tmpv)
554
555	pixSetter := newPixelSetter(dst)
556
557	parallelize(options.Parallelization, srcb.Min.Y, srcb.Max.Y, func(pmin, pmax int) {
558		for y := pmin; y < pmax; y++ {
559			for x := srcb.Min.X; x < srcb.Max.X; x++ {
560				pxh := pixGetterH.getPixel(x, y)
561				pxv := pixGetterV.getPixel(x, y)
562				r := sqrtf32(pxh.r*pxh.r + pxv.r*pxv.r)
563				g := sqrtf32(pxh.g*pxh.g + pxv.g*pxv.g)
564				b := sqrtf32(pxh.b*pxh.b + pxv.b*pxv.b)
565				pixSetter.setPixel(dstb.Min.X+x-srcb.Min.X, dstb.Min.Y+y-srcb.Min.Y, pixel{r, g, b, pxh.a})
566			}
567		}
568	})
569
570}
571
572// Sobel creates a filter that applies a sobel operator to an image.
573func Sobel() Filter {
574	return &hvConvolutionFilter{
575		hkernel: []float32{-1, 0, 1, -2, 0, 2, -1, 0, 1},
576		vkernel: []float32{-1, -2, -1, 0, 0, 0, 1, 2, 1},
577	}
578}