{ RCSid $Id: fovsample.cal,v 1.5 2019/04/09 00:52:55 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:
		forward =	1 if forward conversion, -1 if reverse
		iDx, iDy, iDz =	input sampling vector
	
	Outputs:
		oDx, oDy, oDz =	output sample direction (normalized)
	
	Other variables computed:
		iTheta =	input angle to center (radians)
		oTheta =	output 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);
sq(x) : x*x;
			{ 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 or opposite) }
iDot = nCx*niDx + nCy*niDy + nCz*niDz;
degen0 = iDot - cos(.05*PI/180);
degen1 = cos(179.95*PI/180) - iDot;

iTheta = acos(iDot);
			{ quadratic function coefficients mate to linear ramp }
qa : (PI - PI/M)/sq(PI - rp);
qb : 1/M - 2*PI*r*(M - 1)/sq(PI - rp);
qc : (PI - PI/M)/sq(PI/rp - 1);

oTheta = if (forward,
		if(rp-iTheta, iTheta/M, qa*iTheta*iTheta + qb*iTheta + qc),
		if(r-iTheta, iTheta*M, (sqrt(qb*qb - 4*qa*(qc - iTheta)) - qb)/(2*qa))
	);
			{ 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);
			{ foveated sample direction }
rDx = nCx*rcos + (nUy*nCz - nUz*nCy)*rsin;
rDy = nCy*rcos + (nUz*nCx - nUx*nCz)*rsin;
rDz = nCz*rcos + (nUx*nCy - nUy*nCx)*rsin;
			{ substitute approximation in degenerate case }
oDx = if(degen0, nCx+(niDx-nCx)*if(forward,1/M,M), if(degen1, niDx, rDx));
oDy = if(degen0, nCy+(niDy-nCy)*if(forward,1/M,M), if(degen1, niDy, rDy));
oDz = if(degen0, nCz+(niDz-nCz)*if(forward,1/M,M), if(degen1, niDz, rDz));