Quantized DCT coeficients, tweaking parameters

main.go

@@ -16,11 +16,18 @@ import (
 
 const (
 	BlockSize         int = 4
-	FrequencyThreshold     = 0.2
+	DistanceThreshold     = 0.4
-	OffsetThreshold        = 276
+	OffsetThreshold       = 72
-	ForgeryThreshold       = 305
+	ForgeryThreshold      = 220
 )
 
+var q4x4 = [][]float64{
+	{16.0, 10.0, 24.0, 51.0},
+	{14.0, 16.0, 40.0, 69.0},
+	{18.0, 37.0, 68.0, 103.0},
+	{49.0, 78.0, 103.0, 120.0},
+}
+
 // pixel struct contains the discrete cosine transformation R,G,B,Y values.
 type pixel struct {
 	r, g, b, y float64
@@ -40,14 +47,14 @@ type imageBlock struct {
 type vector struct {
 	xa, ya           int
 	xb, yb           int
-	offsetX, offsetY int
+	offsetX, offsetY float64
 }
 
 // feature struct contains the feature blocks x, y position and their respective values.
 type feature struct {
 	x    int
 	y    int
-	val float64
+	coef float64
 }
 
 var (
@@ -57,26 +64,34 @@ var (
 )
 
 func main() {
+	start := time.Now()
+
 	input, err := os.Open("parade_forged.jpg")
 	defer input.Close()
 
 	if err != nil {
 		fmt.Printf("Error reading the image file: %v", err)
 	}
-	img, _, err := image.Decode(input)
+	src, _, err := image.Decode(input)
 	if err != nil {
 		fmt.Printf("Error decoding the image: %v", err)
 	}
 
-	start := time.Now()
+	img := imgToNRGBA(src)
+	output := image.NewRGBA(img.Bounds())
+	draw.Draw(output, image.Rect(0, 0, img.Bounds().Dx(), img.Bounds().Dy()), img, image.ZP, draw.Src)
+
+	// Blur the image to eliminate the details.
+	blurImg := StackBlur(img, 1)
 
 	// Convert image to YUV color space
-	yuv := convertRGBImageToYUV(img)
+	yuv := convertRGBImageToYUV(blurImg)
 	newImg := image.NewRGBA(yuv.Bounds())
 	draw.Draw(newImg, image.Rect(0, 0, yuv.Bounds().Dx(), yuv.Bounds().Dy()), yuv, image.ZP, draw.Src)
 
 	dx, dy := yuv.Bounds().Max.X, yuv.Bounds().Max.Y
 	bdx, bdy := (dx - BlockSize + 1), (dy - BlockSize + 1)
+	n := math.Max(float64(dx), float64(dy))
 
 	var blocks []imageBlock
 	for i := 0; i < bdx; i++ {
@@ -84,25 +99,16 @@ func main() {
 			r := image.Rect(i, j, i+BlockSize, j+BlockSize)
 			block := newImg.SubImage(r).(*image.RGBA)
 			blocks = append(blocks, imageBlock{x: i, y: j, img: block})
-			draw.Draw(newImg, image.Rect(0, 0, yuv.Bounds().Max.X, yuv.Bounds().Max.Y), block, image.ZP, draw.Src)
+			//draw.Draw(newImg, image.Rect(0, 0, yuv.Bounds().Max.X, yuv.Bounds().Max.Y), block, image.ZP, draw.Src)
 		}
 	}
 
-	fmt.Printf("Len: %d", len(blocks))
+	fmt.Printf("Len: %d\n", len(blocks))
-
-	out, err := os.Create("output.png")
-	if err != nil {
-		fmt.Printf("Error creating output file: %v", err)
-	}
-
-	if err := png.Encode(out, newImg); err != nil {
-		fmt.Printf("Error encoding image file: %v", err)
-	}
 
+	for _, block := range blocks {
 		// Average RGB value.
 		var avr, avg, avb float64
 
-	for _, block := range blocks {
 		b := block.img.(*image.RGBA)
 		i0 := b.PixOffset(b.Bounds().Min.X, b.Bounds().Min.Y)
 		i1 := i0 + b.Bounds().Dx() * 4
@@ -135,36 +141,42 @@ func main() {
 				// normalization
 				alpha := func(a float64) float64 {
 					if a == 0 {
-						return math.Sqrt(1.0 / float64(dx))
+						return math.Sqrt(1.0 / float64(n))
 					} else {
-						return math.Sqrt(2.0 / float64(dy))
+						return math.Sqrt(2.0 / float64(n))
 					}
 				}
 
-				fi, fj := float64(u), float64(v)
+				cu, cv := float64(u), float64(v)
-				cr *= alpha(fi) * alpha(fj)
+				cr *= alpha(cu) * alpha(cv)
-				cg *= alpha(fi) * alpha(fj)
+				cg *= alpha(cu) * alpha(cv)
-				cb *= alpha(fi) * alpha(fj)
+				cb *= alpha(cu) * alpha(cv)
-				cy *= alpha(fi) * alpha(fj)
+				cy *= alpha(cu) * alpha(cv)
 
 				dctPixels[u][v] = pixel{cr, cg, cb, cy}
+
+				// Get the quantized DCT coefficients.
+				dctPixels[u][v].r = (dctPixels[u][v].r / q4x4[u][v])
+				dctPixels[u][v].g = (dctPixels[u][v].g / q4x4[u][v])
+				dctPixels[u][v].b = (dctPixels[u][v].b / q4x4[u][v])
+				dctPixels[u][v].y = (dctPixels[u][v].y / q4x4[u][v])
 			}
 		}
 		avr /= float64(BlockSize * BlockSize)
 		avg /= float64(BlockSize * BlockSize)
 		avb /= float64(BlockSize * BlockSize)
 
-		features = append(features, feature{x: block.x, y: block.y, val: dctPixels[0][0].y})
+		features = append(features, feature{x: block.x, y: block.y, coef: dctPixels[0][0].y})
-		features = append(features, feature{x: block.x, y: block.y, val: dctPixels[0][1].y})
+		features = append(features, feature{x: block.x, y: block.y, coef: dctPixels[0][1].y})
-		features = append(features, feature{x: block.x, y: block.y, val: dctPixels[1][0].y})
+		features = append(features, feature{x: block.x, y: block.y, coef: dctPixels[1][0].y})
-		features = append(features, feature{x: block.x, y: block.y, val: dctPixels[0][0].r})
+		features = append(features, feature{x: block.x, y: block.y, coef: dctPixels[0][0].r})
-		features = append(features, feature{x: block.x, y: block.y, val: dctPixels[0][0].g})
+		features = append(features, feature{x: block.x, y: block.y, coef: dctPixels[0][0].g})
-		features = append(features, feature{x: block.x, y: block.y, val: dctPixels[0][0].b})
+		features = append(features, feature{x: block.x, y: block.y, coef: dctPixels[0][0].b})
 
 		// Append average R,G,B values to the features vector(slice).
-		features = append(features, feature{x: block.x, y: block.y, val: avr})
+		features = append(features, feature{x: block.x, y: block.y, coef: avr})
-		features = append(features, feature{x: block.x, y: block.y, val: avb})
+		features = append(features, feature{x: block.x, y: block.y, coef: avb})
-		features = append(features, feature{x: block.x, y: block.y, val: avg})
+		features = append(features, feature{x: block.x, y: block.y, coef: avg})
 	}
 
 	// Lexicographically sort the feature vectors
@@ -180,7 +192,21 @@ func main() {
 	}
 
 	simBlocks := getSuspiciousBlocks(vectors)
-	_, result := filterOutNeighbors(simBlocks)
+	forgedBlocks, result := filterOutNeighbors(simBlocks)
+
+	for _, bl := range forgedBlocks {
+		background := color.RGBA{255, 0, 0, 255}
+		draw.Draw(output, image.Rect(bl.xa, bl.ya, bl.xa+BlockSize, bl.ya+BlockSize), &image.Uniform{background}, image.ZP, draw.Src)
+	}
+
+	out, err := os.Create("output.png")
+	if err != nil {
+		fmt.Printf("Error creating output file: %v", err)
+	}
+
+	if err := png.Encode(out, output); err != nil {
+		fmt.Printf("Error encoding image file: %v", err)
+	}
 
 	fmt.Println("\n", result)
 
@@ -217,18 +243,18 @@ func analyzeBlocks(blockA, blockB feature) *vector {
 		ya:      blockA.y,
 		xb:      blockB.x,
 		yb:      blockB.y,
-		offsetX: int(math.Abs(dx)),
+		offsetX: math.Abs(dx),
-		offsetY: int(math.Abs(dy)),
+		offsetY: math.Abs(dy),
 	}
 
-	if dist < FrequencyThreshold {
+	if dist < DistanceThreshold {
 		return res
 	}
 	return nil
 }
 
 type offset struct {
-	x, y int
+	x, y float64
 }
 
 type newVector []vector
@@ -261,6 +287,7 @@ func getSuspiciousBlocks(vect []vector) newVector {
 			})
 		}
 	}
+	fmt.Println("Blocks: ", len(suspiciousBlocks))
 	return suspiciousBlocks
 }
 
@@ -271,6 +298,7 @@ func filterOutNeighbors(vect []vector) (newVector, bool) {
 
 	for i := 1; i < len(vect); i++ {
 		blockA, blockB := vect[i-1], vect[i]
+
 		// Continue only if two regions are not neighbors.
 		if blockA.xa != blockB.xa && blockA.ya != blockB.ya {
 			// Calculate the euclidean distance between both regions.
@@ -280,9 +308,7 @@ func filterOutNeighbors(vect []vector) (newVector, bool) {
 
 			// Evaluate the euclidean distance distance between two regions
 			// and make sure the distance is greater than a predefined threshold.
-			// TODO verify threshold value
 			if dist > ForgeryThreshold {
-				fmt.Println(vect[i].xa, vect[i].ya)
 				forgedBlocks = append(forgedBlocks, vector{
 					vect[i].xa, vect[i].ya, vect[i].xb, vect[i].yb, vect[i].offsetX, vect[i].offsetY,
 				})
@@ -318,22 +344,16 @@ func idct(u, v, x, y, w float64) float64 {
 	return dct(u, v, x, y, w) * alpha(u) * alpha(v)
 }
 
-func RGBtoYUV(r, g, b uint32) (uint32, uint32, uint32) {
+// round rounds float number to it's nearest integer part.
-	y := 0.299*float64(r) + 0.587*float64(g) + 0.114*float64(b)
+func round(x float64) float64 {
-	u := (((float64(b) - float64(y)) * 0.493) + 111) / 222 * 255
+	t := math.Trunc(x)
-	v := (((float64(r) - float64(y)) * 0.877) + 155) / 312 * 255
+	if math.Abs(x-t) >= 0.5 {
-
+		return t + math.Copysign(1, x)
-	return uint32(y), uint32(u), uint32(v)
+	}
-}
+	return t
-
-func YUVtoRGB(y, u, v uint32) (uint32, uint32, uint32) {
-	r := float64(y) + (1.13983 * float64(v))
-	g := float64(y) - (0.39465 * float64(u)) - (0.58060 * float64(v))
-	b := float64(y) + (2.03211 * float64(u))
-
-	return uint32(r), uint32(g), uint32(b)
 }
 
+// clamp255 converts a float64 number to uint8.
 func clamp255(x float64) uint8 {
 	if x < 0 {
 		return 0
@@ -371,11 +391,85 @@ type featVec []feature
 func (a featVec) Len() int      { return len(a) }
 func (a featVec) Swap(i, j int) { a[i], a[j] = a[j], a[i] }
 func (a featVec) Less(i, j int) bool {
-	if a[i].val < a[j].val {
+	if a[i].coef < a[j].coef {
 		return true
 	}
-	if a[i].val > a[j].val {
+	if a[i].coef > a[j].coef {
 		return false
 	}
-	return a[i].val < a[j].val
+	return a[i].coef < a[j].coef
+}
+
+func RGBtoYUV(r, g, b uint32) (uint32, uint32, uint32) {
+	y := 0.299*float64(r) + 0.587*float64(g) + 0.114*float64(b)
+	u := (((float64(b) - float64(y)) * 0.493) + 111) / 222 * 255
+	v := (((float64(r) - float64(y)) * 0.877) + 155) / 312 * 255
+
+	return uint32(y), uint32(u), uint32(v)
+}
+
+func YUVtoRGB(y, u, v uint32) (uint32, uint32, uint32) {
+	r := float64(y) + (1.13983 * float64(v))
+	g := float64(y) - (0.39465 * float64(u)) - (0.58060 * float64(v))
+	b := float64(y) + (2.03211 * float64(u))
+
+	return uint32(r), uint32(g), uint32(b)
+}
+
+// Converts any image type to *image.NRGBA with min-point at (0, 0).
+func imgToNRGBA(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
 }

stackblur.go

@@ -0,0 +1,360 @@
+// Go implementation of StackBlur algorithm described here:
+// http://incubator.quasimondo.com/processing/fast_blur_deluxe.php
+
+package main
+
+import (
+	"image"
+)
+
+// blurStack is a linked list containing the color value and a pointer to the next struct.
+type blurStack struct {
+	r, g, b, a uint32
+	next       *blurStack
+}
+
+var mulTable = []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 shgTable = []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,
+}
+
+// NewBlurStack is a constructor function returning a new struct of type blurStack.
+func (bs *blurStack) NewBlurStack() *blurStack {
+	return &blurStack{bs.r, bs.g, bs.b, bs.a, bs.next}
+}
+
+// StackBlur applies a blur filter to the provided image.
+// The radius defines the bluring average.
+func StackBlur(img *image.NRGBA, radius uint32) *image.NRGBA {
+	var stackEnd, stackIn, stackOut *blurStack
+	var width, height = uint32(img.Bounds().Dx()), uint32(img.Bounds().Dy())
+	var (
+		div, widthMinus1, heightMinus1, radiusPlus1, sumFactor uint32
+		x, y, i, p, yp, yi, yw,
+		rSum, gSum, bSum, aSum,
+		rOutSum, gOutSum, bOutSum, aOutSum,
+		rInSum, gInSum, bInSum, aInSum,
+		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
+
+	mulSum := mulTable[radius]
+	shgSum := shgTable[radius]
+
+	for y = 0; y < height; y++ {
+		rInSum, gInSum, bInSum, aInSum, rSum, gSum, bSum, aSum = 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])
+
+		rOutSum = radiusPlus1 * pr
+		gOutSum = radiusPlus1 * pg
+		bOutSum = radiusPlus1 * pb
+		aOutSum = radiusPlus1 * pa
+
+		rSum += sumFactor * pr
+		gSum += sumFactor * pg
+		bSum += sumFactor * pb
+		aSum += 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
+
+			rSum += stack.r * (radiusPlus1 - i)
+			gSum += stack.g * (radiusPlus1 - i)
+			bSum += stack.b * (radiusPlus1 - i)
+			aSum += stack.a * (radiusPlus1 - i)
+
+			rInSum += pr
+			gInSum += pg
+			bInSum += pb
+			aInSum += pa
+
+			stack = stack.next
+		}
+		stackIn = stackStart
+		stackOut = stackEnd
+
+		for x = 0; x < width; x++ {
+			pa = (aSum * mulSum) >> shgSum
+			img.Pix[yi+3] = uint8(pa)
+
+			if pa != 0 {
+				img.Pix[yi] = uint8((rSum * mulSum) >> shgSum)
+				img.Pix[yi+1] = uint8((gSum * mulSum) >> shgSum)
+				img.Pix[yi+2] = uint8((bSum * mulSum) >> shgSum)
+			} else {
+				img.Pix[yi] = 0
+				img.Pix[yi+1] = 0
+				img.Pix[yi+2] = 0
+			}
+
+			rSum -= rOutSum
+			gSum -= gOutSum
+			bSum -= bOutSum
+			aSum -= aOutSum
+
+			rOutSum -= stackIn.r
+			gOutSum -= stackIn.g
+			bOutSum -= stackIn.b
+			aOutSum -= 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])
+
+			rInSum += stackIn.r
+			gInSum += stackIn.g
+			bInSum += stackIn.b
+			aInSum += stackIn.a
+
+			rSum += rInSum
+			gSum += gInSum
+			bSum += bInSum
+			aSum += aInSum
+
+			stackIn = stackIn.next
+
+			pr = stackOut.r
+			pg = stackOut.g
+			pb = stackOut.b
+			pa = stackOut.a
+
+			rOutSum += pr
+			gOutSum += pg
+			bOutSum += pb
+			aOutSum += pa
+
+			rInSum -= pr
+			gInSum -= pg
+			bInSum -= pb
+			aInSum -= pa
+
+			stackOut = stackOut.next
+
+			yi += 4
+		}
+		yw += width
+	}
+
+	for x = 0; x < width; x++ {
+		rInSum, gInSum, bInSum, aInSum, rSum, gSum, bSum, aSum = 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])
+
+		rOutSum = radiusPlus1 * pr
+		gOutSum = radiusPlus1 * pg
+		bOutSum = radiusPlus1 * pb
+		aOutSum = radiusPlus1 * pa
+
+		rSum += sumFactor * pr
+		gSum += sumFactor * pg
+		bSum += sumFactor * pb
+		aSum += 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
+
+			rSum += stack.r * (radiusPlus1 - i)
+			gSum += stack.g * (radiusPlus1 - i)
+			bSum += stack.b * (radiusPlus1 - i)
+			aSum += stack.a * (radiusPlus1 - i)
+
+			rInSum += pr
+			gInSum += pg
+			bInSum += pb
+			aInSum += 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 = (aSum * mulSum) >> shgSum
+			img.Pix[p+3] = uint8(pa)
+
+			if pa > 0 {
+				img.Pix[p] = uint8((rSum * mulSum) >> shgSum)
+				img.Pix[p+1] = uint8((gSum * mulSum) >> shgSum)
+				img.Pix[p+2] = uint8((bSum * mulSum) >> shgSum)
+			} else {
+				img.Pix[p] = 0
+				img.Pix[p+1] = 0
+				img.Pix[p+2] = 0
+			}
+
+			rSum -= rOutSum
+			gSum -= gOutSum
+			bSum -= bOutSum
+			aSum -= aOutSum
+
+			rOutSum -= stackIn.r
+			gOutSum -= stackIn.g
+			bOutSum -= stackIn.b
+			aOutSum -= 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])
+
+			rInSum += stackIn.r
+			gInSum += stackIn.g
+			bInSum += stackIn.b
+			aInSum += stackIn.a
+
+			rSum += rInSum
+			gSum += gInSum
+			bSum += bInSum
+			aSum += aInSum
+
+			stackIn = stackIn.next
+
+			pr = stackOut.r
+			pg = stackOut.g
+			pb = stackOut.b
+			pa = stackOut.a
+
+			rOutSum += pr
+			gOutSum += pg
+			bOutSum += pb
+			aOutSum += pa
+
+			rInSum -= pr
+			gInSum -= pg
+			bInSum -= pb
+			aInSum -= pa
+
+			stackOut = stackOut.next
+
+			yi += width
+		}
+	}
+	return img
+}