{ RCSid $Id: fovsample.cal,v 1.2 2019/04/05 17:17:53 greg Exp $ }
{
	Foveated sampling based on central direction vector,
	magnification factor at center, and central field of view.

	Greg Ward	April 2019

	Preset constants:
		Cx, Cy, Cz :	central view direction
		M :		central FOV magnification factor
		D :		field of view (full angle in degrees)

	Inputs:
		iDx, iDy, iDz =	uniformly generated sampling vector
	
	Outputs:
		oDx, oDy, oDz =	foveated sample direction (normalized)
	
	Other variables computed:
		iTheta =	original angle to center (radians)
		oTheta =	new angle to center (radians)
}

r : PI/360 * D;		{ half-angle of foveal region in radians }
rp : M * r;		{ half-angle of mapped boundary }

len(x,y,z) : sqrt(x*x + y*y + z*z);
			{ normalize central vector }
Clen : len(Cx,Cy,Cz);
nCx : Cx/Clen;
nCy : Cy/Clen;
nCz : Cz/Clen;
			{ normalize input direction }
iDnf = 1/len(iDx,iDy,iDz);
niDx = iDx*iDnf;
niDy = iDy*iDnf;
niDz = iDz*iDnf;
			{ check for degenerate case (center of fovea) }
iDot = nCx*niDx + nCy*niDy + nCz*niDz;
degen = iDot - cos(.05*PI/180);

iTheta = acos(iDot);

gacc : 2;		{ gamma acceleration needed for C1 continuity }

gamma = 1 + ((iTheta-rp)/(PI-rp))^gacc * (log(M)/(log(PI/(PI-r)) - log(M)));

oTheta = if(rp-iTheta, iTheta/M,
		r + (PI-r)*((iTheta/(M*(PI-r)))^gamma - r/(PI-r)));

			{ normalized "up" vector for rotation }
Unf = 1/sqrt(1 - iDot*iDot);
nUx = (nCy*niDz - nCz*niDy)*Unf;
nUy = (nCz*niDx - nCx*niDz)*Unf;
nUz = (nCx*niDy - nCy*niDx)*Unf;
			{ rotate central vector by oTheta to get new output }
rcos = cos(oTheta);
rsin = sqrt(1 - rcos*rcos);
vPx = nUy*nCz - nUz*nCy;
vPy = nUz*nCx - nUx*nCz;
vPz = nUx*nCy - nUy*nCx;
			{ foveated sample direction }
rDx = nCx*rcos + vPx*rsin;
rDy = nCy*rcos + vPy*rsin;
rDz = nCz*rcos + vPz*rsin;
			{ substitute approximation in degenerate case }
oDx = if(degen, nCx+(niDx-nCx)/M, rDx);
oDy = if(degen, nCy+(niDy-nCy)/M, rDy);
oDz = if(degen, nCz+(niDz-nCz)/M, rDz);