#
All entries for Tuesday 21 November 2006

## 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: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…

### Points In Circle

I need to get all valid coloured pixels within a given radius of an origin point. The follow code creates a Form with an image. When clicking the image, all black pixels within a 10 pixel radius of the clicked pixel are coloured red.

Special care is taken not to go off the edge of the image (note the Max and Min function calls).

```
using System
using System.Drawing
using System.Math
using System.Windows.Forms
def GetColoredPointsInCircle(bmp : Bitmap, origin : Point, radius : int, color : Color) : list[Point]
def xLine = $[Max(0, origin.X - radius) .. Min(bmp.Width - 1, origin.X + radius)]
def yLine = $[Max(0, origin.Y - radius) .. Min(bmp.Height - 1, origin.Y + radius)]
def r2 = radius * radius
def InCircle(x, y)
def dX = x - origin.X
def dY = y - origin.Y
(dX * dX + dY * dY) < r2
def CorrectColorAt(x, y)
bmp.GetPixel(x, y).ToArgb() == color.ToArgb()
$[ Point(x, y) | x in xLine, y in yLine, InCircle(x,y) && CorrectColorAt(x,y)]
def f = Form()
def bmp = Bitmap("line.png")
f.Paint += fun(_, e)
e.Graphics.DrawImage(bmp, 0,0)
f.MouseUp += fun(_, e)
def pts = GetColoredPointsInCircle(bmp, Point(e.X, e.Y), 10, Color.Black)
using (g = f.CreateGraphics())
foreach (pt in pts)
g.FillRectangle(Brushes.Red, pt.X, pt.Y, 1, 1)
Application.Run(f)
```

This allows me to visually check that the right pixels would be used when performing a linear regression to approximate the line. The size of the radius, 10 pixels, is arbitrary at the moment. It reflects the width of the line drawn in the image.

Also the choice of using ”<” and not ”<=” to check if a point is in the circle is arbitrary for now.

P.S.

Check out the awesome list comprehensions Nemerle lets me do!

### Getting Back On Track

The project’s progress has somewhat stalled over the past few weeks. I’ve been too busy with work outside of university. Rather than get bogged down panicking about timetables, phases and other management issues, I decided the most productive thing to do is code. Just get something concrete working! A fair amount of my project is about discovering algorithms and refining them. So instead of pretending that I spent hours writing pseudo-code and defining all the mathematics needed, I am using the most expressive medium I know: code!

My current objective is to build a simple line vectorizer. I am building up the code in a very iterative fashion by writing very simple prototypes. These are most often consecutive; the next prototype copies code from the previous.

I think it will be useful to document the evolution of these prototypes so that at the end I know the design decisions I made along the way. So the plan is to dump each prototype in a blog post and jot down some notes about it.