Real EdgeArea(const Vector &p, const Vector &q, const Vector &n)
{ 
	// Numerically stable edge-area calculator. (Including degenerate cases.)
	// Andrew J. Willmott, 1992

	const Real epsilon = 1e-6;
	Real sinFactor, qdotp, projArea;

	qdotp =  dot(p, q);
	sinFactor = sqrlen(p) * sqrlen(q) - sqr(qdotp);

	projArea = -dot(cross(p, q), n);

	// If SinFactor > 0 we may apply the formula without problem 
	
	if (sinFactor > 0)				
	{
		sinFactor = sqrt(sinFactor);
		return(projArea * (kHalfPi - atan(qdotp / sinFactor)) / (2 * kPi * sinFactor));
	}								
	
	// If not, we have three remaining cases... 

	else if (qdotp > epsilon) 		// CASE 1: p and q point in the same direction }
		return(0);
	else if (qdotp < -epsilon) 		// CASE 2: Viewpoint lies on the edge }
		return(0.5);
	else
		return(0.125);				// CASE 3: Viewpoint lies at either end of the edge }
}