November 22, 2006

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.


- No comments Not publicly viewable


Add a comment

You are not allowed to comment on this entry as it has restricted commenting permissions.

Trackbacks

Google Ads

Search this blog

Most recent comments

  • I scratched my eye while i was holding some kind of plastic packaging.. Anyways the corner of the pl… by Ercan on this entry
  • About a year ago my contacts that i was wearing, i guess were fautly, because shortly after they wer… by Jon on this entry
  • I got shower gel in my eye 4 and a half days ago, and becasue i rubbed my eyes a lot, i have scratch… by Chris on this entry
  • This website may help http://www.webmd.com/eye–health/tc/Eye–Injuries–Home–Treatment by S on this entry
  • I somehow got dirt, or debris in my eye. The terrible pain sent me in a tailspin. I was afraid of sa… by Bobbi on this entry

Tags

November 2006

Mo Tu We Th Fr Sa Su
Oct |  Today  | Dec
      1 2 3 4 5
6 7 8 9 10 11 12
13 14 15 16 17 18 19
20 21 22 23 24 25 26
27 28 29 30         

Galleries

Blog archive

Loading…
Not signed in
Sign in

Powered by BlogBuilder
© MMXIV