Tuesday, January 18, 2011

Improved code for concentric map

If you need to warp points on a square to a disk, many people use Ken Chiu's and my code from jgt:

/* seedx, seedy is point on [0,1]^2. x, y is point on radius 1 disk */
void to_unit_disk( double seedx, double seedy, double *x, double *y )
{

double phi, r;

double a = 2*seedx - 1; >/* (a,b) is now on [-1,1]^2 */
double b = 2*seedy - 1;

if (a > -b) { /* region 1 or 2 */
if (a > b) { /* region 1, also |a| > |b| */
r = a;
phi = (M_PI/4 ) * (b/a);
}
else { /* region 2, also |b| > |a| */
phi = (M_PI/4) * (2 - (a/b));
}
}
else { /* region 3 or 4 */
if (a < b) /* region 3, also |a| >= |b|, a != 0
r = -a;
phi = (M_PI/4) * (4 + (b/a));
}
else { /* region 4, |b| >= |a|, but a==0 and b==0 could occur.
r = -b;
if (b != 0)
phi = (M_PI/4) * (6 - (a/b));
else
phi = 0;
}
}

*x = r * cos(phi);
*y = r * sin(phi);

}

Dave Cline recently sent me a neat trick that uses negative radii and I think is correct. Let me know if you try it. (cut and pasted from his mail)

Vector2 ToUnitDisk(Vector2 O) {
float phi,r;
float a = 2*O.x - 1;
float b = 2*O.y - 1;
if (a*a> b*b) { // use squares instead of absolute values
r = a;
phi = (PI/4)*(b/a);
} else {
r = b;
phi = (PI/4)*(a/b) + (PI/2);
}
return Vector2( r*cos(phi), r*sin(phi) );
}

2 comments:

franz said...

Hey,

As is, the improved formulation leads to clumping and alignment issues when fed samples coming from a Hammersley sequence. It's easy to fix: in the second branch of the if, replace

phi = (PI/4)*(a/b) + (PI/2);

by

phi = (PI/2) - (PI/4)*(a/b);

With this modification, the new function will return the exact same results as the original jgt one.

Cheers,
Franz

Greg Ward said...

Nice! Looks like we missed the case where (a==0&&b==0), so we'll have a divide-by-zero in the else section unless we add a test.

BTW, I'm using this heavily in my new BSDF sampling scheme. Very handy!