ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/disk2square.cal
Revision: 2.5
Committed: Tue Aug 13 18:29:23 2019 UTC (3 years, 3 months ago) by greg
Branch: MAIN
CVS Tags: rad5R3, HEAD
Changes since 2.4: +4 -3 lines
Log Message:
Improved handling of edge cases

File Contents

# User Rev Content
1 greg 2.5 { RCSid $Id: disk2square.cal,v 2.4 2015/03/27 18:58:06 greg Exp $ }
2 greg 2.1 {
3     Convert between unit square and disk, using Shirley-Chiu mapping
4    
5     5/20/2011 G. Ward
6    
7     Inputs:
8     in_disk_x Input unit disk Cartesian coordiantes
9     in_disk_y center at (0,0)
10     or:
11     in_square_x Input unit square [0,1]^2
12     in_square_y
13    
14     Corresponding outputs:
15     out_square_x Output unit square [0,1]^2
16     out_square_y
17     or:
18     out_disk_x Output unit disk Cartesian coordinates
19     out_disk_y center at (0,0)
20     }
21     { Compute disk position from square coordinates }
22     in_square_a = 2*in_square_x - 1;
23     in_square_b = 2*in_square_y - 1;
24    
25     in_square_rgn = if(in_square_a + in_square_b,
26     if(in_square_a - in_square_b, 1, 2),
27     if(in_square_b - in_square_a, 3, 4));
28    
29     out_disk_r = select(in_square_rgn, in_square_a, in_square_b,
30     -in_square_a, -in_square_b);
31    
32     out_disk_phi = PI/4 * select(in_square_rgn,
33     in_square_b/in_square_a,
34     2 - in_square_a/in_square_b,
35     4 + in_square_b/in_square_a,
36     if(in_square_b*in_square_b,
37     6 - in_square_a/in_square_b,
38     0));
39    
40     out_disk_x = out_disk_r*cos(out_disk_phi);
41     out_disk_y = out_disk_r*sin(out_disk_phi);
42    
43     { Compute square position from disk coordinates }
44     norm_radians(p) : if(-p - PI/4, p + 2*PI, p);
45     in_disk_r = sqrt(in_disk_x*in_disk_x + in_disk_y*in_disk_y);
46     in_disk_phi = norm_radians(atan2(in_disk_y, in_disk_x));
47    
48     in_disk_rgn = floor((in_disk_phi + PI/4)/(PI/2)) + 1;
49    
50     out_square_a = select(in_disk_rgn,
51     in_disk_r,
52     (PI/2 - in_disk_phi)*in_disk_r/(PI/4),
53     -in_disk_r,
54 greg 2.5 (in_disk_phi - 3*PI/2)*in_disk_r/(PI/4),
55     in_disk_r);
56 greg 2.1
57     out_square_b = select(in_disk_rgn,
58     in_disk_phi*in_disk_r/(PI/4),
59     in_disk_r,
60     (PI - in_disk_phi)*in_disk_r/(PI/4),
61 greg 2.5 -in_disk_r, -in_disk_r);
62 greg 2.1
63     out_square_x = (out_square_a + 1)/2;
64     out_square_y = (out_square_b + 1)/2;
65     {
66     The following forumulas compute Shirley-Chiu bin "scbin" based on:
67    
68 greg 2.4 RHS - right-handed system (-1 for left-handed coords)
69 greg 2.1 Dx,Dy,Dz - Incident direction (normalized, towards surface front)
70 greg 2.3 rNx,rNy,rNz - Surface normal (normalized, away from surface)
71 greg 2.1 Ux,Uy,Uz - Up direction vector (does not need to be normalized)
72    
73     The SCdim variable assigns the square side dimension for bins, which are
74     ordered with the "up" direction changing fastest.
75     }
76 greg 2.4 RHS = 1;
77 greg 2.1 { Compute oriented axis values }
78 greg 2.3 inc_dz = -Dx*rNx-Dy*rNy-Dz*rNz;
79 greg 2.4 inc_rx = -RHS*(Dx*(Uy*rNz-Uz*rNy) + Dy*(Uz*rNx-Ux*rNz) + Dz*(Ux*rNy-Uy*rNx));
80 greg 2.3 inc_ry = -Dx*Ux-Dy*Uy-Dz*Uz - inc_dz*(rNx*Ux+rNy*Uy+rNz*Uz);
81 greg 2.1 inc_den2 = inc_rx*inc_rx + inc_ry*inc_ry;
82 greg 2.3 inc_radf = if(inc_den2-1e-7, sqrt((1 - inc_dz*inc_dz)/inc_den2), 0);
83 greg 2.1 { Pass to formulas in first section }
84     in_disk_x = inc_rx*inc_radf;
85     in_disk_y = inc_ry*inc_radf;
86     { Compute final bin (-1 if behind surface) }
87 greg 2.2 scbin = if(inc_dz,
88 greg 2.1 floor(out_square_x*SCdim)*SCdim + floor(out_square_y*SCdim),
89     -1);