in libraries/hvvr/raycaster/samples.cpp [44:142]
DynamicArray<DirectionalBeam> getEyeSpacePolarFoveatedSamples(size_t& samplesPerRing,
EccentricityMap& emap,
float maxEyeTrackingUncertaintyDegrees,
float maxFOVDegrees,
float marSlope,
float fovealMARDegrees) {
const float m = marSlope;
const float w_0 = fovealMARDegrees;
const float switchPoint1 = maxEyeTrackingUncertaintyDegrees / w_0;
const float S = maxEyeTrackingUncertaintyDegrees;
emap = EccentricityMap(marSlope, fovealMARDegrees, maxEyeTrackingUncertaintyDegrees);
samplesPerRing = 0;
size_t ringCount = 0;
size_t irregularGridSampleCount = 0;
{ // Calculate number of samples so we can allocate before generation
float E = 0.0f;
while (E <= maxFOVDegrees * RadiansPerDegree) {
// Angular distance (in degrees) between samples on this annulus
float w = emap.apply(ringCount + 1.0f) - E;
float ringRadius = sinf(E);
float angularDistanceAroundRing = 2.0f * Pi * ringRadius;
size_t samplesOnAnnulus = (size_t)(std::ceil(angularDistanceAroundRing / w));
printf("New samplesOnAnnulus: %zu = ceil(%f/%f)\n", samplesOnAnnulus, angularDistanceAroundRing, w);
samplesPerRing = std::max(samplesPerRing, samplesOnAnnulus);
irregularGridSampleCount += samplesOnAnnulus;
++ringCount;
E = emap.apply((float)ringCount);
}
}
printf("(%zu*%zu=%d)/%d %f times the minimal sample count\n", samplesPerRing, ringCount,
(int)(samplesPerRing * ringCount), (int)irregularGridSampleCount,
(samplesPerRing * ringCount) / (float)irregularGridSampleCount);
DynamicArray<DirectionalBeam> samples(ringCount * samplesPerRing);
/**
A note on differentials:
The zenith differential is found by taking the nearest point on the next concentric ring outwards, and
projecting it onto
the tangent plane of the current direction (where the direction is on the unit sphere)
The vector from the intersection point of the current direction and the tangent plane to the intersection point
of the
"nearest direction on the next ring" and the tangent plane is the zenith differential
Ray-plane intersection:
Ray: P = P_0 + tV
Plane: dot(P, N) + d = 0
For our purposes, d = 1, P_0 is zero, so we actually have
Ray: P = tV
Plane: dot(P, N) = -1
We'll reverse the convention, so that the normal of the plane faces away from the origin
Ray: P = tV
Plane: dot(P, N) = 1
dot(tV,N) = 1
t * dot(V,N) = 1
t = 1 / dot(V,N)
P = t*V
P = V/dot(V,N) <- all we need
Then the zenith differential is just P-N
The azimuthal differential is taken by finding a perpendicular vector to the zenith differential,
and then scaling it to the distance along the ring to the next sample
*/
// Generate concentric circles of samples with spacing equal or less than MAR of eye at the eccentricity
float E = 0.0f;
for (size_t r = 0; r < ringCount; ++r) {
// Angular distance (in degrees) between samples on this annulus
float E_next = emap.apply(r + 1.0f);
float ringRadius = sinf(E);
float angularDistanceAroundRing = 2.0f * Pi * ringRadius;
vector3 baseVector = normalize(quaternion::fromAxisAngle(vector3(0, 1, 0), E) * vector3(0, 0, -1));
vector3 zenithDiffBaseVector =
normalize(quaternion::fromAxisAngle(vector3(0, 1, 0), E_next) * vector3(0, 0, -1));
for (size_t i = 0; i < samplesPerRing; ++i) {
float rotationRadians = (i + 0.5f) / float(samplesPerRing) * 2.0f * Pi;
vector3 p = normalize(quaternion::fromAxisAngle(vector3(0, 0, -1), rotationRadians) * baseVector);
vector3 zenithDiffDirection =
normalize(quaternion::fromAxisAngle(vector3(0, 0, -1), rotationRadians) * zenithDiffBaseVector);
// Project zenith direction onto tangent plane of p
// P = V/dot(V,N)
vector3 zenithDiffDirectionOnTangentPlane = zenithDiffDirection * (1.0f / dot(zenithDiffDirection, p));
vector3 zenithDiff = zenithDiffDirectionOnTangentPlane - p;
vector3 azimuthalDiffUnit = cross(p, normalize(zenithDiff));
vector3 azimuthalDiff = azimuthalDiffUnit * angularDistanceAroundRing / (float)samplesPerRing;
size_t idx = samplesPerRing * r + i;
samples[idx].centerRay = p;
samples[idx].du = zenithDiff;
samples[idx].dv = azimuthalDiff;
}
E = E_next;
}
return samples;
}