Skip to content

Instantly share code, notes, and snippets.

@castano

castano/hemicube.cpp

Created June 20, 2014 09:46
Show Gist options
  • Save castano/5f7c8847f82e45cb56d0 to your computer and use it in GitHub Desktop.
Save castano/5f7c8847f82e45cb56d0 to your computer and use it in GitHub Desktop.
Download ZIP
Hemicube Integrator
Raw
hemicube.cpp
#include "hemicube.h"
#define PACK_HEMICUBES 1
static void get_hemicube_face_normal(int index, Vector3 *forward, Vector3 *left, Vector3 *up) {
// Unwrapped hemicube with positive-Z in the middle.
switch (index) {
case 0: *forward = Vector3(+1, 0, 0); *left = Vector3( 0, 1, 0); break;
case 1: *forward = Vector3(-1, 0, 0); *left = Vector3( 0, 1, 0); break;
case 2: *forward = Vector3( 0, +1, 0); *left = Vector3( 0, 0, -1); break;
case 3: *forward = Vector3( 0, -1, 0); *left = Vector3( 0, 0, 1); break;
case 4: *forward = Vector3( 0, 0, +1); *left = Vector3( 0, 1, 0); break;
}
*up = cross(*forward, *left);
}
static int face_coordinates(int x, int y, int cube_size, Vector2 * f) {
#if PACK_HEMICUBES == 0
x += cube_size/2;
y += cube_size/2;
// Compute face and translate x, y
int face = -1;
float fx, fy;
if (x < cube_size) {
if (y >= cube_size && y < 2*cube_size) {
face = 2; // positive y
fx = x;
fy = y - cube_size;
}
}
else if (x >= 2*cube_size) {
if (y >= cube_size && y < 2*cube_size) {
face = 3; // negative y
fx = x - 2*cube_size;
fy = y - cube_size;
}
}
else {
if (y < cube_size) {
face = 1; // positive x
fx = x - cube_size;
fy = y;
}
else if (y >= 2*cube_size) {
face = 0; // negative x
fx = x - cube_size;
fy = y - 2*cube_size;
}
else {
face = 4; // positive z
fx = x - cube_size;
fy = y - cube_size;
}
}
if (face == -1) {
// Outside.
return -1;
}
#else
int face = -1;
float fx, fy;
if (x < cube_size) {
face = 4; // positive z
fx = x;
fy = y;
}
else if (x < 2*cube_size) {
if (y < cube_size/2) {
face = 1; // X-
fx = x - cube_size;
fy = y + cube_size/2;
}
else {
face = 0; // X+
fx = x - cube_size;
fy = y - cube_size/2;
}
}
else {
if (x < 2*cube_size + cube_size/2) {
face = 2; // Y+
fx = x - 2 * cube_size + cube_size/2;
fy = y;
}
else {
face = 3; // Y-
fx = x - 2 * cube_size - cube_size/2;
fy = y;
}
}
assert (face != -1);
#endif
assert(fx >= 0 && fx < cube_size);
assert(fy >= 0 && fy < cube_size);
// Transform to [-1, 1] range, offset by 0.5 to point to texel center.
f->x = 2 * (fx + 0.5f) / cube_size - 1;
f->y = 2 * (fy + 0.5f) / cube_size - 1;
return face;
}
static void atlas_coordinate(int x, int y, int face, int cube_size, int * xx, int * yy) {
// Adjust face origin. (We do this first to simplify the wrap around rules)
if (face == 1) y -= cube_size/2;
if (face == 2) x -= cube_size/2;
// Handle wrap around:
if (face == 4) {
if (x == -1) {
face = 2;
x = cube_size/2 - 1;
}
else if (x == cube_size) {
face = 3;
x = 0;
}
else if (y == -1) {
face = 1;
y = cube_size/2 - 1;
}
else if (y == cube_size) {
face = 0;
y = 0;
}
}
else if (face == 0) {
if (x == -1) {
face = 2;
x = cube_size/2 - y - 1;
y = cube_size - 1;
}
else if (x == cube_size) {
face = 3;
x = y;
y = cube_size - 1;
}
}
else if (face == 1) {
if (x == -1) {
face = 2;
x = y;
y = 0;
}
else if (x == cube_size) {
face = 3;
x = cube_size/2 - y - 1;
y = 0;
}
}
// The remaining faces are not allowed to wrap around.
#if PACK_HEMICUBES == 0
if (face == 0) {
assert (x >= 0 && x < cube_size);
assert (y >= 0 && y < cube_size/2);
*xx = x + cube_size/2;
*yy = y + 3*cube_size/2;
}
else if (face == 1) {
assert (x >= 0 && x < cube_size);
assert (y >= 0 && y < cube_size/2);
*xx = x + cube_size/2;
*yy = y + 0;
}
else if (face == 2) {
assert (x >= 0 && x < cube_size/2);
assert (y >= 0 && y < cube_size);
*xx = x + 0;
*yy = y + cube_size/2;
}
else if (face == 3) {
assert (x >= 0 && x < cube_size/2);
assert (y >= 0 && y < cube_size);
*xx = x + 3*cube_size/2;
*yy = y + cube_size/2;
}
else if (face == 4) {
assert (x >= -1 && x <= cube_size);
assert (y >= -1 && y <= cube_size);
*xx = x + cube_size/2;
*yy = y + cube_size/2;
}
assert (*xx >= 0 && x < 2*cube_size);
assert (*yy >= 0 && y < 2*cube_size);
#else
if (face == 0) {
assert (x >= 0 && x < cube_size);
assert (y >= 0 && y < cube_size/2);
*xx = x + cube_size;
*yy = y + cube_size/2;
}
else if (face == 1) {
assert (x >= 0 && x < cube_size);
assert (y >= 0 && y < cube_size/2);
*xx = x + cube_size;
*yy = y + 0;
}
else if (face == 2) {
assert (x >= 0 && x < cube_size/2);
assert (y >= 0 && y < cube_size);
*xx = x + 2*cube_size;
*yy = y + 0;
}
else if (face == 3) {
assert (x >= 0 && x < cube_size/2);
assert (y >= 0 && y < cube_size);
*xx = x + 2*cube_size + cube_size/2;
*yy = y + 0;
}
else if (face == 4) {
assert (x >= 0 && x < cube_size);
assert (y >= 0 && y < cube_size);
*xx = x;
*yy = y;
}
assert (*xx >= 0 && x < 3*cube_size);
assert (*yy >= 0 && y < cube_size);
#endif
}
// Compute direction vector to the hemicube texel.
static Vector3 face_vector(const Vector2 & f, int face) {
Vector3 forward, left, up;
get_hemicube_face_normal(face, (Vector3*)&forward, (Vector3*)&left, (Vector3*)&up);
Vector3 d = normalize(forward - left * f.x - up * f.y);
assert(d.z >= 0);
return d;
}
// Solid angle of an axis aligned quad from (0,0,1) to (x,y,1)
// See: http://www.fizzmoll11.com/thesis/ for a derivation of this formula.
static float area_element(float x, float y) {
return atan2(x*y, sqrtf(x*x + y*y + 1));
}
// Solid angle of a hemicube texel.
static float solid_angle_term(const Vector2 & f, int cube_size) {
#if 1 // Exact solid angle:
float x0 = f.x - (1.0f / cube_size);
float y0 = f.y - (1.0f / cube_size);
float x1 = f.x + (1.0f / cube_size);
float y1 = f.y + (1.0f / cube_size);
float solid_angle = area_element(x0, y0) - area_element(x0, y1) - area_element(x1, y0) + area_element(x1, y1);
assert(solid_angle > 0.0f);
return solid_angle;
#else // This formula is equivalent, but not as precise.
float pixel_area = square(2.0f / cube_size);
float dist_square = 1.0f + square(f.x) + square(f.y);
float cos_theta = 1.0f / sqrt(dist_square);
float cos_theta_d2 = cos_theta / dist_square; // Funny this is just 1/dist^3 or cos(tetha)^3
return pixel_area * cos_theta_d2;
#endif
}
static float clamped_cosine_term(const Vector3 & dir) {
return clamp(dir.z, 0.0f, 1.0f);
}
// This computes the gradients of the cosine weighted solid angle of the hemicube.
static Vector2 static_irradiance_gradients(int face, const Vector3 & dir) {
float dx, dy;
if (face == 0 || face == 1) { // YZ faces:
float yi = dir.y / fabs(dir.x);
float zi = dir.z / fabs(dir.x);
assert(zi >= 0.0f);
float l2 = (1 + yi*yi + zi*zi);
//dx = 4 * zi / (l2*l2*l2) - zi / (l2*l2);
dx = 3 / powf(l2, 5.0f/2.0f) - 1 / powf(l2, 3.0f/2.0f);
if(dir.x < 0) dx = -dx;
//dy = 4 * yi * zi / (l2*l2*l2);
dy = 3 * yi / powf(l2, 5.0f/2.0f);
// sanity checks:
assert(dx/fabs(dx) == dir.x/fabs(dir.x));
}
else if (face == 2 || face == 3) { // XZ faces:
float xi = dir.x / fabs(dir.y);
float zi = dir.z / fabs(dir.y);
assert(zi >= 0.0f);
float l2 = (xi*xi + 1 + zi*zi);
//dx = 4 * xi * zi / (l2*l2*l2);
dx = 3 * xi / powf(l2, 5.0f/2.0f);
//dy = 4 * zi / (l2*l2*l2) - zi / (l2*l2);
dy = 3 / powf(l2, 5.0f/2.0f) - 1 / powf(l2, 3.0f/2.0f);
if(dir.y < 0) dy = -dy;
// sanity checks:
assert(dy/fabs(dy) == dir.y/fabs(dir.y));
}
else /*if (face == 4)*/ { // XY face:
assert(face == 4);
float xi = dir.x / dir.z;
float yi = dir.y / dir.z;
float l2 = (xi*xi + yi*yi + 1);
//dx = 4 * xi / (l2*l2*l2);
//dy = 4 * yi / (l2*l2*l2);
dx = 3 * xi / powf(l2, 5.0f/2.0f);
dy = 3 * yi / powf(l2, 5.0f/2.0f);
// sanity checks:
assert(xi/fabs(xi) == dx/fabs(dx));
assert(yi/fabs(yi) == dy/fabs(dy));
}
dx *= dir.z;
dy *= dir.z;
return Vector2(dx, dy);
}
Hemicube_Data::Hemicube_Data() {
cube_size = 0;
}
void Hemicube_Data::init(const int cube_size) {
if (cube_size == this->cube_size) return;
this->cube_size = cube_size;
#if PACK_HEMICUBES == 0
width = cube_size * 2;
height = cube_size * 2;
#else
width = cube_size * 3;
height = cube_size;
#endif
#if PACK_HEMICUBES == 0
viewport_offset_x[0] = cube_size/2;
viewport_offset_y[0] = 3*cube_size/2;
viewport_offset_x[1] = cube_size/2;
viewport_offset_y[1] = 0;
viewport_offset_x[2] = 0;
viewport_offset_y[2] = cube_size/2;
viewport_offset_x[3] = 3*cube_size/2;
viewport_offset_y[3] = cube_size/2;
viewport_offset_x[4] = cube_size/2;
viewport_offset_y[4] = cube_size/2;
#else
viewport_offset_x[0] = cube_size;
viewport_offset_y[0] = cube_size/2;
viewport_offset_x[1] = cube_size;
viewport_offset_y[1] = 0;
viewport_offset_x[2] = 2*cube_size;
viewport_offset_y[2] = 0;
viewport_offset_x[3] = 5*cube_size/2;
viewport_offset_y[3] = 0;
viewport_offset_x[4] = 0;
viewport_offset_y[4] = 0;
#endif
viewport_size_x[0] = cube_size;
viewport_size_y[0] = cube_size/2;
viewport_size_x[1] = cube_size;
viewport_size_y[1] = cube_size/2;
viewport_size_x[2] = cube_size/2;
viewport_size_y[2] = cube_size;
viewport_size_x[3] = cube_size/2;
viewport_size_y[3] = cube_size;
viewport_size_x[4] = cube_size;
viewport_size_y[4] = cube_size;
texels.resize(width * height);
for (int y = 0; y < height; y++) {
for (int x = 0; x < width; x++) {
Hemicube_Texel & texel = texels[y * width + x];
texel.face = face_coordinates(x, y, cube_size, &texel.f);
if (texel.face != -1) {
texel.dir = face_vector(texel.f, texel.face);
texel.solid_angle = solid_angle_term(texel.f, cube_size);
texel.clamped_cosine = clamped_cosine_term(texel.dir);
//texel.static_gradient = static_irradiance_gradients(texel.face, texel.dir);
}
}
}
// Computing the marginal change in irradiance for this cell reduces to computing the marginal
// change in the two highlighted cell walls with respect to translation
edges.reset();
edges.reserve(cube_size * cube_size * 6);
// Y-walls of top face:
for (int y = 0; y < cube_size; y++) {
for (int x = 0; x <= cube_size; x++) {
Hemicube_Edge edge;
atlas_coordinate(x-1, y, 4, cube_size, &edge.x0, &edge.y0);
atlas_coordinate(x, y, 4, cube_size, &edge.x1, &edge.y1);
// Endpoints of this wall:
float fx = 2 * float(x) / cube_size - 1;
float fy0 = 2 * float(y+0) / cube_size - 1;
float fy1 = 2 * float(y+1) / cube_size - 1;
Vector3 d0 = face_vector(Vector2(fx, fy0), 4);
Vector3 d1 = face_vector(Vector2(fx, fy1), 4);
// Length of the arc described by this wall:
float wall_length = acos(dot(d0, d1));
// Wall normal projected to the XY plane:
Vector3 wall_normal = cross(d0,d1);
Vector2 projected_wall_normal = normalize(wall_normal.xy());
assert(equal(projected_wall_normal, Vector2(0, 1)));
// Length of the cell wall multiplied by the rate of motion of the wall with respect to motion in the normal direction
float fy = 2 * float(y+0.5f) / cube_size - 1;
Vector3 dir = face_vector(Vector2(fx, fy), 4);
float cos_theta = dir.z;
float motion_rate = -cos_theta; // d(theta)/dx = -cos(theta) / depth
float marginal_change = wall_length * motion_rate;
edge.gradient = projected_wall_normal * marginal_change;
edges.append(edge);
}
}
// X-walls for top face:
for (int y = 0; y <= cube_size; y++) {
for (int x = 0; x < cube_size; x++) {
Hemicube_Edge edge;
atlas_coordinate(x, y, 4, cube_size, &edge.x0, &edge.y0);
atlas_coordinate(x, y-1, 4, cube_size, &edge.x1, &edge.y1);
// Endpoints of this wall:
float fx0 = 2 * float(x+0) / cube_size - 1;
float fx1 = 2 * float(x+1) / cube_size - 1;
float fy = 2 * float(y) / cube_size - 1;
Vector3 d0 = face_vector(Vector2(fx0, fy), 4);
Vector3 d1 = face_vector(Vector2(fx1, fy), 4);
// Length of the arc described by this wall:
float wall_length = acos(dot(d0, d1));
// Wall normal projected to the XY plane:
Vector3 wall_normal = cross(d0,d1);
Vector2 projected_wall_normal = normalize(wall_normal.xy());
assert(equal(projected_wall_normal, Vector2(1, 0)));
// Length of the cell wall multiplied by the rate of motion of the wall with respect to motion in the normal direction
float fx = 2 * float(x+0.5f) / cube_size - 1;
Vector3 dir = face_vector(Vector2(fx, fy), 4);
float cos_theta = dir.z;
float motion_rate = -cos_theta; // d(theta)/dx = -cos(theta) / depth
float marginal_change = wall_length * motion_rate;
edge.gradient = projected_wall_normal * marginal_change;
edges.append(edge);
}
}
// X-walls for X+ face.
for (int z = 1; z < cube_size/2; z++) {
for (int y = 0; y < cube_size; y++) {
Hemicube_Edge edge;
atlas_coordinate(y, z, 0, cube_size, &edge.x0, &edge.y0);
atlas_coordinate(y, z-1, 0, cube_size, &edge.x1, &edge.y1);
// Endpoints of this wall:
float fz = 2 * float(z) / cube_size - 1;
float fy0 = 2 * float(y+0) / cube_size - 1;
float fy1 = 2 * float(y+1) / cube_size - 1;
Vector3 d0 = face_vector(Vector2(fy0, fz), 0);
Vector3 d1 = face_vector(Vector2(fy1, fz), 0);
// Length of the arc described by this wall:
float wall_length = acos(dot(d0, d1));
// Wall normal projected to the XY plane:
Vector3 wall_normal = cross(d0,d1);
Vector2 projected_wall_normal = normalize(wall_normal.xy());
assert(equal(projected_wall_normal, Vector2(1, 0)));
// Length of the cell wall multiplied by the rate of motion of the wall with respect to motion in the normal direction
float fy = 2 * float(y+0.5f) / cube_size - 1;
Vector3 dir = face_vector(Vector2(fy, fz), 0);
float cos_theta = dir.z;
float motion_rate = -cos_theta; // d(theta)/dx = -cos(theta) / depth
float marginal_change = wall_length * motion_rate;
edge.gradient = projected_wall_normal * marginal_change;
edges.append(edge);
}
}
// X-walls for X- face.
for (int z = cube_size/2+1; z < cube_size; z++) {
for (int y = 0; y < cube_size; y++) {
Hemicube_Edge edge;
atlas_coordinate(y, z, 1, cube_size, &edge.x0, &edge.y0);
atlas_coordinate(y, z-1, 1, cube_size, &edge.x1, &edge.y1);
// Endpoints of this wall:
float fz = 2 * float(z) / cube_size - 1;
float fy0 = 2 * float(y+0) / cube_size - 1;
float fy1 = 2 * float(y+1) / cube_size - 1;
Vector3 d0 = face_vector(Vector2(fy0, fz), 1);
Vector3 d1 = face_vector(Vector2(fy1, fz), 1);
// Length of the arc described by this wall:
float wall_length = acos(dot(d0, d1));
// Wall normal projected to the XY plane:
Vector3 wall_normal = cross(d0,d1);
Vector2 projected_wall_normal = normalize(wall_normal.xy());
assert(equal(projected_wall_normal, Vector2(1, 0)));
// Length of the cell wall multiplied by the rate of motion of the wall with respect to motion in the normal direction
float fy = 2 * float(y+0.5f) / cube_size - 1;
Vector3 dir = face_vector(Vector2(fy, fz), 1);
float cos_theta = dir.z;
float motion_rate = -cos_theta; // d(theta)/dx = -cos(theta) / depth
float marginal_change = wall_length * motion_rate;
edge.gradient = projected_wall_normal * marginal_change;
edges.append(edge);
}
}
// Y-walls for Y- face.
for (int x = 0; x < cube_size; x++) {
for (int z = 1; z < cube_size/2; z++) {
Hemicube_Edge edge;
atlas_coordinate(z-1, x, 3, cube_size, &edge.x0, &edge.y0);
atlas_coordinate(z, x, 3, cube_size, &edge.x1, &edge.y1);
// Endpoints of this wall:
float fz = 2 * float(z) / cube_size - 1;
float fx0 = 2 * float(x+0) / cube_size - 1;
float fx1 = 2 * float(x+1) / cube_size - 1;
Vector3 d0 = face_vector(Vector2(fz, fx0), 3);
Vector3 d1 = face_vector(Vector2(fz, fx1), 3);
// Length of the arc described by this wall:
float wall_length = acos(dot(d0, d1));
// Wall normal projected to the XY plane:
Vector3 wall_normal = cross(d0,d1);
Vector2 projected_wall_normal = normalize(wall_normal.xy());
assert(equal(projected_wall_normal, Vector2(0, 1)));
// Length of the cell wall multiplied by the rate of motion of the wall with respect to motion in the normal direction
float fx = 2 * float(x+0.5f) / cube_size - 1;
Vector3 dir = face_vector(Vector2(fz, fx), 3);
float cos_theta = dir.z;
float motion_rate = -cos_theta; // d(theta)/dx = -cos(theta) / depth
float marginal_change = wall_length * motion_rate;
edge.gradient = projected_wall_normal * marginal_change;
edges.append(edge);
}
}
// Y-walls for Y+ face.
for (int x = 0; x < cube_size; x++) {
for (int z = cube_size/2+1; z < cube_size; z++) {
Hemicube_Edge edge;
atlas_coordinate(z-1, x, 2, cube_size, &edge.x0, &edge.y0);
atlas_coordinate(z, x, 2, cube_size, &edge.x1, &edge.y1);
// Endpoints of this wall:
float fz = 2 * float(z) / cube_size - 1;
float fx0 = 2 * float(x+0) / cube_size - 1;
float fx1 = 2 * float(x+1) / cube_size - 1;
Vector3 d0 = face_vector(Vector2(fz, fx0), 2);
Vector3 d1 = face_vector(Vector2(fz, fx1), 2);
// Length of the arc described by this wall:
float wall_length = acos(dot(d0, d1));
// Wall normal projected to the XY plane:
Vector3 wall_normal = cross(d0,d1);
Vector2 projected_wall_normal = normalize(wall_normal.xy());
assert(equal(projected_wall_normal, Vector2(0, 1)));
// Length of the cell wall multiplied by the rate of motion of the wall with respect to motion in the normal direction
float fx = 2 * float(x+0.5f) / cube_size - 1;
Vector3 dir = face_vector(Vector2(fz, fx), 2);
float cos_theta = dir.z;
float motion_rate = -cos_theta; // d(theta)/dx = -cos(theta) / depth
float marginal_change = wall_length * motion_rate;
edge.gradient = projected_wall_normal * marginal_change;
edges.append(edge);
}
}
// Z-walls for X+ face.
for (int z = 0; z < cube_size/2; z++) {
for (int y = 0; y <= cube_size; y++) {
Hemicube_Edge edge;
atlas_coordinate(y-1, z, 0, cube_size, &edge.x0, &edge.y0);
atlas_coordinate(y, z, 0, cube_size, &edge.x1, &edge.y1);
// Endpoints of this wall:
float fz0 = 2 * float(z+0) / cube_size - 1;
float fz1 = 2 * float(z+1) / cube_size - 1;
float fy = 2 * float(y) / cube_size - 1;
Vector3 d0 = face_vector(Vector2(fy, fz0), 0);
Vector3 d1 = face_vector(Vector2(fy, fz1), 0);
// Length of the arc described by this wall:
float wall_length = acos(dot(d0, d1));
// Wall normal projected to the XY plane:
Vector3 wall_normal = cross(d0,d1);
Vector2 projected_wall_normal = normalize(wall_normal.xy());
// Length of the cell wall multiplied by the rate of motion of the wall with respect to motion in the normal direction
float fz = 2 * float(z+0.5f) / cube_size - 1;
Vector3 dir = face_vector(Vector2(fy, fz), 0);
float cos_theta = dir.z;
float sin_theta = sqrtf(1 - cos_theta*cos_theta);
float motion_rate = -1 / sin_theta; // d(phi)/dx = -1 / (sin(theta) * depth)
float marginal_change = wall_length * motion_rate;
edge.gradient = projected_wall_normal * marginal_change;
edges.append(edge);
}
}
// Z-walls for X- face.
for (int z = cube_size/2; z < cube_size; z++) {
for (int y = 0; y <= cube_size; y++) {
Hemicube_Edge edge;
atlas_coordinate(y-1, z, 1, cube_size, &edge.x0, &edge.y0);
atlas_coordinate(y, z, 1, cube_size, &edge.x1, &edge.y1);
// Endpoints of this wall:
float fz0 = 2 * float(z+0) / cube_size - 1;
float fz1 = 2 * float(z+1) / cube_size - 1;
float fy = 2 * float(y) / cube_size - 1;
Vector3 d0 = face_vector(Vector2(fy, fz0), 1);
Vector3 d1 = face_vector(Vector2(fy, fz1), 1);
// Length of the arc described by this wall:
float wall_length = acos(dot(d0, d1));
// Wall normal projected to the XY plane:
Vector3 wall_normal = cross(d0,d1);
Vector2 projected_wall_normal = normalize(wall_normal.xy());
// Length of the cell wall multiplied by the rate of motion of the wall with respect to motion in the normal direction
float fz = 2 * float(z+0.5f) / cube_size - 1;
Vector3 dir = face_vector(Vector2(fy, fz), 1);
float cos_theta = dir.z;
float sin_theta = sqrtf(1 - cos_theta*cos_theta);
float motion_rate = -1 / sin_theta; // d(phi)/dx = -1 / (sin(theta) * depth)
float marginal_change = wall_length * motion_rate;
edge.gradient = projected_wall_normal * marginal_change;
edges.append(edge);
}
}
// Z-walls for Y- face.
for (int x = 1; x < cube_size; x++) {
for (int z = 0; z < cube_size/2; z++) {
Hemicube_Edge edge;
atlas_coordinate(z, x, 3, cube_size, &edge.x0, &edge.y0);
atlas_coordinate(z, x-1, 3, cube_size, &edge.x1, &edge.y1);
// Endpoints of this wall:
float fz0 = 2 * float(z+0) / cube_size - 1;
float fz1 = 2 * float(z+1) / cube_size - 1;
float fx = 2 * float(x) / cube_size - 1;
Vector3 d0 = face_vector(Vector2(fz0, fx), 3);
Vector3 d1 = face_vector(Vector2(fz1, fx), 3);
// Length of the arc described by this wall:
float wall_length = acos(dot(d0, d1));
// Wall normal projected to the XY plane:
Vector3 wall_normal = cross(d0,d1);
Vector2 projected_wall_normal = normalize(wall_normal.xy());
// Length of the cell wall multiplied by the rate of motion of the wall with respect to motion in the normal direction
float fz = 2 * float(z+0.5f) / cube_size - 1;
Vector3 dir = face_vector(Vector2(fz, fx), 3);
float cos_theta = dir.z;
float sin_theta = sqrtf(1 - cos_theta*cos_theta);
float motion_rate = -1 / sin_theta; // d(phi)/dx = -1 / (sin(theta) * depth)
float marginal_change = wall_length * motion_rate;
edge.gradient = projected_wall_normal * marginal_change;
edges.append(edge);
}
}
// Z-walls for Y+ face.
for (int x = 1; x < cube_size; x++) {
for (int z = cube_size/2; z < cube_size; z++) {
Hemicube_Edge edge;
atlas_coordinate(z, x, 2, cube_size, &edge.x0, &edge.y0);
atlas_coordinate(z, x-1, 2, cube_size, &edge.x1, &edge.y1);
// Endpoints of this wall:
float fz0 = 2 * float(z+0) / cube_size - 1;
float fz1 = 2 * float(z+1) / cube_size - 1;
float fx = 2 * float(x) / cube_size - 1;
Vector3 d0 = face_vector(Vector2(fz0, fx), 2);
Vector3 d1 = face_vector(Vector2(fz1, fx), 2);
// Length of the arc described by this wall:
float wall_length = acos(dot(d0, d1));
// Wall normal projected to the XY plane:
Vector3 wall_normal = cross(d0,d1);
Vector2 projected_wall_normal = normalize(wall_normal.xy());
// Length of the cell wall multiplied by the rate of motion of the wall with respect to motion in the normal direction
float fz = 2 * float(z+0.5f) / cube_size - 1;
Vector3 dir = face_vector(Vector2(fz, fx), 2);
float cos_theta = dir.z;
float sin_theta = sqrtf(1 - cos_theta*cos_theta);
float motion_rate = -1 / sin_theta; // d(phi)/dx = -1 / (sin(theta) * depth)
float marginal_change = wall_length * motion_rate;
edge.gradient = projected_wall_normal * marginal_change;
edges.append(edge);
}
}
}
Raw
hemicube.h
struct Hemicube_Texel {
int face; // Face index.
Vector2 f; // Texel coordinates within the hemicube face in [-1, 1] range.
Vector3 dir; // Direction vector to center of hemicube texel.
float solid_angle; // Projected solid angle of hemimicube texel.
float clamped_cosine; // max(0, dir.z)
};
struct Hemicube_Edge {
int x0, y0, x1, y1; // Coordinates of adjacent texels.
Vector2 gradient; // Translation gradient contribution of this edge.
};
struct Hemicube_Data {
Hemicube_Data();
void init(int cube_size);
int cube_size;
int width;
int height;
int viewport_offset_x[5];
int viewport_offset_y[5];
int viewport_size_x[5];
int viewport_size_y[5];
Array <Hemicube_Texel> texels;
Array <Hemicube_Edge> edges;
};
Raw
hemicube_integrate.cpp
#include "hemicube.h"
Hemicube_Data hemicube_data;
struct Irradiance_Sample {
Vector3 position;
Basis3 frame;
float horizon_cos;
bool is_valid_irradiance;
Vector3 color;
Vector3 rotation_gradients[3];
Vector3 translation_gradients[3];
};
void integrate(Irradiance_Sample * irradiance_sample, const float * hemicube_ptr)
{
Array<Vector2> solid_angle_gradient;
// Sum edge gradients.
solid_angle_gradient.resize(hemicube_data.width*hemicube_data.height, Vector2(0));
const int hemicube_edge_count = hemicube_data.edges.count();
for (int i = 0; i < hemicube_edge_count; i++) {
const Hemicube_Edge & edge = hemicube_data.edges[i];
const float depth0 = hemicube_ptr[edge.y0 * hemicube_pitch + edge.x0 * 4 + 3];
const float depth1 = hemicube_ptr[edge.y1 * hemicube_pitch + edge.x1 * 4 + 3];
const float min_depth = min(depth0, depth1);
if (min_depth <= 0) continue; // Ignore interior samples.
const float inv_min_depth = 1.0f / fabs(min_depth);
solid_angle_gradient[edge.y0 * hemicube_data.width + edge.x0] -= edge.gradient * inv_min_depth;
solid_angle_gradient[edge.y1 * hemicube_data.width + edge.x1] += edge.gradient * inv_min_depth;
}
float visible_area = 0.0f;
float total_area = 0.0f;
Vector3 color_sum(0);
float depth_min = FLT_MAX;
float depth_inv_sum = 0;
Vector2 rotation_gradient_sum[3] = { Vector2(0), Vector2(0), Vector2(0) };
Vector2 translation_gradient_sum[3] = { Vector2(0), Vector2(0), Vector2(0) };
// Make threshold dependent on angle between face normal and interpolated vertex normal.
float horizon_cos = irradiance_sample->horizon_cos;
for (int y = 0; y < hemicube_data.height; y++) {
const float * scanline = hemicube_ptr + y * hemicube_pitch;
for (int x = 0; x < hemicube_data.width; x++) {
const Hemicube_Texel & texel = hemicube_data.texels[y * hemicube_data.width + x];
int face = texel.face;
if (face == -1) continue;
const float * rgbd = scanline + x*4;
const Vector3 color = Vector3(rgbd[0], rgbd[1], rgbd[2]);
const float depth = rgbd[3];
const float texel_factor = texel.clamped_cosine * texel.solid_angle;
total_area += texel_factor;
// Ignore interior faces.
if (depth > 0) {
visible_area += texel_factor;
color_sum += color * texel_factor;
}
// Depth min and harmonic mean:
const float abs_depth = min(1024.0f, fabsf(depth));
depth_inv_sum += 1.0f / abs_depth; // to compute harmonic mean
if (texel.clamped_cosine > horizon_cos) { // Do not take into account samples near the tangent plane.
depth_min = min(depth_min, abs_depth);
}
// Ignore interior faces.
if (depth > 0) {
Vector2 rotation_vector = Vector2(texel.dir.y, -texel.dir.x) * texel.solid_angle;
rotation_gradient_sum[0] += rotation_vector * color.x;
rotation_gradient_sum[1] += rotation_vector * color.y;
rotation_gradient_sum[2] += rotation_vector * color.z;
Vector2 translation_gradient = solid_angle_gradient[y * hemicube_data.width + x] * texel.clamped_cosine;
translation_gradient_sum[0] += translation_gradient * color.x;
translation_gradient_sum[1] += translation_gradient * color.y;
translation_gradient_sum[2] += translation_gradient * color.z;
}
}
}
const float normalization_factor = (1.0f / total_area);
const float confidence = visible_area * normalization_factor;
const float confidence_threshold = 0.9f;
// Update irradiance sample.
const bool is_valid_sample = (confidence > confidence_threshold);
irradiance_sample->is_valid_irradiance = is_valid_sample;
if (is_valid_sample) {
irradiance_sample->color = color_sum * normalization_factor;
irradiance_sample->depth_min = depth_min;
irradiance_sample->depth_harmonic_mean = hemicube_data.cube_size * hemicube_data.cube_size * 3 / depth_inv_sum;
// Transform rotation gradients to object space.
for (int i = 0; i < 3; i++) {
irradiance_sample->rotation_gradients[i] = irradiance_sample->frame.transform(Vector3(rotation_gradient_sum[i], 0.0f));
irradiance_sample->rotation_gradients[i] *= 1.0f / PI;
}
float max_length = -FLT_MAX;
for (int i = 0; i < 3; i++) {
irradiance_sample->translation_gradients[i] = irradiance_sample->frame.transform(Vector3(translation_gradient_sum[i], 0.0f));
irradiance_sample->translation_gradients[i] *= 1.0f / PI;
max_length = max(max_length, lengthSquared(irradiance_sample->translation_gradients[i]));
}
// Clamp gradients without modifying the hue.
max_length = max_gradient_length / max(max_gradient_length, sqrtf(max_length));
for (int i = 0; i < 3; i++) {
irradiance_sample->translation_gradients[i] *= max_length;
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment