/* (C) 2012 Claudio Rocchini - CC BY 3.0 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <vector>
#include <algorithm>

const double PI = 3.1415926535897932384626433832795;
inline double drand() { return double(rand())/RAND_MAX; }

class point2 {
public:
	double x,y;
	point2() {}
	point2( double nx, double ny ): x(nx),y(ny) {}
};

class seed {
public:
	point2 p;		// seed position
	double a;		// main angle
	double l[2];	// extents
};

class inter {		// intersection
public:
	int i,j;
	double ti; double tj;
	inline inter();
	inline inter( int ni, int nj, double nti, double ntj ) : i(ni),j(nj),ti(nti),tj(ntj) {}
	inline bool operator< ( const inter & ii ) const { return fabs(ti)<fabs(ii.ti); }
};

bool intersects( const seed & s1, const seed & s2, double & t1, double & t2 ) {
	double dx1 = cos(s1.a); double dy1 = sin(s1.a);
	double dx2 = cos(s2.a); double dy2 = sin(s2.a);
	const double de = dx1*dy2-dy1*dx2;
	if(fabs(de)<1e-8) return false;
	double ox1 = s1.p.x; double oy1 = s1.p.y;
	double ox2 = s2.p.x; double oy2 = s2.p.y;
	t1 = (-ox1*dy2+ox2*dy2+dx2*oy1-dx2*oy2)/de;
	t2 = ( dx1*oy1-dx1*oy2-dy1*ox1+dy1*ox2)/de;
	return t1<s1.l[0] && t1>-s1.l[1] && t2<s2.l[0] && t2>-s2.l[1];
}

int main() {
	const double SX = 800; const double SY = 800; const int N = 256;
	int i,j;
	static seed seeds[N];
	std::vector<inter> inters; std::vector<inter>::iterator ii;
	
	for(i=0;i<N;++i) {
		seeds[i].p.x = drand()*SX; seeds[i].p.y = drand()*SY;
		seeds[i].a   = PI*drand(); // ortho version: drand()>= 0.5 ? 0 : PI/2; 
		seeds[i].l[0] = SX*10; seeds[i].l[1] = SX*10;
	}
	for(i=0;i<N-1;++i) for(j=i+1;j<N;++j) {
		double t1,t2;
		if(intersects(seeds[i],seeds[j],t1,t2)) {
			inters.push_back( inter(i,j,t1,t2) );
			inters.push_back( inter(j,i,t2,t1) );
	}	}
	std::sort(inters.begin(),inters.end());
	for(ii=inters.begin();ii!=inters.end();++ii) {	// compute intersection on time
		if(ii->ti < seeds[ii->i].l[0] && ii->ti > -seeds[ii->i].l[1] &&
	       ii->tj < seeds[ii->j].l[0] && ii->tj > -seeds[ii->j].l[1]) {
			if( fabs(ii->ti) > fabs(ii->tj) ) {
				if(ii->ti>0) seeds[ii->i].l[0] =  ii->ti;
				else         seeds[ii->i].l[1] = -ii->ti;
	}   }   }
	
	for(i=0;i<N;++i) {	// page clipping
		seed s;
		s.p.x = 0; s.p.y = 0; s.a = 0; s.l[0] = SX; s.l[1] = 0;
		double t1,t2;
		if(intersects(seeds[i],s,t1,t2)) { if(t1>0) seeds[i].l[0] = t1; else seeds[i].l[1] = -t1; }
		s.p.y = SY;
		if(intersects(seeds[i],s,t1,t2)) { if(t1>0) seeds[i].l[0] = t1; else seeds[i].l[1] = -t1; }
		s.p.x = 0; s.p.y = 0; s.l[0] = SY; s.l[1] = 0; s.a = PI/2;
		if(intersects(seeds[i],s,t1,t2)) { if(t1>0) seeds[i].l[0] = t1; else seeds[i].l[1] = -t1; }
		s.p.x = SX;
		if(intersects(seeds[i],s,t1,t2)) { if(t1>0) seeds[i].l[0] = t1; else seeds[i].l[1] = -t1; }
	}
	FILE * fo = fopen("gilbert.svg","w");	// output the results
	fprintf(fo,
		"<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n"
		"<svg xmlns=\"http://www.w3.org/2000/svg\" version=\"1.0\" width=\"%g\" height=\"%g\" id=\"%s\">\n"
		,SX,SY,"girlbert"
	);
	fprintf(fo,"<g style=\"stroke:#000000;stroke-width:2;stroke-linecap:round;fill:none\">\n");
	for(i=0;i<N;++i){
		fprintf(fo,"<line x1=\"%6.2f\" y1=\"%6.2f\" x2=\"%6.2f\" y2=\"%6.2f\"/>\n"
			,seeds[i].p.x+seeds[i].l[0]*cos(seeds[i].a)
			,seeds[i].p.y+seeds[i].l[0]*sin(seeds[i].a)
			,seeds[i].p.x+seeds[i].l[1]*cos(seeds[i].a+PI)
			,seeds[i].p.y+seeds[i].l[1]*sin(seeds[i].a+PI)
		);
	}
	fprintf(fo,"</g>\n");
	fprintf(fo,"</svg>\n");
	fclose(fo);
	return 0;
}