#
All entries for November 2006

## November 29, 2006

### Point Grouping

Another part of my project is vectorising point objects on an image. The code for this is fairly simple, but I think I may want to re-use it to improve the line follower.

I am thinking about using the closest “point” as a hint to the next heading to take when following a line.

## November 23, 2006

### Line Follower

My initial line follower is working! It was tricky working out which end of a line segment, created by best fitting points, to follow each time. The code is still not perfect since tight hair-pin corners confuse it. The main assumption is that a line only changes direction gradually. So the line follower always try to continue in a direction close to the previous heading.Here is a sample of a vectorized line (shown in red):

To improve performance when turning sharp corners I think the algorithm needs to be able to look at the image at different levels. Perhaps it can start with a small radius to capture points and grow it a few times to get a sample of possible headings. This would probably require the erasure of “used” points from the image.

In fact maintaining multiple possible paths could make the algorithm more adaptive. This would make junction handling easy. Both routes are taken, only one has to succeed.

Here is the NextSegment function so far:

```
public NextSegment (prevPoint : Point, prevLine : Line) : (Point * Point * Line)
def points = GetColoredPointsInCircle(prevPoint)
def (centroid, angle) = BestFit(points)
def rCosA = Cos(angle) * _radius
def rSinA = Sin(angle) * _radius
def startPoint = Point(centroid.X - rCosA, centroid.Y - rSinA)
def endPoint = Point(centroid.X + rCosA, centroid.Y + rSinA)
def line = Line(startPoint, endPoint)
def betterStartPoint = match (line.Intersect(prevLine))
| Line.Intersection.Point (pt) =>
if (PointInLineSegment(pt, startPoint, endPoint))
pt
else
prevPoint
| Line.Intersection.Parallel =>
prevPoint
| Line.Intersection.Coincident =>
prevPoint
def (next,l) = if (Abs(prevLine.AngleTo(line)) > PI/2)
(startPoint, Line(line.Point2, line.Point1))
else
(endPoint, line)
_finished = (next.DistanceTo(_end) < _radius)
(betterStartPoint, next, l)
```

The code to choose the correct next point on the polyline is a bit of hack. Basically it has to flip the best fit “line” as well to ensure the heading is maintained. This block will be changed in future versions to include the multiple branching.

### Line–Line Intersection

I decided to represent lines as two distinct points. Two lines can either:- Intersect at a single point
- Be coincident (i.e. identical)
- Be parallel

I created an Intersection function on my Line type. It takes another Line and returns an Intersection value. “Intersection” is a union type representing the above three cases.

Here’s the code:

```
#pragma indent
using System
using System.Console
using System.Math
using Nemerle.Utility
namespace Prototype
static class DoubleExtensions
public static IsTiny(this d : double) : bool
-double.Epsilon < d && d < double.Epsilon
[Record] struct Point
[Accessor] _x : double
[Accessor] _y : double
public override ToString () : string
$"($_x, $_y)"
public static @:( pt : int * int ) : Point
| (x,y) => Point(x,y)
public static @:( pt : double * double ) : Point
| (x,y) => Point(x,y)
[Record] struct Line
[Accessor] _point1 : Point
[Accessor] _point2 : Point
public Intersect (other : Line) : Intersection
def Det (a,b,c,d)
a * d - b * c
def (x1,x2,x3,x4) = (Point1.X, Point2.X, other.Point1.X, other.Point2.X)
def (y1,y2,y3,y4) = (Point1.Y, Point2.Y, other.Point1.Y, other.Point2.Y)
def d1 = Det(x1, y1, x2, y2)
def d2 = Det(x3, y3, x4, y4)
def denom = Det(x1 - x2, y1 - y2, x3 - x4, y3 - y4)
if (denom.IsTiny())
def x = Det(d1, x1 - x2, d2, x3 - x4)
def y = Det(d1, y1 - y2, d2, y3 - y4)
if (x.IsTiny() && y.IsTiny())
Intersection.Coincident()
else
Intersection.Parallel()
else
def x = Det(d1, x1 - x2, d2, x3 - x4) / denom
def y = Det(d1, y1 - y2, d2, y3 - y4) / denom
Intersection.Point(Point(x,y))
public variant Intersection
| Parallel
| Coincident
| Point
pt : Prototype.Point
using Prototype
def l1 = Line((1,1), (10,10))
def l2 = Line((0,1), (12,13))
def l3 = Line((5,1), (1,5))
def l4 = Line((0,0), (12,12))
def foos = [l1.Intersect(l2), l1.Intersect(l3), l1.Intersect(l4)]
foreach (i in foos)
| Line.Intersection.Point(pt) => WriteLine($"Point $pt")
| Line.Intersection.Coincident => WriteLine("Coincident")
| Line.Intersection.Parallel => WriteLine("Parallel")
```

P.S.

Some Nemerle gems:

1) Implicit cast from (int * int) to Point. This lets me use (1,3) in place of Point(1,3) – cuts down on the visual noise when building up lists of points!

2) match inside foreach

```
foreach (x in xs)
| Foo => ...
| Bar => ...
```

very cool!!

## November 22, 2006

### Line representation

What is the most sensible representation of a line? There are a few possibilities including:- Gradient and intercept: (y=mx+c)
- Two distinct points: (x1,y1) & (x2, y2)
- r = x cos(a) + y sin(a) : r = the perpendicular distance to the line. a = the angle r makes with the x-axis.

I suppose the actual operations I need to perform on lines should dictate the best represention to ease calculations. Using gradients gets tricky when dealing with vertical lines. Angles are preferable since they range between -Pi/2 and Pi/2.

The line following algorithm relies on calculating the intersection of lines. It produces start and end points of line segments. Therefore representing lines as two distinct points may be the best approach. Calculating the intersection of two lines represented this way is described here

### 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.

### Line Following Prototype

Using my previous prototypes I created a line following program. Given two initial points that indicate the start point and heading, it is able to follow a line drawn on my sample image.

The key function is:

def NextPoint(prevGradient : double, prevStartPoint : Point, prevEndPoint : Point, radius : int, bmp : Bitmap) : (double * Point * Point)

- previous gradient
- line segment (start and end points)
- pixel capture radius
- bitmap

and returns - gradient of the next line segment
- the improved end point for the previous line segment
- the end point of the new line segment

Repeatedly calling this function passing in the previous data follows the line. The “improved end point” return is the intersection of the previous and current line segments – this make the polyline fit much better around corners. Without that it will overshoot then swing back each time.

I will not post the code here because frankily it’s getting messy! Also I am having some real difficulty expressing vertical lines as Infinite gradients and an X-intercept. It is now time to introduce some proper data structures to represent lines, etc, and better methods of calculating intercept points, etc.

One more quick idea:

Once the line follower has “used” some pixels should they be “removed” to prevent later confusion? I could not remove after one step, since previous pixels are used to shape the current line segment. But the pixels two-steps back should now be out of range. This would also stop the NextPoint function “turning around” at a dead-end.

## 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.

## November 16, 2006

### VM Crash

I fired up my project’s Virtual PC image only to have it start reporting weird memory errors. I restarted the OS, but then it would not boot into windows. I was seeing a brief blue screen error. So I mounted my Windows ISO image and ran “chkdsk /r”. It seemed to fix some errors. I also ran fixboot.

Still the OS would not boot (some error reading the registry now). I ran chkdsk on the host OS just in case – no errors were found.

I then tried an in-place install of Windows… still wouldn’t boot.

DAMN.

I decided to cut my losses. I launched a different Virtual PC with the broken PCs hard drive attached as secondary. I was then able to copy all the important files.

Now I’m going to re-install Windows from scratch… so much for a productive day!