Talk:Ray-sphere intersection

From DmWiki

It took me awhile to convince myself the optimized code was correct, here are the derivation steps to save others my confusion:

// Un-optimized
float intersectRaySphere(const Ray &ray, const Sphere &sphere) {
	// Assumes ray.d is normalized
	// Uses quadratic formula, A = 1 because ray.d is normalized.
	// For derivation, see:
	// http://www.siggraph.org/education/materials/HyperGraph/raytrace/rtinter1.htm
	Vec dst = ray.o - sphere.o;
	float B = 2.0f * dot(dst, ray.d);
	float C = dot(dst, dst) - sphere.r2;
	float D = B*B - 4.0f * C; // Discriminant.
	if (D < 0)
		return std::numeric_limits<float>::infinity();
	float sqrtD = sqrt(D);
	float t0 = (-B - sqrtD) * 0.5f; // Will always be the smaller, but may be negative
//	float t1 = (-B + sqrtD) * 0.5f; // The second intersection, if you care
	return t0;
}

We pull out the factor of 2 in B, which makes D = 2*B*2*B - 4.0f * C, which we then pull out the factor of 4, since that won't change whether it is negative or not for the test. After the test, we factor that 4 out of the sqrt (the factor becomes 2), which gives t0 = (-2.0f * B - 2.0f * sqrtD) * 0.5f, and the factors cancel out. Very clever optimization, but a derivation would have been nice.

// Optimized
float intersectRaySphere(const Ray &ray, const Sphere &sphere) {
	// Assumes ray.d is normalized
	// Uses quadratic formula, A = 1 because ray.d is normalized.
	// For derivation, see:
	// http://www.siggraph.org/education/materials/HyperGraph/raytrace/rtinter1.htm
	Vec dst = ray.o - sphere.o;
	float B = dot(dst, ray.d); // B really equals 2.0f this value
	float C = dot(dst, dst) - sphere.r2;
	float D = B*B - C;  // Discriminant.  D really equals 4.0f times this value
	if (D < 0)
		return std::numeric_limits<float>::infinity();
	float sqrtD = sqrt(D); // sqrtD really equals 2.0f times this value
//	float t0 = (-2.0f * B - 2.0f * sqrtD) * 0.5f;  // The factors cancel out
	float t0 = (-B - sqrtD); // Will always be the smaller, but may be negative
//	float t1 = (-B + sqrtD); // The second intersection, if you care
	return t0;
}