November 21, 2006

Single point best fit

A while ago I implemented a basic linear regression algorithm to calculate the line of best fit through a set of points. For this stage in prototyping I needed to see what the best fit line looked like at different places on a test image. The “get valid points in circle” function from before is used to get the set of pixels to work with.
An initial use of my regression function highlighted a glaring problem. A set points forming a near vertical “bar” will produce a regression line almost orthogonal to them! This is because a normal linear regression minimizes the Y-error only. (The assumption is that X values are accurate and Y values are experimental results.)
However, I needed both X and Y errors minimized. I did some quick research and it is possible to perform a “perpendicular” regression. This minimizes the perpendicular distance of points to the line.
However, in the interests of getting something working fast, I decided to calculate two regular regression lines. The first as before, the second with X and Y values interchanged. These two lines represent possible best fits. To find the best I calculate the sums of perpendicular distances to each line from every point; taking the smallest sum as best.

Now that I’ve actually written all that down, maybe doing a “perpendicular regression” would be cleaner! (http://www.mathpages.com/home/kmath110.htm has a good description of the process.) For now though the code does produce good looking lines on top of my test image.

Here is the code (omitting the using statements and GetColoredPointsInCircle function from before).

def DistanceBetween(pt : Point, gradient : double, intercept : double)
    def orthogGradient = -1 / gradient

    if (double.IsInfinity(gradient))
        Abs(pt.X - intercept)
    else if (double.IsInfinity(orthogGradient))
        Abs(pt.Y - intercept)
    else
        def orthogIntercept = pt.Y - orthogGradient * pt.X

        def a = gradient
        def b = 1
        def c = orthogGradient
        def d = 1
        def detInv = 1 / (a * d - b * c)
        def aInv = detInv * d
        def bInv = detInv * -b
        def cInv = detInv * -c
        def dInv = detInv * a

        def x = aInv * intercept + bInv * orthogIntercept
        def y = cInv * intercept + dInv * orthogIntercept

        def xDiff = x + pt.X
        def yDiff = y - pt.Y
        Sqrt(xDiff * xDiff + yDiff * yDiff)

def GetLine(points : list[Point]) : (double * int * int)
    def SumWith(f) { points.Map(f).FoldLeft(0, (x,y) => x + y) }

    def xSum : long = SumWith(pt => pt.X)
    def ySum : long = SumWith(pt => pt.Y)
    def xySum : long = SumWith(pt => pt.X * pt.Y)
    def xxSum : long = SumWith(pt => pt.X * pt.X)
    def yySum : long = SumWith(pt => pt.Y * pt.Y)

    def oX = (xSum / points.Length) :> int
    def oY = (ySum / points.Length) :> int

    def denomX : double = (xSum * xSum - points.Length * xxSum)
    def denomY : double = (ySum * ySum - points.Length * yySum)

    def useX = (denomX < -double.Epsilon || double.Epsilon < denomX)
    def useY = (denomY < -double.Epsilon || double.Epsilon < denomY)

    if (useX && useY)
        def gradientX = (xSum * ySum - points.Length * xySum) / denomX
        def interceptX = (xSum * xySum - ySum * xxSum) / denomX
        def gradientY = 1 / ((ySum * xSum - points.Length * xySum) / denomY)
        def interceptY = \
            if(double.IsInfinity(gradientY)) \
                xSum / (points.Length :> double) \
            else \
                -gradientY * (ySum * xySum - xSum * yySum) / denomY

        def countX = points.FoldLeft(0.0, (pt,a) => a + DistanceBetween(pt, gradientX, interceptX))
        def countY = points.FoldLeft(0.0, (pt,a) => a + DistanceBetween(pt, gradientY, interceptY))

        if (countX < countY)
            (gradientX, oX, oY)
        else
            (gradientY, oX, oY)
    else if (useX)
        ((xSum * ySum - points.Length * xySum) / denomX, oX, oY)
    else // useY
        ((ySum * xSum - points.Length * xySum) / denomY, oX, oY)

def f = Form()
def bmp = Bitmap("line.png")

f.Paint += fun(_, e)
    e.Graphics.DrawImage(bmp, 0,0)

f.MouseUp += fun(_, e)
    try
        def radius = 10
        def points = GetColoredPointsInCircle(bmp, Point(e.X, e.Y), radius, Color.Black)
        def (gradient, oX, oY) = GetLine(points)
        using (g = f.CreateGraphics())
            def angle = Atan(gradient)
            def x0 = (0.5 + Cos(angle) * radius + oX) :> int
            def x1 = (0.5 + Cos(angle + PI) * radius + oX) :> int
            def y0 = (0.5 + Sin(angle) * radius + oY) :> int
            def y1 = (0.5 + Sin(angle + PI) * radius + oY) :> int
            g.DrawLine(Pens.Yellow, x0, y0, x1, y1)
    catch
        | ex => Console.WriteLine(ex.ToString())

Application.Run(f)

The amount of special case handling of vertical lines (infinite gradients) does look ugly. That is definately something to look at refactoring later.

Here is an example of the output when clicking on the image in a few places:
Best Fit Demo

It looks pretty good! The only problem to note is when we are inside a huge mass of points (viz. top right of image). Here there is no way to correctly identify a line. In this special case we would probably just continue forward on our current heading.

The next step is to join up these line segments somehow…


- No comments Not publicly viewable


Add a comment

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

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
© MMXXII