Initial commit

carver.go

@@ -0,0 +1,136 @@
+package caire
+
+import (
+	"image"
+	"fmt"
+	"math"
+)
+
+type DPTable struct {
+	width int
+	height int
+	table []float64
+}
+
+type Seam []Point
+
+type Point struct {
+	X int
+	Y int
+}
+
+// Get energy pixel value
+func (dpt *DPTable) get(x, y int) float64 {
+	px := x + y * dpt.width
+	return dpt.table[px]
+}
+
+// Set energy pixel value
+func (dpt *DPTable) set(x, y int, px float64) {
+	idx := x + y * dpt.width
+	dpt.table[idx] = px
+}
+
+// Compute the minimum energy level based on the following logic:
+// 	- traverse the image from the second row to the last row
+// 	  and compute the cumulative minimum energy M for all possible
+//	  connected seams for each entry (i, j).
+//
+//	- the minimum energy level is calculated by summing up the current pixel value
+// 	  with the minimum pixel value of the neighboring pixels from the previous row.
+func (dpt *DPTable) ComputeSeams(img *image.NRGBA) []float64 {
+	width, height := img.Bounds().Max.X, img.Bounds().Max.Y
+
+	for x := 0; x < width; x++ {
+		for y := 0; y < height; y++ {
+			r, _, _, a := img.At(x, y).RGBA()
+			dpt.set(x, y, float64(r) / float64(a))
+		}
+	}
+
+	computeEnergyLevel := func(x, y int) {
+		var left, middle, right float64
+		left, right = dpt.get(x, y-1), dpt.get(x, y-1)
+
+		// Do not compute edge cases: pixels are far left.
+		if x > 0 {
+			left = dpt.get(x-1, y-1)
+		}
+		middle = dpt.get(x, y-1)
+		// Do not compute edge cases: pixels are far right.
+		if x < width-1 {
+			right = dpt.get(x+1, y-1)
+		}
+		// Obtain the minimum pixel value
+		min := math.Min(math.Min(left, middle), right)
+		dpt.set(x, y, dpt.get(x, y) + min)
+	}
+
+	for x := 0; x < width; x++ {
+		for y := 1; y < height; y++ {
+			computeEnergyLevel(x, y)
+		}
+	}
+
+	//fmt.Println(dpt.table)
+	return dpt.table
+}
+
+// Find the lowest vertical energy seam.
+func (dpt *DPTable) FindLowestEnergySeams() []Point {
+	// Find the lowest cost seam from the energy matrix starting from the last row.
+	var min float64 = math.MaxFloat64
+	var px int
+	seams := make([]Point, 0)
+
+	// Find the lowest seam from the bottom row
+	for x := 0; x < dpt.width; x++ {
+		seam := dpt.get(x, dpt.height-1)
+		if seam < min {
+			min = seam
+			px = x
+		}
+	}
+	seams = append(seams, Point{X: px, Y: dpt.height-1})
+
+	// Walk up in the matrix table,
+	// check the immediate three top pixel seam level and
+	// add add the one which has the lowest cumulative energy.
+	for y := dpt.height-2; y >= 0; y-- {
+		left, center, right := math.MaxFloat64, math.MaxFloat64, math.MaxFloat64
+		// Leftmost seam, no child to the left
+		if px == 0 {
+			right = dpt.get(px+1, y)
+			center = dpt.get(px, y)
+			if right < center {
+				px += 1
+			}
+		// Rightmost seam, no child to the right
+		} else if px == dpt.width-1 {
+			left = dpt.get(px-1, y)
+			center = dpt.get(px, y)
+			if left < center {
+				px -= 1
+			}
+		} else {
+			left = dpt.get(px-1, y)
+			center = dpt.get(px, y)
+			right = dpt.get(px+1, y)
+			min := math.Min(math.Min(left, center), right)
+
+			if min == left {
+				px -= 1
+			} else if min == right {
+				px += 1
+			}
+		}
+		seams = append(seams, Point{X: px, Y: y})
+	}
+	fmt.Println(seams)
+	return seams
+}
+
+// Remove image pixels based on energy seams level
+func (dpt *DPTable) RemoveSeams(seams Seam) {
+
+}

cmd/caire/main.go

@@ -0,0 +1,145 @@
+package main
+
+import (
+	"flag"
+	"fmt"
+	"io/ioutil"
+	"log"
+	"os"
+	"path"
+	"path/filepath"
+	"strings"
+	"time"
+
+	"github.com/esimov/caire"
+)
+
+var (
+	// Flags
+	source         = flag.String("in", "", "Source")
+	destination    = flag.String("out", "", "Destination")
+	blurRadius     = flag.Int("blur", 1, "Blur radius")
+	sobelThreshold = flag.Int("sobel", 1, "Sobel filter threshold")
+)
+
+func main() {
+	flag.Parse()
+
+	if len(*source) == 0 || len(*destination) == 0 {
+		log.Fatal("Usage: caire -in input.jpg -out out.jpg")
+	}
+
+	fs, err := os.Stat(*source)
+	if err != nil {
+		log.Fatalf("Unable to open source: %v", err)
+	}
+
+	toProcess := make(map[string]string)
+
+	p := &caire.Processor{
+		BlurRadius:     *blurRadius,
+		SobelThreshold: *sobelThreshold,
+	}
+
+	switch mode := fs.Mode(); {
+	case mode.IsDir():
+		// Supported image files.
+		extensions := []string{".jpg", ".png"}
+
+		// Read source directory.
+		files, err := ioutil.ReadDir(*source)
+		if err != nil {
+			log.Fatalf("Unable to read dir: %v", err)
+		}
+
+		// Read destination file or directory.
+		dst, err := os.Stat(*destination)
+		if err != nil {
+			log.Fatalf("Unable to get dir stats: %v", err)
+		}
+
+		// Check if the image destination is a directory or a file.
+		if dst.Mode().IsRegular() {
+			log.Fatal("Please specify a directory as destination!")
+			os.Exit(2)
+		}
+		output, err := filepath.Abs(filepath.Base(*destination))
+		if err != nil {
+			log.Fatalf("Unable to get absolute path: %v", err)
+		}
+
+		// Range over all the image files and save them into a slice.
+		images := []string{}
+		for _, f := range files {
+			ext := filepath.Ext(f.Name())
+			for _, iex := range extensions {
+				if ext == iex {
+					images = append(images, f.Name())
+				}
+			}
+		}
+
+		// Process images from directory.
+		for _, img := range images {
+			// Get the file base name.
+			name := strings.TrimSuffix(img, filepath.Ext(img))
+			dir := strings.TrimRight(*source, "/")
+			out := output + "/" + name + ".png"
+			in := dir + "/" + img
+
+			toProcess[in] = out
+		}
+
+	case mode.IsRegular():
+		toProcess[*source] = *destination
+	}
+
+	for in, out := range toProcess {
+		file, err := os.Open(in)
+		if err != nil {
+			log.Fatalf("Unable to open source file: %v", err)
+		}
+		defer file.Close()
+
+		s := new(spinner)
+		s.start("Process start...")
+		start := time.Now()
+		_, processErr := p.Process(file, out)
+		s.stop()
+
+		if err == nil {
+			fmt.Printf("\nRescaled in: \x1b[92m%.2fs\n", time.Since(start).Seconds())
+			fmt.Printf("Saved as: %s \x1b[92m✓\n\n", path.Base(out))
+		} else {
+			fmt.Printf("\nError rescaling image: %s: %s", file.Name(), processErr.Error())
+		}
+	}
+}
+
+type spinner struct {
+	stopChan chan struct{}
+}
+
+// Start process
+func (s *spinner) start(message string) {
+	s.stopChan = make(chan struct{}, 1)
+
+	go func() {
+		for {
+			for _, r := range `-\|/` {
+				select {
+				case <-s.stopChan:
+					return
+				default:
+					fmt.Printf("\r%s%s %c%s", message, "\x1b[92m", r, "\x1b[39m")
+					time.Sleep(time.Millisecond * 100)
+				}
+			}
+		}
+	}()
+}
+
+// End process
+func (s *spinner) stop() {
+	s.stopChan <- struct{}{}
+}

grayscale.go

@@ -0,0 +1,21 @@
+package caire
+
+import (
+	"image"
+	"image/color"
+)
+
+// Convert image to grayscale.
+func Grayscale(src *image.NRGBA) *image.NRGBA {
+	dx, dy := src.Bounds().Max.X, src.Bounds().Max.Y
+	dst := image.NewNRGBA(src.Bounds())
+	for x := 0; x < dx; x++ {
+		for y := 0; y < dy; y++ {
+			r, g, b, _ := src.At(x, y).RGBA()
+			lum := float32(r)*0.299 + float32(g)*0.587 + float32(b)*0.114
+			pixel := color.Gray{uint8(lum / 256)}
+			dst.Set(x, y, pixel)
+		}
+	}
+	return dst
+}

process.go

@@ -0,0 +1,111 @@
+package caire
+
+import (
+	"image"
+	"image/color"
+	"image/png"
+	"io"
+	"os"
+	"github.com/fogleman/gg"
+)
+
+// Processor options
+type Processor struct {
+	BlurRadius      int
+	SobelThreshold  int
+}
+
+// Process : Convert image
+func (p *Processor) Process(file io.Reader, output string) (*os.File, error) {
+	src, _, err := image.Decode(file)
+	if err != nil {
+		return nil, err
+	}
+
+	width, height := src.Bounds().Dx(), src.Bounds().Dy()
+	ctx := gg.NewContext(width, height)
+	ctx.DrawRectangle(0, 0, float64(width), float64(height))
+	ctx.SetRGBA(1, 1, 1, 1)
+	ctx.Fill()
+
+	img := toNRGBA(src)
+
+	blur := Stackblur(img, uint32(width), uint32(height), uint32(p.BlurRadius))
+	gray := Grayscale(blur)
+	sobel := SobelFilter(gray, float64(p.SobelThreshold))
+
+	dpTable := &DPTable{width, height, make([]float64, width*height)}
+	dpTable.ComputeSeams(sobel)
+	dpTable.FindLowestEnergySeams()
+
+	fq, err := os.Create(output)
+	if err != nil {
+		return nil, err
+	}
+	defer fq.Close()
+
+	if err = png.Encode(fq, sobel); err != nil {
+		return nil, err
+	}
+
+	return fq, err
+}
+
+// toNRGBA converts any image type to *image.NRGBA with min-point at (0, 0).
+func toNRGBA(img image.Image) *image.NRGBA {
+	srcBounds := img.Bounds()
+	if srcBounds.Min.X == 0 && srcBounds.Min.Y == 0 {
+		if src0, ok := img.(*image.NRGBA); ok {
+			return src0
+		}
+	}
+	srcMinX := srcBounds.Min.X
+	srcMinY := srcBounds.Min.Y
+
+	dstBounds := srcBounds.Sub(srcBounds.Min)
+	dstW := dstBounds.Dx()
+	dstH := dstBounds.Dy()
+	dst := image.NewNRGBA(dstBounds)
+
+	switch src := img.(type) {
+	case *image.NRGBA:
+		rowSize := srcBounds.Dx() * 4
+		for dstY := 0; dstY < dstH; dstY++ {
+			di := dst.PixOffset(0, dstY)
+			si := src.PixOffset(srcMinX, srcMinY+dstY)
+			for dstX := 0; dstX < dstW; dstX++ {
+				copy(dst.Pix[di:di+rowSize], src.Pix[si:si+rowSize])
+			}
+		}
+	case *image.YCbCr:
+		for dstY := 0; dstY < dstH; dstY++ {
+			di := dst.PixOffset(0, dstY)
+			for dstX := 0; dstX < dstW; dstX++ {
+				srcX := srcMinX + dstX
+				srcY := srcMinY + dstY
+				siy := src.YOffset(srcX, srcY)
+				sic := src.COffset(srcX, srcY)
+				r, g, b := color.YCbCrToRGB(src.Y[siy], src.Cb[sic], src.Cr[sic])
+				dst.Pix[di+0] = r
+				dst.Pix[di+1] = g
+				dst.Pix[di+2] = b
+				dst.Pix[di+3] = 0xff
+				di += 4
+			}
+		}
+	default:
+		for dstY := 0; dstY < dstH; dstY++ {
+			di := dst.PixOffset(0, dstY)
+			for dstX := 0; dstX < dstW; dstX++ {
+				c := color.NRGBAModel.Convert(img.At(srcMinX+dstX, srcMinY+dstY)).(color.NRGBA)
+				dst.Pix[di+0] = c.R
+				dst.Pix[di+1] = c.G
+				dst.Pix[di+2] = c.B
+				dst.Pix[di+3] = c.A
+				di += 4
+			}
+		}
+	}
+
+	return dst
+}

sobel.go

@@ -0,0 +1,105 @@
+package caire
+
+import (
+	"image"
+	"math"
+)
+
+type kernel [][]int32
+
+var (
+	kernelX kernel = kernel{
+		{-1, 0, 1},
+		{-2, 0, 2},
+		{-1, 0, 1},
+	}
+
+	kernelY kernel = kernel{
+		{-1, -2, -1},
+		{0, 0, 0},
+		{1, 2, 1},
+	}
+)
+
+// Detect image edges.
+// See https://en.wikipedia.org/wiki/Sobel_operator
+func SobelFilter(img *image.NRGBA, threshold float64) *image.NRGBA {
+	var sumX, sumY int32
+	dx, dy := img.Bounds().Max.X, img.Bounds().Max.Y
+	dst := image.NewNRGBA(img.Bounds())
+
+	// Get 3x3 window of pixels because image data given is just a 1D array of pixels
+	maxPixelOffset := dx*dy + len(kernelX) - 1
+
+	data := getImageData(img)
+	length := len(data) - maxPixelOffset
+	magnitudes := make([]int32, length)
+
+	for i := 0; i < length; i++ {
+		// Sum each pixel with the kernel value
+		sumX, sumY = 0, 0
+		for x := 0; x < len(kernelX); x++ {
+			for y := 0; y < len(kernelY); y++ {
+				px := data[i+(dx*y)+x]
+				if len(px) > 0 {
+					r := px[0]
+					// We are using px[0] (i.e. R value) because the image is grayscale anyway
+					sumX += int32(r) * kernelX[y][x]
+					sumY += int32(r) * kernelY[y][x]
+				}
+			}
+		}
+		magnitude := math.Sqrt(float64(sumX*sumX) + float64(sumY*sumY))
+		// Check for pixel color boundaries
+		if magnitude < 0 {
+			magnitude = 0
+		} else if magnitude > 255 {
+			magnitude = 255
+		}
+
+		// Set magnitude to 0 if doesn't exceed threshold, else set to magnitude
+		if magnitude > threshold {
+			magnitudes[i] = int32(magnitude)
+		} else {
+			magnitudes[i] = 0
+		}
+	}
+
+	dataLength := dx * dy * 4
+	edges := make([]int32, dataLength)
+
+	// Apply the kernel values.
+	for i := 0; i < dataLength; i++ {
+		if i%4 != 0 {
+			m := magnitudes[i/4]
+			if m != 0 {
+				edges[i-1] = m
+			}
+		}
+	}
+
+	// Generate the new image with the sobel filter applied.
+	for idx := 0; idx < len(edges); idx += 4 {
+		dst.Pix[idx] = uint8(edges[idx])
+		dst.Pix[idx+1] = uint8(edges[idx+1])
+		dst.Pix[idx+2] = uint8(edges[idx+2])
+		dst.Pix[idx+3] = 255
+	}
+	return dst
+}
+
+// Group pixels into 2D array, each one containing the pixel RGB value.
+func getImageData(img *image.NRGBA) [][]uint8 {
+	dx, dy := img.Bounds().Max.X, img.Bounds().Max.Y
+	pixels := make([][]uint8, dx*dy*4)
+
+	for i := 0; i < len(pixels); i += 4 {
+		pixels[i/4] = []uint8{
+			img.Pix[i],
+			img.Pix[i+1],
+			img.Pix[i+2],
+			img.Pix[i+3],
+		}
+	}
+	return pixels
+}

stackblur.go

@@ -0,0 +1,357 @@
+// Go implementation of StackBlur algorithm described here:
+// http://incubator.quasimondo.com/processing/fast_blur_deluxe.php
+
+package caire
+
+import (
+	"image"
+)
+
+type blurstack struct {
+	r, g, b, a uint32
+	next       *blurstack
+}
+
+var mul_table []uint32 = []uint32{
+	512, 512, 456, 512, 328, 456, 335, 512, 405, 328, 271, 456, 388, 335, 292, 512,
+	454, 405, 364, 328, 298, 271, 496, 456, 420, 388, 360, 335, 312, 292, 273, 512,
+	482, 454, 428, 405, 383, 364, 345, 328, 312, 298, 284, 271, 259, 496, 475, 456,
+	437, 420, 404, 388, 374, 360, 347, 335, 323, 312, 302, 292, 282, 273, 265, 512,
+	497, 482, 468, 454, 441, 428, 417, 405, 394, 383, 373, 364, 354, 345, 337, 328,
+	320, 312, 305, 298, 291, 284, 278, 271, 265, 259, 507, 496, 485, 475, 465, 456,
+	446, 437, 428, 420, 412, 404, 396, 388, 381, 374, 367, 360, 354, 347, 341, 335,
+	329, 323, 318, 312, 307, 302, 297, 292, 287, 282, 278, 273, 269, 265, 261, 512,
+	505, 497, 489, 482, 475, 468, 461, 454, 447, 441, 435, 428, 422, 417, 411, 405,
+	399, 394, 389, 383, 378, 373, 368, 364, 359, 354, 350, 345, 341, 337, 332, 328,
+	324, 320, 316, 312, 309, 305, 301, 298, 294, 291, 287, 284, 281, 278, 274, 271,
+	268, 265, 262, 259, 257, 507, 501, 496, 491, 485, 480, 475, 470, 465, 460, 456,
+	451, 446, 442, 437, 433, 428, 424, 420, 416, 412, 408, 404, 400, 396, 392, 388,
+	385, 381, 377, 374, 370, 367, 363, 360, 357, 354, 350, 347, 344, 341, 338, 335,
+	332, 329, 326, 323, 320, 318, 315, 312, 310, 307, 304, 302, 299, 297, 294, 292,
+	289, 287, 285, 282, 280, 278, 275, 273, 271, 269, 267, 265, 263, 261, 259,
+}
+
+var shg_table []uint32 = []uint32{
+	9, 11, 12, 13, 13, 14, 14, 15, 15, 15, 15, 16, 16, 16, 16, 17,
+	17, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 18, 18, 18, 18, 19,
+	19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 20, 20, 20,
+	20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 21,
+	21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21,
+	21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 22, 22, 22, 22, 22, 22,
+	22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22,
+	22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 23,
+	23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23,
+	23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23,
+	23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23,
+	23, 23, 23, 23, 23, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24,
+	24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24,
+	24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24,
+	24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24,
+	24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24,
+}
+
+func (bs *blurstack) NewBlurStack() *blurstack {
+	return &blurstack{bs.r, bs.g, bs.b, bs.a, bs.next}
+}
+
+func Stackblur(img *image.NRGBA, width, height, radius uint32) *image.NRGBA {
+	var stackEnd, stackIn, stackOut *blurstack
+	var (
+		div, widthMinus1, heightMinus1, radiusPlus1, sumFactor uint32
+		x, y, i, p, yp, yi, yw,
+		r_sum, g_sum, b_sum, a_sum,
+		r_out_sum, g_out_sum, b_out_sum, a_out_sum,
+		r_in_sum, g_in_sum, b_in_sum, a_in_sum,
+		pr, pg, pb, pa uint32
+	)
+
+	div = radius + radius + 1
+	widthMinus1 = width - 1
+	heightMinus1 = height - 1
+	radiusPlus1 = radius + 1
+	sumFactor = radiusPlus1 * (radiusPlus1 + 1) / 2
+
+	bs := blurstack{}
+	stackStart := bs.NewBlurStack()
+	stack := stackStart
+
+	for i = 1; i < div; i++ {
+		stack.next = bs.NewBlurStack()
+		stack = stack.next
+		if i == radiusPlus1 {
+			stackEnd = stack
+		}
+	}
+	stack.next = stackStart
+
+	mul_sum := mul_table[radius]
+	shg_sum := shg_table[radius]
+
+	for y = 0; y < height; y++ {
+		r_in_sum, g_in_sum, b_in_sum, a_in_sum, r_sum, g_sum, b_sum, a_sum = 0, 0, 0, 0, 0, 0, 0, 0
+
+		pr = uint32(img.Pix[yi])
+		pg = uint32(img.Pix[yi+1])
+		pb = uint32(img.Pix[yi+2])
+		pa = uint32(img.Pix[yi+3])
+
+		r_out_sum = radiusPlus1 * pr
+		g_out_sum = radiusPlus1 * pg
+		b_out_sum = radiusPlus1 * pb
+		a_out_sum = radiusPlus1 * pa
+
+		r_sum += sumFactor * pr
+		g_sum += sumFactor * pg
+		b_sum += sumFactor * pb
+		a_sum += sumFactor * pa
+
+		stack = stackStart
+
+		for i = 0; i < radiusPlus1; i++ {
+			stack.r = pr
+			stack.g = pg
+			stack.b = pb
+			stack.a = pa
+			stack = stack.next
+		}
+
+		for i = 1; i < radiusPlus1; i++ {
+			var diff uint32
+			if widthMinus1 < i {
+				diff = widthMinus1
+			} else {
+				diff = i
+			}
+			p = yi + (diff << 2)
+			pr = uint32(img.Pix[p])
+			pg = uint32(img.Pix[p+1])
+			pb = uint32(img.Pix[p+2])
+			pa = uint32(img.Pix[p+3])
+
+			stack.r = pr
+			stack.g = pg
+			stack.b = pb
+			stack.a = pa
+
+			r_sum += stack.r * (radiusPlus1 - i)
+			g_sum += stack.g * (radiusPlus1 - i)
+			b_sum += stack.b * (radiusPlus1 - i)
+			a_sum += stack.a * (radiusPlus1 - i)
+
+			r_in_sum += pr
+			g_in_sum += pg
+			b_in_sum += pb
+			a_in_sum += pa
+
+			stack = stack.next
+		}
+		stackIn = stackStart
+		stackOut = stackEnd
+
+		for x = 0; x < width; x++ {
+			pa = (a_sum * mul_sum) >> shg_sum
+			img.Pix[yi+3] = uint8(pa)
+
+			if pa != 0 {
+				pa = 255 / pa
+				img.Pix[yi] = uint8((r_sum * mul_sum) >> shg_sum)
+				img.Pix[yi+1] = uint8((g_sum * mul_sum) >> shg_sum)
+				img.Pix[yi+2] = uint8((b_sum * mul_sum) >> shg_sum)
+			} else {
+				img.Pix[yi] = 0
+				img.Pix[yi+1] = 0
+				img.Pix[yi+2] = 0
+			}
+
+			r_sum -= r_out_sum
+			g_sum -= g_out_sum
+			b_sum -= b_out_sum
+			a_sum -= a_out_sum
+
+			r_out_sum -= stackIn.r
+			g_out_sum -= stackIn.g
+			b_out_sum -= stackIn.b
+			a_out_sum -= stackIn.a
+
+			p = x + radius + 1
+
+			if p > widthMinus1 {
+				p = widthMinus1
+			}
+			p = (yw + p) << 2
+
+			stackIn.r = uint32(img.Pix[p])
+			stackIn.g = uint32(img.Pix[p+1])
+			stackIn.b = uint32(img.Pix[p+2])
+			stackIn.a = uint32(img.Pix[p+3])
+
+			r_in_sum += stackIn.r
+			g_in_sum += stackIn.g
+			b_in_sum += stackIn.b
+			a_in_sum += stackIn.a
+
+			r_sum += r_in_sum
+			g_sum += g_in_sum
+			b_sum += b_in_sum
+			a_sum += a_in_sum
+
+			stackIn = stackIn.next
+
+			pr = stackOut.r
+			pg = stackOut.g
+			pb = stackOut.b
+			pa = stackOut.a
+
+			r_out_sum += pr
+			g_out_sum += pg
+			b_out_sum += pb
+			a_out_sum += pa
+
+			r_in_sum -= pr
+			g_in_sum -= pg
+			b_in_sum -= pb
+			a_in_sum -= pa
+
+			stackOut = stackOut.next
+
+			yi += 4
+		}
+		yw += width
+	}
+
+	for x = 0; x < width; x++ {
+		r_in_sum, g_in_sum, b_in_sum, a_in_sum, r_sum, g_sum, b_sum, a_sum = 0, 0, 0, 0, 0, 0, 0, 0
+
+		yi = x << 2
+		pr = uint32(img.Pix[yi])
+		pg = uint32(img.Pix[yi+1])
+		pb = uint32(img.Pix[yi+2])
+		pa = uint32(img.Pix[yi+3])
+
+		r_out_sum = radiusPlus1 * pr
+		g_out_sum = radiusPlus1 * pg
+		b_out_sum = radiusPlus1 * pb
+		a_out_sum = radiusPlus1 * pa
+
+		r_sum += sumFactor * pr
+		g_sum += sumFactor * pg
+		b_sum += sumFactor * pb
+		a_sum += sumFactor * pa
+
+		stack = stackStart
+
+		for i = 0; i < radiusPlus1; i++ {
+			stack.r = pr
+			stack.g = pg
+			stack.b = pb
+			stack.a = pa
+			stack = stack.next
+		}
+
+		yp = width
+
+		for i = 1; i <= radius; i++ {
+			yi = (yp + x) << 2
+			pr = uint32(img.Pix[yi])
+			pg = uint32(img.Pix[yi+1])
+			pb = uint32(img.Pix[yi+2])
+			pa = uint32(img.Pix[yi+3])
+
+			stack.r = pr
+			stack.g = pg
+			stack.b = pb
+			stack.a = pa
+
+			r_sum += stack.r * (radiusPlus1 - i)
+			g_sum += stack.g * (radiusPlus1 - i)
+			b_sum += stack.b * (radiusPlus1 - i)
+			a_sum += stack.a * (radiusPlus1 - i)
+
+			r_in_sum += pr
+			g_in_sum += pg
+			b_in_sum += pb
+			a_in_sum += pa
+
+			stack = stack.next
+
+			if i < heightMinus1 {
+				yp += width
+			}
+		}
+
+		yi = x
+		stackIn = stackStart
+		stackOut = stackEnd
+
+		for y = 0; y < height; y++ {
+			p = yi << 2
+			pa = (a_sum * mul_sum) >> shg_sum
+			img.Pix[p+3] = uint8(pa)
+
+			if pa > 0 {
+				pa = 255 / pa
+				img.Pix[p] = uint8((r_sum * mul_sum) >> shg_sum)
+				img.Pix[p+1] = uint8((g_sum * mul_sum) >> shg_sum)
+				img.Pix[p+2] = uint8((b_sum * mul_sum) >> shg_sum)
+			} else {
+				img.Pix[p] = 0
+				img.Pix[p+1] = 0
+				img.Pix[p+2] = 0
+			}
+
+			r_sum -= r_out_sum
+			g_sum -= g_out_sum
+			b_sum -= b_out_sum
+			a_sum -= a_out_sum
+
+			r_out_sum -= stackIn.r
+			g_out_sum -= stackIn.g
+			b_out_sum -= stackIn.b
+			a_out_sum -= stackIn.a
+
+			p = y + radiusPlus1
+
+			if p > heightMinus1 {
+				p = heightMinus1
+			}
+			p = (x + (p * width)) << 2
+
+			stackIn.r = uint32(img.Pix[p])
+			stackIn.g = uint32(img.Pix[p+1])
+			stackIn.b = uint32(img.Pix[p+2])
+			stackIn.a = uint32(img.Pix[p+3])
+
+			r_in_sum += stackIn.r
+			g_in_sum += stackIn.g
+			b_in_sum += stackIn.b
+			a_in_sum += stackIn.a
+
+			r_sum += r_in_sum
+			g_sum += g_in_sum
+			b_sum += b_in_sum
+			a_sum += a_in_sum
+
+			stackIn = stackIn.next
+
+			pr = stackOut.r
+			pg = stackOut.g
+			pb = stackOut.b
+			pa = stackOut.a
+
+			r_out_sum += pr
+			g_out_sum += pg
+			b_out_sum += pb
+			a_out_sum += pa
+
+			r_in_sum -= pr
+			g_in_sum -= pg
+			b_in_sum -= pb
+			a_in_sum -= pa
+
+			stackOut = stackOut.next
+
+			yi += width
+		}
+	}
+	return img
+}