Triangle class

using System;
using Sharp3D.Math.Core;

namespace PathTracer
{
	public class Triangle:Intersectable
	{
		Vector3F pointA, pointB, pointC;
		Vector3F normA, normB, normC;
		float xmin, ymin, zmin, xmax, ymax, zmax;
		Plane surfacePlane;
		Material surfaceMaterial;
		bool isLightSource;

		public Triangle (Vector3F pointA, Vector3F pointB, Vector3F pointC)
		{
			this.pointA = pointA;
			this.pointB = pointB;
			this.pointC = pointC;
			Vector3F normal = getNormal ();
			surfacePlane = new Plane (normal, pointA);
			computeBoundaries ();
		}

		public Triangle (Vector3F pointA, Vector3F pointB, Vector3F pointC, Material mat)
		{
			this.pointA = pointA;
			this.pointB = pointB;
			this.pointC = pointC;
			surfaceMaterial = mat;
			Vector3F normal = getNormal ();
			surfacePlane = new Plane (normal, pointA);
		}

		public Triangle (Vector3F a, Vector3F b, Vector3F c, Vector3F normA,
			Vector3F normB, Vector3F normC)
		{
			this.pointA = a;
			this.pointB = b;
			this.pointC = c;
			this.normA = normA;
			this.normB = normB;
			this.normC = normC;
			normA.Normalize ();
			normB.Normalize ();
			normC.Normalize ();
			Vector3F normal = getNormal ();
			surfacePlane = new Plane (normal, a);
			computeBoundaries ();
		}

		public Triangle (Vector3F a, Vector3F b, Vector3F c, Vector3F normA,
			Vector3F normB, Vector3F normC, Material mat)
		{
			this.pointA = a;
			this.pointB = b;
			this.pointC = c;
			this.normA = normA;
			this.normB = normB;
			this.normC = normC;
			normA.Normalize ();
			normB.Normalize ();
			normC.Normalize ();
			Vector3F normal = getNormal ();
			surfacePlane = new Plane (normal, a);
			computeBoundaries ();
			surfaceMaterial = mat;
		}

		private Vector3F getNormal ()
		{
			Vector3F ab = Vector3F.Subtract (pointB, pointA);
			ab.Normalize ();
			Vector3F ac = Vector3F.Subtract (pointC, pointA);
			ac.Normalize ();
			Vector3F normal = Vector3F.CrossProduct (ab, ac);
			return normal;
		}
		
		public Vector3F getNormal(float beta, float gamma) {
		// if the vertex normals have not been set, return the triangle plane's
		// normal
		Vector3F ret;
		if (normA == null) {
			return surfacePlane.getNormal();
		} else { // we return the linear interpolation based on the beta and
			// gamma factors of the normals.
			Vector3F Bfactor = new Vector3F(normB);
			Bfactor = Vector3F.Multiply(Bfactor,beta);
			Vector3F Cfactor = new Vector3F(normC);
			Cfactor = Vector3F.Multiply(Cfactor,gamma);
			Vector3F Afactor = new Vector3F(normA);
			Afactor = Vector3F.Multiply(Afactor,1 - beta - gamma);
			ret = new Vector3F(Afactor);
			ret = Vector3F.Add(ret,Bfactor);
			ret = Vector3F.Add(ret,Cfactor);
			ret.Normalize();
		}
		return ret;
	}

		public void computeBoundaries ()
		{
			xmin = pointA.X;
			ymin = pointA.Y;
			zmin = pointA.Z;
			xmax = pointA.X;
			ymax = pointA.Y;
			zmax = pointA.Z;
			if (pointB.X < xmin) {
				xmin = pointB.X;
			} else if (pointC.X < xmin) {
				xmin = pointC.X;
			}
			if (pointB.Y < ymin) {
				ymin = pointB.Y;
			} else if (pointC.Y < ymin) {
				ymin = pointC.Y;
			}
			if (pointB.Z < zmin) {
				zmin = pointB.Z;
			} else if (pointC.Z < zmin) {
				zmin = pointC.Z;
			}
			if (pointB.X > xmax) {
				xmax = pointB.X;
			} else if (pointC.X > xmax) {
				xmax = pointC.X;
			}
			if (pointB.Y > ymax) {
				ymax = pointB.Y;
			} else if (pointC.Y > ymax) {
				ymax = pointC.Y;
			}
			if (pointB.Z > zmax) {
				zmax = pointB.Z;
			} else if (pointC.Z > zmax) {
				zmax = pointC.Z;
			}
		}

		#region Intersectable implementation
		
		Hitrecord Intersectable.Intersect (Ray r)
		{
			Matrix3F A, Ax;
			float beta, gamma, t;
			Hitrecord hit;
			Vector3F hitPoint;
			Vector3F sol = new Vector3F (pointA.X - r.getOrigin ().X, pointA.Y
				- r.getOrigin ().Y, pointA.Z - r.getOrigin ().Z);

			A = new Matrix3F ();
			A.M11 = pointA.X - pointB.X; A.M12 = pointA.X - pointC.X; A.M13 = r.getDir().X;
			A.M21 = pointA.Y - pointB.Y; A.M22 = pointA.Y - pointC.Y; A.M23 = r.getDir().Z;
			A.M31 = pointA.Z - pointB.Z; A.M32 = pointA.Z - pointC.Z; A.M33 = r.getDir().Z;
			
			if (A.Determinant() < -Material.EPSILON
				|| A.Determinant () > Material.EPSILON) { // solving using
				// Cramer's rule
				float Adet = A.Determinant ();
				Ax = new Matrix3F (A);
				Ax.M11  = sol.X; Ax.M21 = sol.Y; Ax.M31 = sol.Z;
				beta = Ax.Determinant () / Adet;
				Ax = new Matrix3F (A);
				Ax.M21 = sol.X; Ax.M22=sol.Y; Ax.M23 = sol.Z;
				gamma = Ax.Determinant () / Adet;
				Ax = new Matrix3F (A);
				Ax.M31 = sol.X; Ax.M32 = sol.Y; Ax.M33 = sol.Z;
				t = Ax.Determinant () / Adet; // if beta + gamma are
				// between 0 and 1, we're in
				// the triangle.
				if (beta + gamma < 1 + Material.EPSILON && beta > -Material.EPSILON
					&& gamma > -Material.EPSILON) {
					hitPoint = new Vector3F (r.getDir ());
					hitPoint = Vector3F.Multiply (hitPoint, t);
					hitPoint = Vector3F.Add (hitPoint, r.getOrigin ());
					hit = new Hitrecord (hitPoint, t, getNormal (beta, gamma),
						((Intersectable)this).getMaterial (), this);
					return hit;
				}
			}
			return null;
		}

		public Vector3F getA ()
		{
			return pointA;
		}

		public Vector3F getB ()
		{
			return pointB;
		}

		public Vector3F getC ()
		{
			return pointC;
		}

		public Vector3F getNormA ()
		{
			return normA;
		}

		public Vector3F getNormB ()
		{
			return normB;
		}

		public Vector3F getNormC ()
		{
			return normC;
		}

		public float getXmax ()
		{
			return xmax;
		}

		public float getXmin ()
		{
			return xmin;
		}

		public float getYmax ()
		{
			return ymax;
		}

		public float getYmin ()
		{
			return ymin;
		}

		public float getZmax ()
		{
			return zmax;
		}

		public float getZmin ()
		{
			return zmin;
		}

		public bool isAggregate ()
		{
			return false;
		}

		Material Intersectable.getMaterial ()
		{
			return surfaceMaterial;
		}

		bool Intersectable.isLightSource ()
		{
			return isLightSource;
		}

		#endregion
	}
}