CIMG problem

Hallo,

ein guter Freund aus Amerika bat mich letztens eine Software zu schreiben die aus einem Bild die Ellipsen erkennt. Daraufhin habe ich den Algoritmus von Xi und Ji gefunden, der aber leider in Matlab geschrieben ist. Nach weiteren suchen habe ich dann gefunden, was der Algoritmus macht:

The following algorithm compares every white pixel to every other white pixel.  We’ll call the first Xi and the white pixel compared to Xj. 
- Calculate the distance between Xi and Xj divided by two and store it in ‘a’. If ‘a’ is less than a pre-determined ellipse ‘a’ minimum (30) than move on to the next pixel, same goes for a greater than threshold (45).
- Calculate the angle of the two pixels and store it in alpha.
- Calculate the midpoint and store it in ‘m’.
- Now we once again look at every other white pixel (hence the O(n^3)). We’ll call this one Xk. For each Xk calculate the distance between it and ‘m’ and store that in ‘d’.  Compare it to the ‘b’ threshold of 15 min/40 max, if it fails move onto the next one.
- Next we calculate the angle between Xk and ‘m’ and store it in theta.
- Now for the heart of the algorithm, calculate the ‘b’ value of the ellipse for these points using a bit of algebra against ellipse equations you can arrive at b=(a^2*d^2*sin^2(theta))/(a^2-d^2*cos^2(theta))^.5
- Increment the ‘b’ count for the found ‘b’.  ie if ‘b’ is ‘2’ then record that ‘2’ has shown up once, if ‘2’ appears from another Xk then it has appeared twice etc.
- Once all of the Xk’s have been exhausted find the ‘b’ value that has appeared the most often and compare this value against every other Xi/Xj pair’s highest ‘b’ count and whichever has the largest amount of values that create the most amount of similar ‘b’ values is the best fit ellipse for the image.
- Now that you have the best fit ellipse values you can easily calculate everything needed to draw the ellipse.

Jetzt hab ich mal selber versucht es zu coden, aber ich bekomme irgendwie sehr viele Werte und weiß nicht mehr weiter. Vielleicht könnt ihr mir ja helfen.
Hier mal mein bisheriger Code:

#include <stdio.h>
#define cimg_debug 2
#include <string>
#include "CImg.h"
#include <map>
#include <iostream>
#include <algorithm>
#include <math.h>
#include <vector>
#include <iterator>
class Pixel
	{
	public:
		Pixel(){}
		Pixel(int x, int y): x(x), y(y) {}

		friend std::ostream& operator<<(std::ostream& ostr, const Pixel& p)
		{
			ostr << p.x << " " << p.y << "\n";
			return ostr;
		}

		int GetX(){return x;}
		int GetY(){return y;}

	private:
		int x;
		int y;
	};

	class a
	{
	public:
		a(){}
		a(int xa, int ya, int xb, int yb, double distanz, double Alpha): xa(xa), ya(ya), xb(xb), yb(yb), distanz(distanz), Alpha(Alpha) {}
		int GetXa(){ return xa; }
		int GetYa(){ return ya; }
		int GetXb(){ return xb; }
		int GetYb(){ return yb; }
		int GetDistanz(){ return distanz; }
		int GetWinkel(){ return Alpha; }
	private:
		int xa, ya, xb, yb;
		double distanz, Alpha;
	};
	typedef std::vector<a> aa;
	typedef std::vector<Pixel> W_pixel;

int main( int argc, const char* argv[] )
{ 
	Map map;
	aa A,D;
	W_pixel pixel;
	W_pixel m;
	Weiß white;
	Distance distanz;
	CImg<unsigned char> image("1.jpg");
	int width = image.width();
	int height = image.height();

double distanz2, Alpha, Theta;
	Pixel aaa, b, c, d;
	for(W_pixel::iterator i = pixel.begin(); i != pixel.end() - 1; ++i)
	{
        for(W_pixel::iterator y = i + 1; y != pixel.end(); ++y) 
		{
			aaa = *i;
			b = *y;
			distanz2 = sqrt((double)((aaa.GetX() - b.GetX())*(aaa.GetX() - b.GetX())) + ((aaa.GetY() - b.GetY())*(aaa.GetY() - b.GetY()))); 
			if(distanz2 > 30 && distanz2 < 45)
			{
				Alpha = math::toDegree( std::atan((double)(aaa.GetY() - b.GetY()) / (aaa.GetX() - b.GetX()))); 
				m.push_back(Pixel((0.5*(aaa.GetX() + b.GetX())), (0.5*(aaa.GetY() + b.GetY()))));
				A.push_back(a(aaa.GetX(), aaa.GetY(), b.GetX(), b.GetY(), distanz2, Alpha));
			}
		}
	}
	double distanz3;
	double b;
	for(W_pixel::iterator k = pixel.begin(); k != pixel.end(); ++k)
	{
        for(W_pixel::iterator l = m.begin(); l != m.end(); ++l) 
		{
			c = *k;
			d = *l;
			distanz3 = sqrt((double)((c.GetX() - d.GetX())*(c.GetX() - d.GetX())) + ((c.GetY() - d.GetY())*(c.GetY() - d.GetY()))); 
			if(distanz3 > 15 && distanz3 < 45)
			{
				Theta = math::toDegree( std::atan((double)(c.GetY() - d.GetY()) / (c.GetX() - d.GetX()))); 
				D.push_back(a(c.GetX(), c.GetY(), d.GetX(), d.GetY(), distanz3, Theta));
				for(aa::iterator i = A.begin(); i != A.end() - 1; ++i)
				{
//b?
				}
			}
		}
	}

Danke schon mal!

P.S.: Mir und mein amerikanischer Freund ist bewusst, dass dieses Programm dann sehr lang brauchen kann um Ellipsen zu finden. Dieses Programm soll auch nur Ellipsen aus Schwarz/Weiß Bilder (also es gibt nur Pixel mit 255 (weiß) oder 0 (schwarz)).