forked from: ransac straight line

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

// forked from makc3d's ransac straight line
// http://en.wikipedia.org/wiki/RANSAC
// (params are not optimal, however, sometimes no solution found :)
package  
{

    import flash.display.Sprite;
    import flash.events.MouseEvent;
    import flash.geom.Point;
    
    /**
         * @author http://www.dreamincode.net/forums/user/112276-atik97/
         *             
        //http://www.dreamincode.net/code/snippet2911.htm
        
     * @author Nicolas Barradeau
     * http://en.nicoptere.net
     */
    public class LeastSquares extends Sprite 
    {
        private var points:Vector.<Point> = new Vector.<Point>();
        
        public function LeastSquares() 
        {
            stage.addEventListener( MouseEvent.MOUSE_DOWN, reset );
        }
        
        private function reset(e:MouseEvent):void 
        {
            graphics.clear();
            graphics.lineStyle( 0 );

            var p:Point = new Point( mouseX, mouseY );
            points.push( p );
            for each( p in points )
            {
                graphics.drawCircle( p.x, p.y, 2 );
            }
            
                        var x:int = stage.stageWidth;
                        
            p = randsac( points );//least squares computation
            
            graphics.lineStyle( 0, 0xFF0000 );
            graphics.moveTo( 0, p.x * 0 + p.y );
            graphics.lineTo( x, p.x * x + p.y );
            
        }

        private function randsac ( data:Vector.<Point> ):Point
        {
            var iterations:int = 0;
            var best_model:Point = null;
            var best_consensus_set:Vector.<Point> = new Vector.<Point>;
            var best_error:Number = Number.MAX_VALUE;
            
            while (iterations++ < 100) {
                // n randomly selected values from data
                var maybe_inliers:Vector.<Point> = new Vector.<Point>;
                for (var n:int = Math.min (data.length, 4); n > 0; n--)
                    maybe_inliers.push (data [Math.round (Math.random () * (data.length - 1))]);

                var maybe_model:Point = leastSquares (maybe_inliers);
                var consensus_set:Vector.<Point> = maybe_inliers.slice ();

                for each (var point:Point in data)
                if (maybe_inliers.indexOf (point) < 0) {
                    var error:Number = Math.abs (maybe_model.x * point.x + maybe_model.y - point.y);
                    if (error < 50) consensus_set.push (point);
                }

                if (consensus_set.length >= Math.min (data.length, 10)) {
                    // this implies that we may have found a good model
                    // now test how good it is
                    var better_model:Point = leastSquares (consensus_set);
                    var this_error:Number = 0;
                    for each (point in consensus_set)
                        this_error += Math.abs (better_model.x * point.x + better_model.y - point.y); this_error /= consensus_set.length;

                    if (this_error < best_error) {
                        // we have found a model which is better than any of the previous ones,
                        // keep it until a better one is found
                        best_model = better_model;
                        best_consensus_set = consensus_set;
                        best_error = this_error;
                    }
                }
            }

            return best_model ? best_model : leastSquares (data);
        }
        
        private function leastSquares( points:Vector.<Point> ):Point
        {
            
            var n:int = points.length;
            var sum_x:Number = 0, sum_y:Number = 0, sum_xx:Number = 0, sum_xy:Number = 0;
            
            for each( var p:Point in points )
            {
                sum_x += p.x;
                sum_y += p.y;
                sum_xx += ( p.x * p.x );
                sum_xy += ( p.x * p.y );
            }
            
            var a:Number = ( -sum_x * sum_xy + sum_xx * sum_y ) / ( n * sum_xx - sum_x * sum_x );
            var b:Number = ( -sum_x * sum_y + n * sum_xy ) / ( n * sum_xx - sum_x * sum_x );
            
            return new Point( b, a );
        }
        
    }
}