forked from: Fitting gaussian to image histogram (3)
Fitting gaussian to image histogram.
How do we limit fit to dominant gaussian?
@author makc
/**
* Copyright makc3d ( http://wonderfl.net/user/makc3d )
* MIT License ( http://www.opensource.org/licenses/mit-license.php )
* Downloaded from: http://wonderfl.net/c/mtxi
*/
package {
import flash.display.BitmapData;
import flash.display.Sprite;
import flash.events.Event;
import flash.events.MouseEvent;
import flash.filters.ColorMatrixFilter;
import flash.geom.Rectangle;
import flash.media.Camera;
import flash.media.Video;
/**
* Fitting gaussian to image histogram.
* How do we limit fit to dominant gaussian?
* @author makc
*/
[SWF(width=465,height=465,backgroundColor='#003F00')]
public class BellCurveFitting extends Sprite {
public var video:Video;
public var frame:BitmapData;
public var grayscale:ColorMatrixFilter;
public var rectangle:Rectangle;
public var rectSprite:Sprite;
public function BellCurveFitting () {
video = new Video; addChild (video);
video.attachCamera (Camera.getCamera ());
frame = new BitmapData (320, 240, false, 0);
grayscale = new ColorMatrixFilter (
[ 0.299, 0.587, 0.114, 0, 0,
0.299, 0.587, 0.114, 0, 0,
0.299, 0.587, 0.114, 0, 0,
0, 0, 0, 1, 0 ]
);
rectangle = frame.rect;
addChild (rectSprite = new Sprite);
stage.addEventListener (Event.ENTER_FRAME, loop);
stage.addEventListener (MouseEvent.MOUSE_DOWN, rectStart);
}
public function loop (e:Event):void {
frame.draw (video);
frame.applyFilter (frame, frame.rect, frame.rect.topLeft, grayscale);
var area:Number = Math.abs (rectangle.width * rectangle.height);
var h:Vector.<Number> = frame.histogram (getPositiveRect (rectangle)) [0];
// find max
var h_max:Number = 0;
for (var i:int = 0; i < 256; i++) {
if (h [i] > h_max) h_max = h [i];
}
// graph h[i]
graphics.clear ();
for (i = 0; i < 256; i++) {
var v:int = 200 * h [i] / h_max;
graphics.beginFill (0x10101 * i);
graphics.drawRect (i, 465 - v, 1, v);
}
graphics.endFill ();
// fit and graph
for (var j:int = 0; j < 3; j++) {
var g:Gaussian = fitAndSubtract (h);
graphics.lineStyle (1, [0xFF0000, 0xFF00, 0xFF00FF][j]);
for (i = 0; i < 256; i++) {
v = 200 * g.f (i) / h_max;
graphics [(i < 1) ? "moveTo" : "lineTo"] (i, 465 - v);
}
}
// rect
rectSprite.graphics.clear ();
rectSprite.graphics.lineStyle (1, 0xFF0000);
rectSprite.graphics.drawRect (rectangle.x, rectangle.y, rectangle.width, rectangle.height);
}
public function fitAndSubtract (h:Vector.<Number>):Gaussian {
var area:Number = 0;
// re-sample to 8 values, find max
var h8_max:Number = 0, p_max:int = 0;
for (var p:int = 0; p < 8; p++) {
var h8:Number = 0;
for (var q:int = 0; q < 32; q++) {
h8 += h [32 * p + q];
}
if (h8_max < h8) {
h8_max = h8; p_max = p;
}
area += h8;
}
// seed our lloyd-like process with p_max-based values
var g:Gaussian = new Gaussian;
g.area = area;
g.mean = 32 * p_max + 16; g.sigma2 = 32;
for (var i:int = 0; i < 4; i++) {
// calculate raw moments
var m1:Number = 0, m2:Number = 0, n:Number = 1;
var from:int = Math.max (0, g.mean - 2 * g.sigma2);
var to:int = Math.min (256, g.mean + 2 * g.sigma2);
for (var j:int = from; j < to; j++) {
// try to filter distant outliers
var h_ideal:Number = g.f (j);
var coef:Number = 0.2 + 0.8 * Math.exp ( -Math.abs (j - g.mean) / g.sigma2);
if (Math.abs (h [j] - h_ideal) > coef * h_ideal) continue;
var j_hj:Number = j * h [j];
m1 += j_hj;
m2 += j_hj * j;
n += h [j];
}
m1 /= n; m2 /= n;
// re-estimate
g.mean = m1; g.sigma2 = m2 - m1 * m1;
}
// oh, and subtract g.f(j) from h[j]
for (j = 0; j < 256; j++) {
h [j] = Math.max (0, h [j] - g.f (j));
}
return g;
}
public function rectStart (e:MouseEvent):void {
rectangle.x = Math.min (320, Math.max (0, mouseX)); rectangle.width = 0;
rectangle.y = Math.min (240, Math.max (0, mouseY)); rectangle.height = 0;
stage.addEventListener (MouseEvent.MOUSE_MOVE, rectDrag);
stage.addEventListener (MouseEvent.MOUSE_UP, rectEnd);
}
public function rectDrag (e:MouseEvent):void {
rectangle.right = Math.min (320, Math.max (0, mouseX));
rectangle.bottom = Math.min (240, Math.max (0, mouseY));
}
public function rectEnd (e:MouseEvent):void {
stage.removeEventListener (MouseEvent.MOUSE_MOVE, rectDrag);
stage.removeEventListener (MouseEvent.MOUSE_UP, rectEnd);
rectangle = getPositiveRect (rectangle);
}
public function getPositiveRect (rectangle:Rectangle):Rectangle {
var r:Rectangle = rectangle;
if ((r.width < 0) || (r.height < 0)) {
r = new Rectangle (
Math.min (rectangle.x, rectangle.right),
Math.min (rectangle.y, rectangle.bottom),
Math.abs (rectangle.width),
Math.abs (rectangle.height)
);
}
return r;
}
}
}
class Gaussian {
public var area:Number;
public var mean:Number;
public var sigma2:Number;
public function f (i:int):Number {
var a:Number = i - mean, b:Number = 2 * sigma2;
return area * Math.exp ( -a * a / b) / Math.sqrt (b * Math.PI);
}
}