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