gistfile1.txt

/*
Need cost function to divide by z, so that the errors we minimize are _projected_ errors.
Without it, 1px error might be 1mm if the ball is close up, or 100mm if the ball is further away. 
Thus to minimize error the optimizer would move the ball closer. By optimizing projected error there's no such bias.
*/
float CostFunction(vec3 dir, float r, vec3 pos)
{
	// Compute distance to the sphere surface
	vec3 posOnSphere = dir*dot(dir, pos);
	vec3 dirToRay = posOnSphere - pos;
	float distToSphere = length(dirToRay) - r;
	return distToSphere / posOnSphere.z;
}

vec3 Gradient(vec3 dir, float r, vec3 pos)
{	
	mat3 id(1);
	mat3 m = id - outerProduct(dir,dir);
	
	float dirDotPos = dot(dir, pos);
	vec3 posOnSphere = dir*dirDotPos;
	float distToSphereCenter = length(posOnSphere - pos);
	vec3 d_lengthDist = pos * (1.0f/distToSphereCenter); // length'(dist(pos))
	
	vec3 d_distToSphere = m * d_lengthDist; // Gradient for distance to sphere surface
	
	// oneOverSphereposz = 1/dir.z * dot(dir, pos)
	// oneOverSphereposz' = (dx/dz*dot(dir,pos), dy/dz*dot(dir,pos), 1/dot(dir,pos)^2)
	// Gradient for distance to sphere over the z value of the closet point using product rule:
	// (distToSphere/posOnSphere.z)' = distToSphere' / posOnSphere.z + distToSphere * (1/posOnSphere.z)'

	float distToSphere = distToSphereCenter - r;	
	vec3 d_oneOverClosestPosZ = vec3(dir.x / (dir.z*dirDotPos), dir.y / (dir.z*dirDotPos), 1 / (dirDotPos*dirDotPos));
	vec3 grad = d_distToSphere / posOnSphere.z + distToSphere * d_oneOverClosestPosZ;
				
	return grad;
}