path: root/graytrace.cpp

#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include <string>
#include <sstream>
#include <cmath>
#include <cassert>
#include "lodepng.h"

#define SCREEN_WIDTH (960)
#define SCREEN_HEIGHT (640)

#define FLOAT_SMALL_NUM (1e-8)

using namespace std;

/*
Scene description file (.rtsc, Ray Trace SCene)
===============================================

Z           The positive Y-axis points away from you in this picture.
^  _Y
|  /'
| /
|/
#------> X

A "number" is a floating-point number, unless stated otherwise.
A "coordinate" is three numbers, separated by a space.
Between numbers, any amount of whitespace is tolerated.

First line: coordinate of camera, the coordinate of point the camera looks at
  directly and a number in [0,2π) indicating the amount the camera is rotated
  clockwise around the axis between the camera position and the point the camera
  looks at. Z normally faces up in this right-handed coordinate system.
Second line: coordinate of the point light source, followed by the intensity in
  [0,1].
Further lines: a string in {triangle,sphere} giving the object shape, a
  coordinate giving the centre of the object and further parameters for the
  object. For a sphere, give the size (radius); for a triangle give the
  coordinates of the two other vertices; for a plane, give the normal the
  position is any point on the plane.

Example file:
1 -15 2  0 0 0  0
-2 -2 6  1
sphere 0 0 0  1
sphere 1 10 -3  1.5
*/

bool samesign(const double a,const double b){
	return (a>0&&b>0)||(a<0&&b<0);
}

class Vec3{
public:
	double x,y,z;
	Vec3(void){x=y=z=0;}
	Vec3(double _x,double _y,double _z){x=_x;y=_y;z=_z;}
	Vec3& operator=(const Vec3 &other){
		x=other.x;
		y=other.y;
		z=other.z;
		return *this;
	}
	Vec3 operator+(const Vec3 other) const{return Vec3(x+other.x,y+other.y,z+other.z);}
	Vec3 operator-(const Vec3 other) const{return Vec3(x-other.x,y-other.y,z-other.z);}
	Vec3 operator*(const double other) const{return Vec3(x*other,y*other,z*other);}
	double norm(void) const{return sqrt(x*x+y*y+z*z);}
	double dot(const Vec3 &other) const{
		return x*other.x+y*other.y+z*other.z;
	}
	Vec3 cross(const Vec3 &other) const{
		return Vec3(y*other.z-z*other.y,z*other.x-x*other.z,x*other.y-y*other.x);
	}
	bool hasnan(void) const{return isnan(x)||isnan(y)||isnan(z);}
	Vec3 normalised(void) const{return *this*(1/norm());}
	void normalise(void){
		double n=norm();
		x/=n; y/=n; z/=n;
	}
	/*The Vec3::rotated[XYZ] functions rotate around that axis, anti-clockwise
	  if looking at the rotation when said axis is pointer towards the observer.
	  Matrices from en.wikipedia.org/wiki/Rotation_matrix#In_three_dimensions */
	Vec3 rotatedX(double theta) const{
		return Vec3(x,y*cos(theta)-z*sin(theta),y*sin(theta)+z*cos(theta));
	}
	Vec3 rotatedY(double theta) const{
		return Vec3(x*cos(theta)+z*sin(theta),y,-x*sin(theta)+z*cos(theta));
	}
	Vec3 rotatedZ(double theta) const{
		return Vec3(x*cos(theta)-y*sin(theta),x*sin(theta)+y*cos(theta),z);
	}
	string str(void) const{
		stringstream ss;
		ss<<'['<<x<<' '<<y<<' '<<z<<']';
		return ss.str();
	}
};

double angleBetween(const Vec3 &a,const Vec3 &b){
	//          /  u . v  \
	//theta=acos|---------|
	//          \|u| * |v|/
	double angle=acos(a.dot(b)/(a.norm()*b.norm()));
	//if(angle<0)cerr<<"angleBetween(Vec3("<<a.x<<','<<a.y<<','<<a.z<<"),Vec3("<<b.x<<','<<b.y<<','<<b.z<<"))="<<angle<<endl;
	return angle;
}

class Camera{
public:
	Vec3 pos,lookat;
	double rot;
	Camera(void){pos=Vec3(0,0,0);lookat=Vec3(1,0,0);rot=0;}
	Camera(Vec3 p,Vec3 l,double r){pos=p;lookat=l;rot=r;}
	Camera(double px,double py,double pz,double lx,double ly,double lz,double r){
		pos=Vec3(px,py,pz);
		lookat=Vec3(lx,ly,lz);
		rot=r;
	}
};

class Primitive{
public:
	bool isLight;
	Primitive(void):isLight(false){}
	virtual ~Primitive(void){}
	/* A line is represented by the formula:
	    X = P + D * t
	   for any t in R, where P is a point on the line, D is the direction vector
	   of the line and X is any point on the line determined by t.
	   This method finds the intersection of the line with this object which is
	   closest to P. It returns Vec3(nan,nan,nan) if no intersections.*/
	virtual Vec3 lineIntersect(const Vec3 &P,const Vec3 &D)=0;
	virtual Vec3 normalat(const Vec3 &P)=0;
};

class Light:public Primitive{
public:
	Vec3 pos;
	double intensity;
	Light(void){isLight=true;}
	Vec3 lineIntersect(const Vec3 &P,const Vec3 &D){
		return Vec3(NAN,NAN,NAN);
	}
	Vec3 normalat(const Vec3 &P){
		return Vec3(NAN,NAN,NAN);
	}
};

class Sphere:public Primitive{
public:
	double size;
	Vec3 pos;
	/* A sphere is represented by the formula:
	    r^2 = abs(X - C)^2
	   where the radius is r, C is the centre of the sphere and X is any
	   resulting point on the sphere's surface.
	   The t of the intersection is given by: (en.wikipedia.org/wiki/Line–sphere_intersection)
	   t=-(D * (P-C)) plusminus sqrt((D * (P-C))^2 - (P-C)^2 + r^2),
	   iff |D|=1, and thus the point X of the intersection by:
	   X = P + D * t.*/
	Vec3 lineIntersect(const Vec3 &P,const Vec3 &D){
		double discr,DPminC;
		Vec3 unitD,PminC;
		unitD=D.normalised(); //to make sure |D|=1
		PminC=P-pos;
		DPminC=unitD.dot(PminC);
		discr=DPminC*DPminC-PminC.dot(PminC)+size*size;
		if(discr<0)return Vec3(NAN,NAN,NAN); //no intersection
		else if(discr<FLOAT_SMALL_NUM)return P+unitD*-DPminC; //one intersection
		else { //two intersections, choose that for which abs(t) is smallest
			if(abs(-DPminC-sqrt(discr))<abs(-DPminC+discr))return P+unitD*(-DPminC-sqrt(discr));
			else return P+unitD*(-DPminC+sqrt(discr));
		}
	}
	//returns unit vector *from origin*, only indicating direction.
	Vec3 normalat(const Vec3 &P){
		return (P-pos).normalised();
	}
};

class Triangle:public Primitive{
private:
	Vec3 normal;
public:
	Vec3 pos0,pos1,pos2;
	Triangle(void):normal(NAN,NAN,NAN){}
	/*Follows the algorithm at http://geomalgorithms.com/a06-_intersect-2.html ;
	  anything that's unclear should be found there!*/
	Vec3 lineIntersect(const Vec3 &P,const Vec3 &D){
		Vec3 P1(P+D); //P0==P
		double ri,si,ti,udotv,udotu,vdotv,wdotu,wdotv;
		Vec3 normal=normalat(pos0);
		ri=normal.dot(P1-P); //denominator actually
		//cerr<<ri<<' ';
		if(abs(ri)<FLOAT_SMALL_NUM)return Vec3(NAN,NAN,NAN); //parallel to plane
		ri=normal.dot(pos0-P)/ri;
		if(ri<0)return Vec3(NAN,NAN,NAN);
		Vec3 intersection=P+D*ri; //intersection with plane, maybe with triangle
		Vec3 u=pos1-pos0,v=pos2-pos0,w=intersection-pos0;
		udotv=u.dot(v); udotu=u.dot(u); vdotv=v.dot(v); wdotu=w.dot(u); wdotv=w.dot(v);
		ti=udotv*udotv-udotu*vdotv; //actually denominator for si and ti
		si=(udotv*wdotv-vdotv*wdotu)/ti;
		if(si<0||si>1)return Vec3(NAN,NAN,NAN); //if si>1, si+ti>1 as well. Saves calculating ti for real
		ti=(udotv*wdotu-udotu*wdotv)/ti;
		if(ti<0||si+ti>1)return Vec3(NAN,NAN,NAN);
		return intersection; //apparently, the ray intersects the triangle as well. Yay!
	}
	Vec3 normalat(const Vec3 &P){
		if(!normal.hasnan())return normal;
		return normal=(pos1-pos0).cross(pos2-pos0).normalised();
	}
};

class Plane:public Primitive{
public:
	Vec3 pos,normal;
	//Also follows the algorithm from http://geomalgorithms.com/a06-_intersect-2.html
	Vec3 lineIntersect(const Vec3 &P,const Vec3 &D){
		Vec3 P1(P+D); //P0==P
		double ri;
		ri=normal.dot(P1-P); //denominator actually
		//cerr<<ri<<' ';
		if(abs(ri)<FLOAT_SMALL_NUM)return Vec3(NAN,NAN,NAN); //if parallel to plane
		ri=normal.dot(pos-P)/ri;
		if(ri<0)return Vec3(NAN,NAN,NAN);
		return P+D*ri;
	}
	Vec3 normalat(const Vec3 &P){
		return normal;
	}
};

/*x,y are screen coordinates in [0,1], FOV is the horizontal field of view in
  radians, N is the number of objects.*/
double calcColour(double x,double y,double FOV,int N,Camera &camera,vector<Primitive*> objs){
	int i,minatidx=0; //=0 for compiler soothing
	double mindist,dist,angle;
	Vec3 intersection,midD,midDrot,minatinter;
	midD=(camera.lookat-camera.pos).normalised();
	mindist=HUGE_VALF;
	//first we rotate midD around Z such that it lies on the XZ plane
	angle=atan2(midD.y,midD.x); //get angle relative to XZ plane, around Z axis (thus, horizontal)
	midDrot=midD.rotatedZ(-angle);
	//then we rotate it around the Y axis (thus, vertically) by the needed amount
	midDrot=midDrot.rotatedY(-FOV*((double)1-2*y*SCREEN_HEIGHT/SCREEN_WIDTH)); //negative in front because it rotates the wrong way round. See axes at top of source and description of Vec3::rotated[XYZ] functions.
	//then we rotate it back around Z to its previous angle
	midDrot=midDrot.rotatedZ(angle);
	//then we apply the necessary horizontal rotation (around Z)
	midDrot=midDrot.rotatedZ(FOV*((double)1-2*x));
	for(i=0;i<N;i++){
		intersection=objs[i]->lineIntersect(camera.pos,midDrot);
		if(!intersection.hasnan()){ //if intersection was found
			//now go test whether it's BEHIND the camera...
			if((intersection.x-camera.pos.x)*midDrot.x<0||(intersection.y-camera.pos.y)*midDrot.y<0||(intersection.z-camera.pos.z)*midDrot.z<0)continue;
			dist=(intersection-camera.pos).norm();
			if(dist<mindist){
				mindist=dist;
				minatidx=i;
				minatinter=intersection;
			}
		}
	}
	if(mindist==HUGE_VALF){
		//cout<<". ";
		return 0;
	}

	if(objs[minatidx]->isLight){ //if intersection was found
		double cross=(camera.pos-intersection).normalised().cross(objs[minatidx]->normalat(intersection)).norm(); //cross==sin(angle between those vectors)
		return 1-cross/6; //this looks reasonable
	}
	//cout<<minatidx<<' ';
	Light *light=NULL;
	for(i=0;i<N;i++)if(objs[i]->isLight){light=(Light*)objs[i];break;}
	if(light==NULL)return 0.1;
	Vec3 DlightToInter=minatinter-light->pos;
	double distLightToInter=DlightToInter.norm();
	for(i=0;i<N;i++){
		if(objs[i]->isLight)continue;
		intersection=objs[i]->lineIntersect(light->pos,DlightToInter);
		dist=abs((intersection-light->pos).norm());
		if(!intersection.hasnan()&&abs((intersection-minatinter).norm())>FLOAT_SMALL_NUM&&dist<distLightToInter){
			//cout<<abs((intersection-minatinter).norm())<<' ';
			break;
		}
	}
	if(i==N){
		//cout<<"N ";
		//return 1;
		//return 1-abs((light->pos-minatinter).normalised().cross((camera.pos-minatinter).normalised()).norm());
		Vec3 cam=(camera.pos-minatinter).normalised(),nor=objs[minatidx]->normalat(minatinter),lig=(light->pos-minatinter).normalised();
		/*double diff=asin(cam.cross(nor).norm())-asin(nor.cross(lig).norm()); //ANCIENT
		if(diff<0)return 0.1;
		else return 0.1+diff*0.9;*/
		//if(abs(factor)>M_PI/2||!samesign(angleBetween(cam,lig)-angleBetween(cam,nor),angleBetween(nor,lig)))return 0.1;
		//else return 0.1+(1-abs(factor))*0.9;
		/*double ligCamAngle=angleBetween(cam,lig);
		double ligNorAngle=angleBetween(nor,lig);
		double camNorAngle=angleBetween(cam,nor);
		angle=abs(ligCamAngle-2*camNorAngle); //how far, in radians, from the perfect reflection angle
		if(angle>M_PI/2)return 0.1;
		else return 0.1+(1-angle/M_PI)*0.9;*/

		/*double factor=angleBetween(cam,nor)-angleBetween(nor,lig); //for spheres
		if(factor>M_PI/2)return 0.5+(double)rand()/RAND_MAX*0.5;
		if(factor<0)return 0.1;
		return 0.1+factor/(M_PI/2)*0.9;*/
		/*double factor=angleBetween(cam,lig)-2*angleBetween(cam,nor); //for cubes
		if(abs(factor)>M_PI)return 0.1;
		return 0.1+(1-abs(factor)/M_PI)*0.9;*/
		Vec3 reflected=nor*cam.dot(nor)*2-cam; // paulbourke.net/geometry/reflected/
		double diff=angleBetween(lig,reflected); //always >=0
		if(diff>M_PI)return 0.1;
		//else return 0.1+0.9*(pow(M_E,-diff)-diff*0.01375541740409698); //y=e^(-x)-x*e^(-pi)/pi gives an exponential thing through (0,1) and (pi,0)
		//else return 0.1+0.9*cos(diff/2);
		//else return 0.1+0.9*sin(0.5*M_PI*nor.dot(lig))*(1-diff/M_PI);
		else return 0.1+0.9*0.5*(sin(M_PI*(nor.dot(lig)-0.5))+1)*(1-diff/M_PI);
	} else {
		//cout<<abs((intersection-minatinter).norm())<<' ';
		return 0.1;
	}
	//return 1/(mindist+1);
	//return (minatinter.z+1)/2;
}

void makeimage(vector<uint8_t> &image,Camera &camera,vector<Primitive*> objs,int N){
	double x,y;
	//double minclr=10000,maxclr=-10000;
	//cerr<<"camera.pos="<<camera.pos.str()<<endl;
	for(y=0;y<SCREEN_HEIGHT;y++){
		for(x=0;x<SCREEN_WIDTH;x++){
			//cout<<".abcdefghijklmnopqrstuvwxyz"[lround(calcColour((double)x/SCREEN_WIDTH,(double)y/SCREEN_HEIGHT,M_PI/4,N,camera,objs)*26)]<<flush;
			//cout<<calcColour((double)x/SCREEN_WIDTH,(double)y/SCREEN_HEIGHT,M_PI/2,N,camera,objs)<<" "<<flush;
			image[4*(SCREEN_WIDTH*y+x)+0]=
			image[4*(SCREEN_WIDTH*y+x)+1]=
			image[4*(SCREEN_WIDTH*y+x)+2]=calcColour((double)x/SCREEN_WIDTH,(double)y/SCREEN_HEIGHT,M_PI/4,N,camera,objs)*255;
			image[4*(SCREEN_WIDTH*y+x)+3]=255;
			//if(image[4*(SCREEN_WIDTH*y+x)+0]/255.0<minclr)minclr=image[4*(SCREEN_WIDTH*y+x)+0]/255.0;
			//if(image[4*(SCREEN_WIDTH*y+x)+0]/255.0>maxclr)maxclr=image[4*(SCREEN_WIDTH*y+x)+0]/255.0;
			//cout<<hex<<setw(2)<<(int)image[4*(SCREEN_WIDTH*y+x)+0]<<' ';
		}
		//cout<<endl;
	}
	//cout<<"colour: ["<<minclr<<','<<maxclr<<']'<<endl;
}

int main(int argc,char **argv){
	if(argc!=2){
		cout<<"Please give the scene description file as a command-line parameter."<<endl;
		return 1;
	}
	ifstream in(argv[1]);
	if(!in.good()){
		cout<<"Could not open file '"<<argv[1]<<"'."<<endl;
		return 1;
	}
	int N,i;
	string type;
	Camera camera;
	vector<Primitive*> objs;
	in>>camera.pos.x>>camera.pos.y>>camera.pos.z>>camera.lookat.x>>camera.lookat.y>>camera.lookat.z>>camera.rot;
	if(camera.rot!=0){
		cout<<"No camera rotation is currently supported; the given value of "<<camera.rot<<" should be 0."<<endl;
		return 1;
	}
	objs.push_back(new Light);
	in>>((Light*)objs[0])->pos.x>>((Light*)objs[0])->pos.y>>((Light*)objs[0])->pos.z>>((Light*)objs[0])->intensity;
	objs[0]->isLight=true;
	N=1;
	for(i=1;;i++){
		in>>type;
		if(!in.good())break;
		if(type=="sphere"){
			objs.push_back(new Sphere);
			in>>((Sphere*)objs[i])->pos.x>>((Sphere*)objs[i])->pos.y>>((Sphere*)objs[i])->pos.z>>((Sphere*)objs[i])->size;
			//((Sphere*)objs[i])->printName();
			N++;
		} else if(type=="triangle"){
			objs.push_back(new Triangle);
			in>>((Triangle*)objs[i])->pos0.x>>((Triangle*)objs[i])->pos0.y>>((Triangle*)objs[i])->pos0.z
			  >>((Triangle*)objs[i])->pos1.x>>((Triangle*)objs[i])->pos1.y>>((Triangle*)objs[i])->pos1.z
			  >>((Triangle*)objs[i])->pos2.x>>((Triangle*)objs[i])->pos2.y>>((Triangle*)objs[i])->pos2.z;
			//((Triangle*)objs[i])->printName();
			N++;
		} else if(type=="cuboid"){
			Vec3 p1,p2;
			in>>p1.x>>p1.y>>p1.z>>p2.x>>p2.y>>p2.z;
			objs.push_back(new Triangle);
			((Triangle*)objs[i])->pos0.x=p1.x; ((Triangle*)objs[i])->pos0.y=p1.y; ((Triangle*)objs[i])->pos0.z=p1.z;
			((Triangle*)objs[i])->pos1.x=p2.x; ((Triangle*)objs[i])->pos1.y=p1.y; ((Triangle*)objs[i])->pos1.z=p1.z;
			((Triangle*)objs[i])->pos2.x=p2.x; ((Triangle*)objs[i])->pos2.y=p1.y; ((Triangle*)objs[i])->pos2.z=p2.z;
			i++; N++; objs.push_back(new Triangle);
			((Triangle*)objs[i])->pos0.x=p1.x; ((Triangle*)objs[i])->pos0.y=p1.y; ((Triangle*)objs[i])->pos0.z=p1.z;
			((Triangle*)objs[i])->pos1.x=p2.x; ((Triangle*)objs[i])->pos1.y=p1.y; ((Triangle*)objs[i])->pos1.z=p2.z;
			((Triangle*)objs[i])->pos2.x=p1.x; ((Triangle*)objs[i])->pos2.y=p1.y; ((Triangle*)objs[i])->pos2.z=p2.z;

			i++; N++; objs.push_back(new Triangle);
			((Triangle*)objs[i])->pos0.x=p1.x; ((Triangle*)objs[i])->pos0.y=p1.y; ((Triangle*)objs[i])->pos0.z=p2.z;
			((Triangle*)objs[i])->pos1.x=p2.x; ((Triangle*)objs[i])->pos1.y=p1.y; ((Triangle*)objs[i])->pos1.z=p2.z;
			((Triangle*)objs[i])->pos2.x=p2.x; ((Triangle*)objs[i])->pos2.y=p2.y; ((Triangle*)objs[i])->pos2.z=p2.z;
			i++; N++; objs.push_back(new Triangle);
			((Triangle*)objs[i])->pos0.x=p1.x; ((Triangle*)objs[i])->pos0.y=p1.y; ((Triangle*)objs[i])->pos0.z=p2.z;
			((Triangle*)objs[i])->pos1.x=p2.x; ((Triangle*)objs[i])->pos1.y=p2.y; ((Triangle*)objs[i])->pos1.z=p2.z;
			((Triangle*)objs[i])->pos2.x=p1.x; ((Triangle*)objs[i])->pos2.y=p2.y; ((Triangle*)objs[i])->pos2.z=p2.z;

			i++; N++; objs.push_back(new Triangle);
			((Triangle*)objs[i])->pos0.x=p2.x; ((Triangle*)objs[i])->pos0.y=p1.y; ((Triangle*)objs[i])->pos0.z=p1.z;
			((Triangle*)objs[i])->pos1.x=p2.x; ((Triangle*)objs[i])->pos1.y=p2.y; ((Triangle*)objs[i])->pos1.z=p1.z;
			((Triangle*)objs[i])->pos2.x=p2.x; ((Triangle*)objs[i])->pos2.y=p2.y; ((Triangle*)objs[i])->pos2.z=p2.z;
			i++; N++; objs.push_back(new Triangle);
			((Triangle*)objs[i])->pos0.x=p2.x; ((Triangle*)objs[i])->pos0.y=p1.y; ((Triangle*)objs[i])->pos0.z=p1.z;
			((Triangle*)objs[i])->pos1.x=p2.x; ((Triangle*)objs[i])->pos1.y=p2.y; ((Triangle*)objs[i])->pos1.z=p2.z;
			((Triangle*)objs[i])->pos2.x=p2.x; ((Triangle*)objs[i])->pos2.y=p1.y; ((Triangle*)objs[i])->pos2.z=p2.z;

			i++; N++; objs.push_back(new Triangle);
			((Triangle*)objs[i])->pos0.x=p1.x; ((Triangle*)objs[i])->pos0.y=p2.y; ((Triangle*)objs[i])->pos0.z=p1.z;
			((Triangle*)objs[i])->pos1.x=p2.x; ((Triangle*)objs[i])->pos1.y=p2.y; ((Triangle*)objs[i])->pos1.z=p1.z;
			((Triangle*)objs[i])->pos2.x=p2.x; ((Triangle*)objs[i])->pos2.y=p2.y; ((Triangle*)objs[i])->pos2.z=p2.z;
			i++; N++; objs.push_back(new Triangle);
			((Triangle*)objs[i])->pos0.x=p1.x; ((Triangle*)objs[i])->pos0.y=p2.y; ((Triangle*)objs[i])->pos0.z=p1.z;
			((Triangle*)objs[i])->pos1.x=p2.x; ((Triangle*)objs[i])->pos1.y=p2.y; ((Triangle*)objs[i])->pos1.z=p2.z;
			((Triangle*)objs[i])->pos2.x=p1.x; ((Triangle*)objs[i])->pos2.y=p2.y; ((Triangle*)objs[i])->pos2.z=p2.z;

			i++; N++; objs.push_back(new Triangle);
			((Triangle*)objs[i])->pos0.x=p1.x; ((Triangle*)objs[i])->pos0.y=p1.y; ((Triangle*)objs[i])->pos0.z=p1.z;
			((Triangle*)objs[i])->pos1.x=p2.x; ((Triangle*)objs[i])->pos1.y=p1.y; ((Triangle*)objs[i])->pos1.z=p1.z;
			((Triangle*)objs[i])->pos2.x=p2.x; ((Triangle*)objs[i])->pos2.y=p2.y; ((Triangle*)objs[i])->pos2.z=p1.z;
			i++; N++; objs.push_back(new Triangle);
			((Triangle*)objs[i])->pos0.x=p1.x; ((Triangle*)objs[i])->pos0.y=p1.y; ((Triangle*)objs[i])->pos0.z=p1.z;
			((Triangle*)objs[i])->pos1.x=p2.x; ((Triangle*)objs[i])->pos1.y=p2.y; ((Triangle*)objs[i])->pos1.z=p1.z;
			((Triangle*)objs[i])->pos2.x=p1.x; ((Triangle*)objs[i])->pos2.y=p2.y; ((Triangle*)objs[i])->pos2.z=p1.z;

			i++; N++; objs.push_back(new Triangle);
			((Triangle*)objs[i])->pos0.x=p1.x; ((Triangle*)objs[i])->pos0.y=p1.y; ((Triangle*)objs[i])->pos0.z=p1.z;
			((Triangle*)objs[i])->pos1.x=p1.x; ((Triangle*)objs[i])->pos1.y=p2.y; ((Triangle*)objs[i])->pos1.z=p1.z;
			((Triangle*)objs[i])->pos2.x=p1.x; ((Triangle*)objs[i])->pos2.y=p2.y; ((Triangle*)objs[i])->pos2.z=p2.z;
			i++; N++; objs.push_back(new Triangle);
			((Triangle*)objs[i])->pos0.x=p1.x; ((Triangle*)objs[i])->pos0.y=p1.y; ((Triangle*)objs[i])->pos0.z=p1.z;
			((Triangle*)objs[i])->pos1.x=p1.x; ((Triangle*)objs[i])->pos1.y=p2.y; ((Triangle*)objs[i])->pos1.z=p2.z;
			((Triangle*)objs[i])->pos2.x=p1.x; ((Triangle*)objs[i])->pos2.y=p1.y; ((Triangle*)objs[i])->pos2.z=p2.z;
			//((Triangle*)objs[i])->printName();
			N++;
		} else if(type=="plane"){
			objs.push_back(new Plane);
			in>>((Plane*)objs[i])->pos.x>>((Plane*)objs[i])->pos.y>>((Plane*)objs[i])->pos.z
			  >>((Plane*)objs[i])->normal.x>>((Plane*)objs[i])->normal.y>>((Plane*)objs[i])->normal.z;
			((Plane*)objs[i])->normal.normalise();
			N++;
		} else {
			cout<<"Unknown object type '"<<type<<"'!"<<endl;
			return 1;
		}
	}
	cout<<"Read objects: "<<N<<" objects"<<endl;

	string fnamebase;
	if(strrchr(argv[1],'/')==NULL)fnamebase=argv[1];
	else fnamebase=strrchr(argv[1],'/');
	vector<uint8_t> image;
	image.resize(4*SCREEN_WIDTH*SCREEN_HEIGHT);
	int error;
#if 0
	char fname[fnamebase.length()+8];
	cout<<"Rendered frames: "<<endl;
	for(i=0;i<36;i++){
		makeimage(image,camera,objs,N);
		cout<<i<<' '<<flush;
		snprintf(fname,sizeof(fname),"%s%02d.png",fnamebase.c_str(),i);
		error=lodepng::encode(fname,image,SCREEN_WIDTH,SCREEN_HEIGHT);
		if(error){
			cout<<"Lodepng encoder error "<<error<<": "<<lodepng_error_text(error)<<endl;
			return 1;
		}
		//camera.lookat=(camera.lookat-camera.pos).rotated(M_PI/9,0)+camera.pos;
		camera.pos=camera.pos.rotatedZ(M_PI/18);
	}
	cout<<endl<<"Written images"<<endl;
#else
	fnamebase+=".png";
	//camera.pos=camera.pos.rotatedZ(M_PI/2);
	makeimage(image,camera,objs,N);
	cout<<"Rendered objects"<<endl;
	error=lodepng::encode(fnamebase.c_str(),image,SCREEN_WIDTH,SCREEN_HEIGHT);
	if(error){
		cout<<"Lodepng encoder error "<<error<<": "<<lodepng_error_text(error)<<endl;
		return 1;
	}
	cout<<"Written image"<<endl;
#endif
	return 0;
}