diff --git a/main.go b/main.go index ce38c3c..d81e441 100644 --- a/main.go +++ b/main.go @@ -15,12 +15,19 @@ import ( ) const ( - BlockSize int = 4 - FrequencyThreshold = 0.2 - OffsetThreshold = 276 - ForgeryThreshold = 305 + BlockSize int = 4 + DistanceThreshold = 0.4 + OffsetThreshold = 72 + 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 + x int + y int + 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,36 +99,27 @@ 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)) - - 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) - } - - // Average RGB value. - var avr, avg, avb float64 + fmt.Printf("Len: %d\n", len(blocks)) for _, block := range blocks { + // Average RGB value. + var avr, avg, avb float64 + b := block.img.(*image.RGBA) i0 := b.PixOffset(b.Bounds().Min.X, b.Bounds().Min.Y) - i1 := i0 + b.Bounds().Dx()*4 + i1 := i0 + b.Bounds().Dx() * 4 - dctPixels := make(dctPx, BlockSize*BlockSize) + dctPixels := make(dctPx, BlockSize * BlockSize) for u := 0; u < BlockSize; u++ { dctPixels[u] = make([]pixel, BlockSize) for v := 0; v < BlockSize; v++ { for i := i0; i < i1; i += 4 { // Get the YUV converted image pixels - yc, uc, vc, _ := b.Pix[i+0], b.Pix[i+2], b.Pix[i+2], b.Pix[i+3] + yc, uc, vc, _ := b.Pix[i + 0], b.Pix[i + 2], b.Pix[i + 2], b.Pix[i + 3] // Convert YUV to RGB and obtain the R value r, g, b := color.YCbCrToRGB(yc, uc, vc) @@ -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) - cr *= alpha(fi) * alpha(fj) - cg *= alpha(fi) * alpha(fj) - cb *= alpha(fi) * alpha(fj) - cy *= alpha(fi) * alpha(fj) + cu, cv := float64(u), float64(v) + cr *= alpha(cu) * alpha(cv) + cg *= alpha(cu) * alpha(cv) + cb *= alpha(cu) * alpha(cv) + 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, val: 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, val: 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, val: dctPixels[0][0].b}) + 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, coef: dctPixels[0][1].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, coef: dctPixels[0][0].r}) + 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, 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, val: avb}) - features = append(features, feature{x: block.x, y: block.y, val: avg}) + features = append(features, feature{x: block.x, y: block.y, coef: avr}) + features = append(features, feature{x: block.x, y: block.y, coef: avb}) + 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)), - offsetY: int(math.Abs(dy)), + offsetX: math.Abs(dx), + 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) { - 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) +// round rounds float number to it's nearest integer part. +func round(x float64) float64 { + t := math.Trunc(x) + if math.Abs(x-t) >= 0.5 { + return t + math.Copysign(1, x) + } + return t } +// 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 } diff --git a/stackblur.go b/stackblur.go new file mode 100644 index 0000000..3923f6f --- /dev/null +++ b/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 +}