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