Better Best Fit

As part of my code tidy up I have been thinking again about fitting a line to a set of points.
The new code uses the method described here

#pragma indent

using System
using System.Console
using System.Math
using System.Windows.Forms
using Nemerle.Utility
using Nemerle.Assertions

[Record] struct Point
    [Accessor] _x : double
    [Accessor] _y : double

    public this (x : int, y : int)
        _x = x
        _y = y

    public Offset (x : double, y : double) : Point
        Point(_x + x, _y + y)

    public override ToString () : string
        $"($_x, $_y)" 

module ListPointExts
    private static Sum[T] (items : list[T], f : T -> double) : double
        items.FoldLeft(0.0, (pt, a) => f(pt) + a)

    private static Centroid (points : list[Point]) : Point \
        requires points.Length > 0
        Point(Sum(points, pt => pt.X) / points.Length, Sum(points, pt => pt.Y) / points.Length)

    public static BestFit (this points : list[Point]) : (Point * double)
        def centroid = Centroid(points)
        def pts = points.Map(pt => pt.Offset(-centroid.X, -centroid.Y))

        def xySum = Sum(pts, pt => pt.X * pt.Y)
        def x2_y2Sum = Sum(pts, pt => pt.X * pt.X - pt.Y * pt.Y)

        if (-double.Epsilon < xySum && xySum < double.Epsilon)
            if (-double.Epsilon < x2_y2Sum && x2_y2Sum < double.Epsilon)
                // Points all evenly spaced around centroid
                // HACK: Return horizontal line for now
                (centroid, 0.0)
            else if (x2_y2Sum < 0)
                // Points lie on vertical line
                (centroid, PI / 2)
            else
                // Points lie on horizontal line
                (centroid, 0.0)
        else
            def b = x2_y2Sum / xySum
            def a1 = Atan( (-b + Sqrt(b*b + 4)) / 2 )
            def a2 = Atan( (-b - Sqrt(b*b + 4)) / 2 )

            def SumDistances(a)        
                def cos_a = Cos(a)
                def sin_a = Sin(a)
                def RotateY2(pt)
                    def y = (-pt.X * sin_a) + (pt.Y * cos_a)
                    y * y
                Sum(pts, RotateY2)

            if (SumDistances(a1) < SumDistances(a2))
                (centroid, a1)
            else
                (centroid, a2)

def f = Form()
mutable pts : list[Point] = []

f.MouseUp += fun (_,e)
    if (e.Button == MouseButtons.Left)
        pts = Point(e.X, e.Y) :: pts
    else
        pts = []
    f.Invalidate()

f.Paint += fun (_,e)
    when (pts.Length > 0)
        def (c, a) = pts.BestFit()

        def x1 = c.X - Cos(a) * 100
        def y1 = c.Y - Sin(a) * 100
        def x2 = c.X + Cos(a) * 100
        def y2 = c.Y + Sin(a) * 100
        e.Graphics.DrawLine(System.Drawing.Pens.Red, x1 :> int, y1:> int, x2:> int, y2:> int)

    foreach (pt in pts)
        e.Graphics.FillRectangle(System.Drawing.Brushes.Black, pt.X :> int, pt.Y :> int,1,1)

Application.Run(f)

For now the function returns the centroid of the points and the angle of the line. One special case is when there is no “best” fit line because all points are evenly spaced around the centroid. When line following we would want to simply ignore the bogus angle and instead put a line through the centroid from the previous segment end point.
So the BestFit function really needs some kind of tagged union (variant in Nemerle) data type to return. I will create this once the line follower is formalized.

: 22 Nov 2006 22:51 |  Tags: Project |  Comments (0) |  Report a problem