PicoGA: Evolving a 3D Projection Matrix

In my previous blog post PicoGA  – a tiny GA for C# I talk about PicoGA and I demonstrate a few problems that can be solved using it. One of the examples is evolving a 3D Projection Matrix from a few pixel to world coordinate mappings taken from a given photograph.

The example didn’t use “height” in the picture because it contained no heights that I could be sure of. In this example, I’m using the heights of the goal posts in the picture below to show that it  can correctly handle World (x,y,z) to Screen(x,y) mapping.

From this picture I’ve identified 9 points that I’ve determined the World(x,y,z) coordinates (the football field is 105×68 meters, the goal is 2.44 meters high and 7.32 meters wide). The coordinates are specified below the picture.

new-pitch-with numbers

I simply use MS Paint (yay!) to determine the pixel coordinate, and I use the given dimensions of the football field (105×68) to determine the world coordinates;

    new WorldToScreenCase(688, 349, 0,0 ,0), // 1
    new WorldToScreenCase(454, 155, 52.5,0,0), // 2
    new WorldToScreenCase(353, 70,  105,0,0), // 3
    new WorldToScreenCase(127, 169, 52.5, 34,0), // 4
    new WorldToScreenCase(155, 411, 0, 34,0), // 5
    new WorldToScreenCase(220, 364, 0, 34 - 7.32*0.5, 2.44), // 6
    new WorldToScreenCase(80,  378, 0, 34+7.32*0.5,2.44), // 7
    new WorldToScreenCase(92,  57,  105, 34 + 7.32*0.5, 2.44), // 8
    new WorldToScreenCase(145, 56,  105, 34-7.32*0.5,2.44), // 9

It takes a while, here’s a video of evolution in action;

but the program is successful in evolving an appropriate matrix;

image

As you can see, the evolved matrix (the red dots) hits the measured points (blue dots) with quite high accuracy. This matrix can now be used for animating something moving across the football field!

As I stated in my previous post, there are analytical methods for doing this, but I’ve been unable to locate any, and this is a cool use for GA…

This is the code that I use to map from world coordinates to screen coordinates, using my evolved matrix;

private static Point3D ProjectCase(Point3D point, List<double> l)
{
	Point3D p = new Point3D(
		point.X - l[12] * 100,
		point.Y - l[13] * 100,
		point.Z - l[14] * 100);

	Point3D result = new Point3D(
		l[0] * p.X + l[1] * p.Y + l[2] * p.Z + l[3],
		l[4] * p.X + l[5] * p.Y + l[6] * p.Z + l[7],
		l[8] * p.X + l[9] * p.Y + l[10] * p.Z + l[11]);

	if (result.Z != 0)
	{
		result.X *= l[15] / result.Z;
		result.Y *= l[15] / result.Z;
	}

	return result;
}

 

You can find the source code among the PicoGA demos.

PicoGA – a tiny GA for C#

Have you ever found that you really want an extremely tiny and simple Genetic Algorithm method for optimizing a set of doubles (the genotype) given a fitness function? Me too!!! So, I wrote one, because it seemed fun.

I’m using it to evolve two toy problems, later I hope to show it with a slightly less toyish problem.

It’s really small (224 lines of code in a single file). It would be very easy to make it smaller – but then it’d be less readable. It would be equally simple to add new functionality, but then it would be larger and would no longer qualify for the name PicoGA. You can find PicoGA here http://picoga.codeplex.com/

What can you use a GA for?

First, two toy problems;

Find three doubles that sum up to five!

Not very difficult, but as I said, it’s a toy problem. First we must specify;

  • a fitness function
  • a population size (how many individuals (I just picked 20 for this simple problem))
  • a genome size / number of genes (3).

That comes out as one line of code;

      GA ga = new GA(
        20, // Individuals
        3, // Genotype length (number of doubles)
        individual => // Fitness function
          Math.Abs(5 - individual.Genotype.Sum(val => val)) 
        );

Simple! Now we need to run that for a couple of generations – but PicoGA takes care of that also. For each generation, we publish the results, if they’ve improved since the previous generation;

private static void SumTo5()
{
  Console.WriteLine("Sum to 5:...");
  GA ga = new GA(
    20, // Individuals
    3, // Genotype length (number of doubles)
    individual => // Fitness function
      Math.Abs(5 - individual.Genotype.Sum(val => val))
    );

  ga.RunEpoch(
    500, // Number of generations to run for
    null, // Action to perform for each generation
    () => // Action to perform once fitness has improved
    {
      Console.WriteLine(
        "Gen {2}: Fit={1}, Genotype={0}",
        string.Join(
          " ",
          ga.BestIndividual.Genotype.Select(val => val.ToString("0.00"))),
        ga.BestIndividual.Fitness.ToString("0.00"),
        ga.CurrentEpochGeneration);
    });

  Console.WriteLine("Sum to 5: done!");
  Console.WriteLine("");
}

 

Voilá, that’s all there is to it! For my test run, at generation 45 a solution is found;

PicoGA, demos

Sum to 5:…

Gen 0: Fit=3,39, Genotype=0,84 0,50 0,26

Gen 1: Fit=3,26, Genotype=0,97 0,50 0,26

Gen 43: Fit=0,04, Genotype=2,45 1,38 1,13

Gen 44: Fit=0,01, Genotype=2,48 1,38 1,13

Gen 45: Fit=0,00, Genotype=2,49 1,38 1,13

Sum to 5: done!

Find 1 2 3 4 5

For my second toy problem, I want the GA to find the values 1 2 3 4 5 (in that exact order) for it’s five genes. Again, not a very useful thing to do, possibly even less useful than my first example, but it’s a toy problem;

private static void Find12345()
{
    Console.WriteLine("Find 1 2 3 4 5:...");
    GA ga = new GA(
        200, // Number of individuals
        5, // Number of genes in the genotype
        individual => // Fitness function
        Math.Abs(individual.Genotype[0] - 1) +
        Math.Abs(individual.Genotype[1] - 2) +
        Math.Abs(individual.Genotype[2] - 3) +
        Math.Abs(individual.Genotype[3] - 4) +
        Math.Abs(individual.Genotype[4] - 5));

    ga.RunEpoch(500, null, () =>
        {
            Console.WriteLine(
                "Gen {2}: Fit={1}, Genotype={0}",
                string.Join(
                " ", 
                ga.BestIndividual.Genotype.Select(
                    val => val.ToString("0.00"))),
                ga.BestIndividual.Fitness.ToString("0.00"),
                ga.CurrentEpochGeneration);
        });

    Console.WriteLine("Find 1 2 3 4 5: done!");
    Console.WriteLine("");
}

 

At generation 105 a perfect solution is found. The reason it takes so long is that genes only mutated slowly and reaching 5 took a while. I’ve since upped mutation rate.

Gen 105: Fit=0,00, Genotype=1,00 2,00 3,00 4,00 5,00

Ok, but what can I really use it for? No more toy problems!

Well, lots and lots of stuff!

I once had an image of a soccer field and I needed a 3D perspective projection matrix that matched the image. I needed to render football players as they were moving around the pitch, for which I had pitch coordinates. But I needed to convert those pitch coordinates to screen coordinates, which isn’t as simple as it sounds.

The 3d Perspective Projection Matrix contains all information about the setup of the camera and target as the image was shot. Camera Field Of View, position of camera relative to subject, orientation of camera with regard to subject etc. For this I used a genotype of 18 genes (in effect 18 doubles).

The data we need

The pitch looked something like this;

image

First I measure points on the I could pinpoint pixels in the image and determine where they where in soccer field space and tie that to image space coordinates.

image

I compiled a list like this (a soccer pitch is 105 x 68 meters)

  1. Pixel 165,178=> World 0, 0, 0
  2. Pixel 377,120=> World 52.5, 0, 0
  3. Pixel 513, 85=> World 105, 0, 0
  4. Pixel 473, 157=>World 52.5, 34, 0
  5. Pixel 611, 208 =>World 52.5, 68, 0
  6. Pixel 735, 136 => World 105, 68, 0

There are analytical methods of retrieving the matrix given only that information, but I don’t know them and I couldn’t find them by googling. But I know GA, and using GA I evolved the matrix!

The fitness function

The fitness function is basically that we feed the six world coordinates into the 3d transformation matrix and sum up the errors of the pixel coordinates it suggests as compared to the pixel coordinates we previously measured. Trivial. Not really, because the formula for transforming from 3D space to 2D screen space isn’t trivial. This is the code I wound up using for this;

private static Point3D ProjectCase(Point3D point, List<double> l)
{
    Point3D p = new Point3D(
        point.X - l[12] * 100,
        point.Y - l[13] * 100,
        point.Z - l[14] * 100);

    Point3D result = new Point3D(
        l[0] * p.X + l[1] * p.Y + l[2] * p.Z + l[3],
        l[4] * p.X + l[5] * p.Y + l[6] * p.Z + l[7],
        l[8] * p.X + l[9] * p.Y + l[10] * p.Z + l[11]);

    if (result.Z != 0)
    {
        result.X *= l[15] / result.Z;
        result.Y *= l[15] / result.Z;
    }

    result.X += l[16] * 100;
    result.Y += l[17] * 100;

    return result;
}

 

Some of the times, PicoGA quickly (in a matter of minutes) finds a matrix that it likes. I’ve plotted in the pixels that the generated matrix produces in the picture below. They’re not perfect, but neither are the measurements that were used for generating the matrix, so they can never be perfect. The red dots are where the evolved matrix suggested;

image

Here’s the program for evolving the matrix;

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Windows.Media.Media3D;

namespace PicoGA
{
    public static class FindMatrix
    {
        internal class WorldToScreenCase
        {
            public WorldToScreenCase(double screenX, double screenY, double worldX, double worldY)
            {
                Screen = new Point3D(screenX, screenY, 0);
                World = new Point3D(worldX, worldY, 0);
            }

            internal Point3D Screen { get; set; }
            internal Point3D World { get; set; }
        }

        public static void FindProjectionMatrix()
        {
            Console.WriteLine("Find ProjectionMatrix:...");
            List<WorldToScreenCase> cases =
                new List<WorldToScreenCase>
                {						
                    // Old Image
                    new WorldToScreenCase(165,178 ,       0,      0),
                    new WorldToScreenCase(377,120 ,       52.5,      0),
                    new WorldToScreenCase(513,85 ,       105,      0),
                    new WorldToScreenCase(473,157 ,       52.5,      34),
                    new WorldToScreenCase(611,208 ,       52.5,      68),
                    new WorldToScreenCase(735,136 ,       105,      68),
                };

            GA ga = new GA(
                2000, // Number of individuals
                18, // Number of genes in the genotype
                individual => // Fitness function
                {
                    double errorSum = 0;
                    foreach (WorldToScreenCase test in cases)
                    {
                        Point3D testScreen = ProjectCase(test.World, individual.Genotype);
                        double sqrError = (testScreen - test.Screen).LengthSquared;
                        errorSum += sqrError;
                    }

                    return errorSum;
                });

            ga.RunEpoch(10000,
                () =>
                {
                    ga.BreakEpochRun = ga.BestIndividual.Fitness <= 1.0 || Console.KeyAvailable;
                },
                () =>
                {
                    Console.WriteLine(
                        "Gen {2}: Fit={1}, Genotype={0}",
                        string.Join(
                        " ",
                        ga.BestIndividual.Genotype.Take(5).Select(
                            val => val.ToString("0.00"))),
                        ga.BestIndividual.Fitness.ToString("0.00"),
                        ga.CurrentEpochGeneration);
                });

            if (Console.KeyAvailable)
            {
                Console.ReadKey();
            }

            Console.WriteLine("Results for training set:");
            foreach (WorldToScreenCase test in cases)
            {
                ShowTestResult(ga, test);
            }

            Console.WriteLine("");
            Console.WriteLine("Additional tests:");
            ShowTestResult(ga, new WorldToScreenCase(120, 73, 105, 34));

            Console.WriteLine("Find ProjectionMatrix: done!");
            Console.WriteLine("");
        }

        private static void ShowTestResult(GA ga, WorldToScreenCase test)
        {
            Point3D p = ProjectCase(test.World, ga.BestIndividual.Genotype);
            Console.WriteLine(
                " World ({0:0.00}, {1:0.00}) => \n" +
                "   Wanted:({2:0.00}, {3:0.00}) \n" +
                "   Received:({4:0.00}, {5:0.00})",
                test.World.X,
                test.World.Y,
                test.Screen.X,
                test.Screen.Y,
                p.X,
                p.Y);
        }

        private static Point3D ProjectCase(Point3D point, List<double> l)
        {
            Point3D p = new Point3D(
                point.X - l[12] * 100,
                point.Y - l[13] * 100,
                point.Z - l[14] * 100);

            Point3D result = new Point3D(
                l[0] * p.X + l[1] * p.Y + l[2] * p.Z + l[3],
                l[4] * p.X + l[5] * p.Y + l[6] * p.Z + l[7],
                l[8] * p.X + l[9] * p.Y + l[10] * p.Z + l[11]);

            if (result.Z != 0)
            {
                result.X *= l[15] / result.Z;
                result.Y *= l[15] / result.Z;
            }

            result.X += l[16] * 100;
            result.Y += l[17] * 100;

            return result;
        }
    }
}

 

You can find PicoGA here http://picoga.codeplex.com/.

Genetic Programming in C#, LambdaGp

I’ve been a fan of Artificial Intelligence for a long time and I’ve been worked both with Neural Networks and Genetic Programming. Since I’ve migrated to C# I’ve lost access to all my old tools, so I decided to start work on a C# based Genetic Programming engine. I call it LambdaGp, because it makes extensive use of Lambda expressions.

After about four hours of coding (this is my fourth GP engine, so I do have some ideas about the design from the start) I have an engine that can solve my first demo problems.

I hope to publish my LambdaGp engine on CodePlex sometime in the future, but here are my two first test problems.

Symbolic regression of a constant

Applying symbolic regression to find a function that evolves to a constant may seem like overkill, but it’s a time-honored “Hello World” test-application in GP circles. Well, I always use it anyway. I’ve stated the problem as; find an expression that evaluates to a specific value, using the functions +,-,/,- and the constant terminals 1, 3 and 5

        private static Population SetupFindConstantTest(double targetValue)
        {
            Population population = new Population();
            population.GenomeDefinition.AddFunctions(
                new Add(),
                new Mul(),
                new Div(),
                new Sub(),
                new Constant(1),
                new Constant(3),
                new Constant(5));
            population.GenomeDefinition.MaxDepth = 3;
            population.PopulationSize = 750;

            population.FitnessFunction =
                ind =>
                Math.Abs(targetValue - ind.Evaluate());
            return population;
        }

Note how the fitness function simply returns the absolute distance between what the expression evaluates to, and the sought value. A perfect hit would return zero.

Below I’ll show how I use the population generated above, but when I ask it to find an expression that evaluates to 13,5 it goes through these iterations;

Generation=1,
New champion:
(-
  (/
    (* 5 3)
    (/ 3 3))
  (+
    (/ 5 5)
    (/ 3 5))),
Fitness=0,25
.
.
.
Generation=43,
New champion:
(-
  (* 3 5)
  (/ 3
    (- 5 3))),
Fitness=0,09
Reached target fitness!
Done

In other words, 3*5-3/(5-3). The reason that the fitness isn’t zero (lower fitness is better) is because of “parsimony pressure”; larger expressions have lower fitness. Fitness is punished depending on the size of the found expression.

Symbolic regression of an actual function

My second demo problem is finding an expression that matches a lambda function, for instance (5 * x^2) + (x * 3) + 5.

You’ll notice how I’ve used lambda for the fitness function instead of using a delegate. I think it’s readable enough, but not everyone would agree. But I could easily use a method delegate instead.

The fitness function iterates from –2 * pi to + 2 * pi in steps of 0.3, and calculates the expected value and the value that the expression suggests. If then keeps a running sum of squared errors. Again, a perfect solution would return zero. Anyway, this method struggled like crazy – to find my equation. Turns out I forgot to add multiplication, a fairly useful function…

The Ephemeral function is like a constant, only it randomizes new constants every time it’s used, in the range –5 to +5.

With multiplication, the new C# code looks like this;

private static Population SetupFindFunctionTest(Func<double, double> func)
        {
            Population population = new Population();
            population.GenomeDefinition.MaxSize = 30;
            population.GenomeDefinition.MaxDepth = 5;
            population.GenomeDefinition.AddFunctions(
                new Add(),
                new Div(),
                new Sub(),
                new Mul(),
                new Pow(),
                new Ephemeral(-5, +5));

            Variable x = population.GenomeDefinition.AddFunction(new Variable("x", 0.1));
            population.GenomeDefinition.MaxDepth = 5;
            population.PopulationSize = 2000;
            population.MaxGenerations = 500;

            population.FitnessFunction =
                ind =>
                {
                    x.Value = -2 * Math.PI;
                    double errorSum = 0;

                    while (x.Value < 2 * Math.PI)
                    {
                        double targetValue = func(x.Value);
                        double error = targetValue - ind.Evaluate();
                        errorSum += error * error;

                        x.Value += 0.3;
                    }

                    return errorSum;
                };

            return population;
        }

And it comes up with this code;

Generation=309,
New champion:
(+ x
  (+ 0,019
    (+
      (+
        (+ x x)
        (+ 4,977
          (*
            (+ x x)
            (+ x x))))
      (* x x)))),
Fitness=0,210671999999998

Again, the solution is actually “perfect”, it’s only parsimony pressure that gives it a non-zero fitness.

Symbolic regression of Sinus

To give it a trickier problem, I’m going have it try find sinus. This a complicated function to solve by symbolic regression – the Taylor expansion is infinite and the first iterations are quite significant.

sin(x) = X - X3/ 3! + X5/ 5! - ... + (-1)(n+1) * X(2*n-1)/ (2n-1)!

This proves to be much harder, as expected. I re-wrote the fitness function to report average squared error instead of the summed squared error, this gives a better feel for the actual result. I found the Taylor expansion for sinus here, and on that page you can also find the information that the fifth power expansion is accurate between –pi/2 to +pi/2 – and since we’re testing from –pi*2 to +pi*2, we need a solution of a higher power than five!

After a while, it closes in and finally finds this solution

Generation=516,
New champion:
(/
  (+ -10,02609
    (* x x))
  (-
    (-
      (/ x -4,28915)
      (*
        (/ -4,64585 x)
        (+ 1,84747 -4,15677)))
    (/
      (/
        (*
          (* x x)
          (/ x -6,95817))
        (/ -4,73401 x))
      (+
        (/ x -4,54048)
        (*
          (/ -3,55989 x)
          (* 1,65336 -1,95395)))))),
Fitness=0,00918493349840388 (before parsimony 0,00918493349840388)

This can of course be further reduced, but I’ll leave that as an exercise for the reader…

For completeness, the code to execute population run looks like this;

public static void Main(string[] args)
        {
            Population pop =
                //SetupFindConstantTest(13.5);
                //SetupFindFunctionTest(x => (x * x * 5) + (x * 3) + 5);
                SetupFindFunctionTest(Math.Sin);

            pop.AfterGeneration +=
                (sender, e) =>
                {
                    if (e.DifferentFromPreviousGeneration)
                    {
                        Console.WriteLine(
                            "Generation={0}, \n\rNew champion: \n\r{1}, \n\rFitness={2} (before parsimony {3})",
                            e.GenerationCount,
                            e.Champion.ToPrettyString(),
                            e.Champion.Fitness,
                            e.Champion.FitnessBeforeParsimony);
                    }
                    if (e.TargetFitnessReached)
                    {
                        Console.WriteLine("Reached target fitness!");
                    }
                };

            while (pop.GenerationCount < pop.MaxGenerations && pop.TargetFitnessReached == false)
            {
                pop.RunGeneration();
                if (Console.KeyAvailable)
                {
                    ConsoleKeyInfo key = Console.ReadKey();
                    if (key.KeyChar == 'R' || key.KeyChar == 'r')
                    {
                        pop.ResetPopulation();
                    }
                }
            }

            Console.WriteLine("Done");
            Console.ReadKey();
        }