| 1 |
/*
|
| 2 |
* Disk2Square.c
|
| 3 |
*
|
| 4 |
* Code by Peter Shirley and Kenneth Chiu
|
| 5 |
*
|
| 6 |
* Associated paper in ~/Documents/Others/Shirley+Chiu_JGT-97.pdf
|
| 7 |
*
|
| 8 |
* Modified interface slightly (G. Ward)
|
| 9 |
*/
|
| 10 |
|
| 11 |
#include <math.h>
|
| 12 |
|
| 13 |
/*
|
| 14 |
This transforms points on [0,1]^2 to points on unit disk centered at
|
| 15 |
origin. Each "pie-slice" quadrant of square is handled as a seperate
|
| 16 |
case. The bad floating point cases are all handled appropriately.
|
| 17 |
The regions for (a,b) are:
|
| 18 |
|
| 19 |
phi = pi/2
|
| 20 |
-----*-----
|
| 21 |
|\ /|
|
| 22 |
| \ 2 / |
|
| 23 |
| \ / |
|
| 24 |
phi=pi* 3 * 1 *phi = 0
|
| 25 |
| / \ |
|
| 26 |
| / 4 \ |
|
| 27 |
|/ \|
|
| 28 |
-----*-----
|
| 29 |
phi = 3pi/2
|
| 30 |
|
| 31 |
change log:
|
| 32 |
26 feb 2004. bug fix in computation of phi (now this matches the paper,
|
| 33 |
which is correct). thanks to Atilim Cetin for catching this.
|
| 34 |
*/
|
| 35 |
|
| 36 |
|
| 37 |
/* Map a [0,1]^2 square to a unit radius disk */
|
| 38 |
void
|
| 39 |
SDsquare2disk(double ds[2], double seedx, double seedy)
|
| 40 |
{
|
| 41 |
|
| 42 |
double phi, r;
|
| 43 |
double a = 2.*seedx - 1; /* (a,b) is now on [-1,1]^2 */
|
| 44 |
double b = 2.*seedy - 1;
|
| 45 |
|
| 46 |
if (a > -b) { /* region 1 or 2 */
|
| 47 |
if (a > b) { /* region 1, also |a| > |b| */
|
| 48 |
r = a;
|
| 49 |
phi = (M_PI/4.) * (b/a);
|
| 50 |
}
|
| 51 |
else { /* region 2, also |b| > |a| */
|
| 52 |
r = b;
|
| 53 |
phi = (M_PI/4.) * (2. - (a/b));
|
| 54 |
}
|
| 55 |
}
|
| 56 |
else { /* region 3 or 4 */
|
| 57 |
if (a < b) { /* region 3, also |a| >= |b|, a != 0 */
|
| 58 |
r = -a;
|
| 59 |
phi = (M_PI/4.) * (4. + (b/a));
|
| 60 |
}
|
| 61 |
else { /* region 4, |b| >= |a|, but a==0 and b==0 could occur. */
|
| 62 |
r = -b;
|
| 63 |
if (b != 0.)
|
| 64 |
phi = (M_PI/4.) * (6. - (a/b));
|
| 65 |
else
|
| 66 |
phi = 0.;
|
| 67 |
}
|
| 68 |
}
|
| 69 |
|
| 70 |
ds[0] = r * cos(phi);
|
| 71 |
ds[1] = r * sin(phi);
|
| 72 |
|
| 73 |
}
|
| 74 |
|
| 75 |
/* Map point on unit disk to a unit square in [0,1]^2 range */
|
| 76 |
void
|
| 77 |
SDdisk2square(double sq[2], double diskx, double disky)
|
| 78 |
{
|
| 79 |
double r = sqrt( diskx*diskx + disky*disky );
|
| 80 |
double phi = atan2( disky, diskx );
|
| 81 |
double a, b;
|
| 82 |
|
| 83 |
if (phi < -M_PI/4) phi += 2*M_PI; /* in range [-pi/4,7pi/4] */
|
| 84 |
if ( phi < M_PI/4) { /* region 1 */
|
| 85 |
a = r;
|
| 86 |
b = phi * a / (M_PI/4);
|
| 87 |
}
|
| 88 |
else if ( phi < 3*M_PI/4 ) { /* region 2 */
|
| 89 |
b = r;
|
| 90 |
a = -(phi - M_PI/2) * b / (M_PI/4);
|
| 91 |
}
|
| 92 |
else if ( phi < 5*M_PI/4 ) { /* region 3 */
|
| 93 |
a = -r;
|
| 94 |
b = (phi - M_PI) * a / (M_PI/4);
|
| 95 |
}
|
| 96 |
else { /* region 4 */
|
| 97 |
b = -r;
|
| 98 |
a = -(phi - 3*M_PI/2) * b / (M_PI/4);
|
| 99 |
}
|
| 100 |
|
| 101 |
sq[0] = (a + 1) * 0.5;
|
| 102 |
sq[1] = (b + 1) * 0.5;
|
| 103 |
}
|
