Fitting gaussian to image histogram

/**
 * Copyright makc3d ( http://wonderfl.net/user/makc3d )
 * MIT License ( http://www.opensource.org/licenses/mit-license.php )
 * Downloaded from: http://wonderfl.net/c/3NvK
 */

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.
	 * @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 h:Vector.<Number> = frame.histogram (getPositiveRect (rectangle)) [0];
			// try lloyd-like process
			var mean:Number = 128, sigma2:Number = 64;
			for (var i:int = 0; i < 7; i++) {
				// calculate raw moments
				var m1:Number = 0, m2:Number = 0, n:Number = 1;
				var from:int = Math.max (0, mean - 2 * sigma2);
				var to:int = Math.min (256, mean + 2 * sigma2);
				for (var j:int = from; j < to; j++) {
					var j_hj:Number = j * h [j];
					m1 += j_hj;
					m2 += j_hj * j;
					n += h [j];
				}
				m1 /= n; m2 /= n;
				// re-estimate
				mean = m1; sigma2 = m2 - m1 * m1;
			}
			// find max
			var h_max:Number = 0;
			for (i = 0; i < 256; i++) {
				if (h [i] > h_max) h_max = h [i];
			}
			// graph
			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 ();
			graphics.lineStyle (1, 0xFF0000);
			var area:Number = Math.abs (rectangle.width * rectangle.height);
			for (i = 0; i < 256; i++) {
				var f:Number = area * Math.exp ( -(i - mean) * (i - mean) / (2 * sigma2)) / Math.sqrt (2 * Math.PI * sigma2);
				v = 200 * f / 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 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;
		}
	}
}