void analyticPipeline( const std::vector<Splat> &splats, ABuffer *abuffer, const Matrix4x4 &objectToCamera, const Matrix4x4 &projection, const Vector4 &viewport, float front, float back )
{
	SS_ASSERT(abuffer);
	SS_ASSERT(subpixels == 1);	//the analytic pipeline supports only 1 sample per pixel

	const int			exptablesize = 256;
	static float		exptable[exptablesize];
	static bool			exptableinitialized = false;
	const float			exptablescale = float(exptablesize)/F;
	const float			oosqrtF = 1.0f / sqrt(F);

	if( !exptableinitialized )
	{
		for(int i=0;i<exptablesize;++i)
		{
			exptable[i] = exp(-0.5f * F * float(i)/float(exptablesize));
		}
	}

	int width = abuffer->getWidth();
	int	height = abuffer->getHeight();
	int	subpixels = abuffer->getSubpixels();

	int scissorxmin = intFloor(viewport.x);
	int scissorymin = intFloor(viewport.y);
	int scissorxmax = intFloor(viewport.z);
	int scissorymax = intFloor(viewport.w);
	if( scissorxmin < 0 ) scissorxmin = 0;
	if( scissorymin < 0 ) scissorymin = 0;
	if( scissorxmax > width ) scissorxmax = width;
	if( scissorymax > height ) scissorymax = height;

	Vector2 vpScale( 0.5f * (viewport.z - viewport.x), 0.5f * (viewport.w - viewport.y) );
	Vector2 vpMid( 0.5f * (viewport.z + viewport.x), 0.5f * (viewport.w + viewport.y) );

	Vector2 oovpScale(1.0f/vpScale.x, 1.0f/vpScale.y);

	//prefilter
	Matrix2x2 Vp( prefiltersize*oovpScale.x*oovpScale.x, 0,
				  0, prefiltersize*oovpScale.y*oovpScale.y );

	Matrix4x4 objectToScreen;
	objectToScreen = objectToCamera;
	objectToScreen *= projection;

	Matrix4x4 screenToObject;
	screenToObject = objectToScreen;
	bool invertible = screenToObject.invert();
	SS_ASSERT(invertible);
	if( screenToObject[2].w != 0.0f ) screenToObject[2] *= 1.0f/screenToObject[2].w;
	else screenToObject[2] = -screenToObject[2];
	Vector4 c = screenToObject[2];

	Vector4 light(10,10,-20,0);
	Matrix4x4 cameraToObject = objectToCamera;
	cameraToObject.invert();
	light *= cameraToObject;
	Vector3 lightInObject(light.x, light.y, light.z);
	lightInObject.normalize();

	//clip plane setup
	Vector4 clipOrigin[6],clipNormal[6];
	int clipPlanes = 6;
	{
		float s = (float)sqrt(2.0)/2.0f;
		clipOrigin[0].set(0,0,0,0);
		clipNormal[0].set(0,0,1,0);

		clipOrigin[1].set(0,0,back,0);
		clipNormal[1].set(0,0,-1,0);

		clipOrigin[2].set(0,0,0,0);
		clipNormal[2].set(-s,0,0,s);

		clipOrigin[3].set(0,0,0,0);
		clipNormal[3].set(s,0,0,s);

		clipOrigin[4].set(0,0,0,0);
		clipNormal[4].set(0,-s,0,s);

		clipOrigin[5].set(0,0,0,0);
		clipNormal[5].set(0,s,0,s);
	}

	int subpixelsTested = 0;
	int subpixelsPassed = 0;

	for(int i=0;i<splats.size();++i)
	{
		Vector3 p0 = splats[i].m_position;
		Vector3 n = splats[i].m_normal;

		//backface cull (in object-space)
		Vector4 pleq(n.x, n.y, n.z, -dot(p0,n));
		int dotSign = extractSign( dot(c, pleq) );
		if( dotSign ) continue;

		//set up the conic matrix
		Matrix2x2 Vr( splats[i].m_size, 0,
					  0,                splats[i].m_size );
		Matrix2x2 VrT = Vr;
		VrT.transpose();
		Matrix2x2 Q = Vr * VrT;
		Q.invert();

		//make the tangent plane (in object-space)
		Vector3 tu = perpendicular( n );
		tu.normalize();
		Vector3 tv = cross(n, tu);

		//transform the tangents and
		// the splat centre to the screen-space
		Vector4 tuh(tu.x, tu.y, tu.z, 0);
		Vector4 tvh(tv.x, tv.y, tv.z, 0);
		Vector4 p0h(p0.x, p0.y, p0.z, 1);
		tuh *= objectToScreen;
		tvh *= objectToScreen;
		p0h *= objectToScreen;

		//viewport culling
		int cp;
		for(cp=0;cp<clipPlanes;cp++) {
			float a = dot(tuh, clipNormal[cp]);
			float b = dot(tvh, clipNormal[cp]);
			float c = dot(p0h, clipNormal[cp]) - dot(clipOrigin[cp], clipNormal[cp]);
			if( a == 0.0f && b == 0.0f ) { //the planes are parallel
				if( c < 0.0f ) break; //behind the clip plane
			} else {
				if( F * (Q[0][0]*b*b - 2.0f*Q[0][1]*a*b + Q[1][1]*a*a) - c*c*Q.det() >= 0.0f )
				{ //the clip plane intersects the splat
					if( cp == 0 ) break; //front plane discards the splat
				} else { //the clip plane doesn't intersect the splat
					if( c < 0.0f ) break; //behind the clip plane
				}
			}
		}
		if( cp < clipPlanes ) continue;	//the splat was behind some of the planes

		//mapping from the splat to the screen-space
		Matrix3x3 M(tuh.x, tuh.y, tuh.w,
				    tvh.x, tvh.y, tvh.w,
				    p0h.x, p0h.y, p0h.w);

		//mapping from the screen-space to the splat
		if( M.det() == 0.0f ) continue; //don't draw if the matrix is singular
		Matrix3x3 Mi = M;
		Mi.invert();

		//discard the splat if the mapping matrix's condition number is too high
		float conditionNumber = norm_inf(M) * norm_inf(Mi);
		if( conditionNumber > conditionThreshold ) continue;

		//affine approximation
		Matrix3x3 Qh( Q[0][0], Q[0][1], 0,
					  Q[1][0], Q[1][1], 0,
					  0,       0,       -F );
		Matrix3x3 MiT = Mi;
		MiT.transpose();
		Matrix3x3 Qhs = Mi * Qh * MiT;

		float delta = Qhs[0][0] * Qhs[1][1] - Qhs[0][1] * Qhs[1][0];
		if( delta <= 0.0f ) continue; //draw only ellipses
		float oodelta = 1.0f/delta;

		float xt = (Qhs[0][1] * Qhs[1][2] - Qhs[1][1] * Qhs[0][2]) * oodelta;
		float yt = (Qhs[0][1] * Qhs[0][2] - Qhs[0][0] * Qhs[1][2]) * oodelta;

		Matrix2x2 Qs( Qhs[0][0], Qhs[0][1],
					  Qhs[1][0], Qhs[1][1] );

		float f = -Qhs[2][2] - Qhs[0][2] * xt - Qhs[1][2] * yt;
		if( f <= 0.0f ) continue;
		Qs *= F / f;

		//make the resampling filter, the code below is equivalent to
		float detQs = Qs.det();
		if( detQs <= 0.0f ) continue;
		Qs[0][0] = (Qs[0][0] + detQs*Vp[1][1]);
		Qs[0][1] = (Qs[0][1] - detQs*Vp[0][1]);
		Qs[1][0] = (Qs[1][0] - detQs*Vp[1][0]);
		Qs[1][1] = (Qs[1][1] + detQs*Vp[0][0]);
		Qs *= detQs/Qs.det();

		//the bounding rectangle
		delta = Qs.det();
		if( delta <= 0.0f ) continue;
		oodelta = 1.0f / delta;
		float bx = sqrt(Qs[1][1] * F * oodelta);
		float by = sqrt(Qs[0][0] * F * oodelta);
		float xmin = xt - bx, xmax = xt + bx;
		float ymin = yt - by, ymax = yt + by;

		//determine the integer bounding rectangle and do scissoring
		xmin = xmin * vpScale.x + vpMid.x;
		xmax = xmax * vpScale.x + vpMid.x;
		ymin = ymin * vpScale.y + vpMid.y;
		ymax = ymax * vpScale.y + vpMid.y;

		int sx = intFloor( xmin ), ex = intCeil( xmax );
		if( sx >= width || ex <= 0 ) continue;
		sx = max2(sx,scissorxmin);
		ex = min2(ex,scissorxmax);
		if( ex - sx == 0 ) continue;

		int sy = intFloor( ymin ), ey = intCeil( ymax );
		if( sy >= height || ey <= 0 ) continue;
		sy = max2(sy,scissorymin);
		ey = min2(ey,scissorymax);
		if( ey - sy == 0 ) continue;

		//depth interpolation setup
		float zmin = p0h.z,zmax = p0h.z;
		float zx = 0.0f,zy = 0.0f;
		float zt = p0h.z;
		float l = Q[0][0]*tuh.z*tuh.z + 2.0f*Q[0][1]*tuh.z*tvh.z + Q[1][1]*tvh.z*tvh.z;
		if( l != 0.0f ) {
			//calculate zmin and zmax
			l = sqrt(max2(0.0f,F/l));
			zmin = p0h.z - l * (tuh.z*tuh.z + tvh.z*tvh.z);
			zmax = p0h.z + l * (tuh.z*tuh.z + tvh.z*tvh.z);
			zt = (zmin + zmax)*0.5f;

			//calculate the depth gradient, not needed for our surface reconstruction
			zx = -tvh.z * M[0][1] + tuh.z * M[1][1];
			zy =  tvh.z * M[0][0] - tuh.z * M[1][0];
			float dz = F*(zy*zy*Qs[0][0] - 2*zy*zx*Qs[0][1] + zx*zx*Qs[1][1]);
			if( dz != 0.0f ) dz = (zmax-zt) * sqrt(max2(0.0f,delta / dz));
			zx *= dz;
			zy *= dz;
		}

		//call the splat shader here
		Vector4 color(0,0,0,1);
		{
			//------------------- template ---------------------------------
			//diffuse
			float d = dot(splats[i].m_normal, lightInObject);
			d = max2(0.5f,d);
			color = splats[i].m_color;
			color *= d;

			//specular
			Vector3 halfway(cameraToObject[3].x, cameraToObject[3].y, cameraToObject[3].z);
			halfway -= splats[i].m_position;
			halfway.normalize();
			halfway += lightInObject;
			halfway.normalize();
			d = dot(halfway, splats[i].m_normal);
			if( d > 0.0f )
			{
				d *= d;
				d *= d;
				d *= d;
				d *= d;
				d *= d;
				color.x += d;
				color.y += d;
				color.z += d;
			}
			color.w = 1.0f;
			//------------------- template ---------------------------------
		}

		if( rasterization )
		{
			//rasterization loop
			ABuffer::Fragment fragment;
			fragment.zmin = zmin;
			fragment.zmax = zmax;
			fragment.dotSign = dotSign;
			fragment.color = color;
			for(int y=sy;y<ey;++y)
			{
				for(int x=sx;x<ex;++x)
				{
					Vector2 xp(((float)x - vpMid.x)*oovpScale.x - xt,
								 ((float)y - vpMid.y)*oovpScale.y - yt);
					subpixelsTested++;
					float radius = Qs[0][0]*xp.x*xp.x + 2.0f*Qs[0][1]*xp.x*xp.y + Qs[1][1]*xp.y*xp.y;
					if( radius < F )
					{
						subpixelsPassed++;
						float z = zt + xp.x * zx + xp.y * zy; //needed, if the surface reconstruction uses z
						fragment.weight = exptable[max2(0,(int)intFloor(radius*exptablescale))];
						fragment.z = z;
						abuffer->appendFragment(x,y,0,fragment);
					}
				}
			}
		}

		if( pointvisualization )
		{
			float w0 = 1.0f / p0h.w;
			float midx = p0h.x * w0;
			float midy = p0h.y * w0;
			glPointSize(3.0f);
			glBegin(GL_POINTS);
			glVertex3f(midx*vpScale.x*figsx, midy*vpScale.y*figsy, -1);
			glEnd();
		}

		if( reconstructionvisualization )
		{
/*			//draw bounding box
			glBegin(GL_LINE_LOOP);
			glVertex3f((xmin-vpMid.x)*figsx, (ymax-vpMid.y)*figsy, -1);
			glVertex3f((xmax-vpMid.x)*figsx, (ymax-vpMid.y)*figsy, -1);
			glVertex3f((xmax-vpMid.x)*figsx, (ymin-vpMid.y)*figsy, -1);
			glVertex3f((xmin-vpMid.x)*figsx, (ymin-vpMid.y)*figsy, -1);
			glEnd();
*/
			//draw splat boundary
			float w0 = 1.0f / p0h.w;
			float midx = p0h.x * w0;
			float midy = p0h.y * w0;

			int segments = 25;
			float a = 0.0f;
			float aa = 2.0f * 3.1415f / float(segments);
			Matrix2x2 J( Mi[0][0], Mi[0][1],
						 Mi[1][0], Mi[1][1] );
			Matrix2x2 Jt = J;
			Jt.transpose();
			J *= Q;
			J *= Jt;
			float sk = sqrt(F);
			glBegin(GL_LINE_LOOP);
			for(int i=0;i<segments;++i)
			{
				Vector2 uv;
				float f,g,t;

				//perspective correct position
				// a*x^2+2*b*x*y+c*y^2 = k * (oow0 + dwx*x + dwy*y)^2
				// sqrt(a*x^2+2*b*x*y+c*y^2) = sqrt(k) * (oow0 + dwx*x + dwy*y)
				// y = f/g*x =>
				// sqrt(a*x^2+2*b*x*f/g*x+c*f/g*x*f/g*x) = sqrt(k) * (oow0 + dwx*x + dwy*f/g*x)
				// sqrt(a*x^2*g^2+2*b*x^2*f*g+c*f^2*x^2) = sqrt(k) * (oow0*g + dwx*x*g + dwy*f*x)
				// sqrt(a*g^2+2*b*f*g+c*f^2) * x - sqrt(k)*dwx*g*x - sqrt(k)*dwy*f*x = sqrt(k)*oow0*g
				// (sqrt(a*g^2+2*b*f*g+c*f^2) - sqrt(k)*dwx*g - sqrt(k)*dwy*f) * x = sqrt(k)*oow0*g
				// x = sqrt(k)*oow0*g / (sqrt(a*g^2+2*b*f*g+c*f^2) - sqrt(k)*dwx*g - sqrt(k)*dwy*f)
				// y = sqrt(k)*oow0*f / (sqrt(a*g^2+2*b*f*g+c*f^2) - sqrt(k)*dwx*g - sqrt(k)*dwy*f)

				f = sin(a);
				g = cos(a);
				t = sk * w0 / (sqrt(J[0][0] * g * g + 2.0f * J[0][1] * f * g + J[1][1] * f * f) - sk * (g * Mi[0][2] + f * Mi[1][2]));
				glVertex3f(((midx + g * t)*vpScale.x)*figsx, ((midy + f * t)*vpScale.y)*figsy, -1);

				a += aa;
			}
			glEnd();
		}
	}
}
