ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/dircode.c
Revision: 2.11
Committed: Wed Mar 4 02:55:43 2020 UTC (4 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.10: +6 -6 lines
Log Message:
Made encodings for X, Y, and Z-plane normals exact

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.11 static const char RCSid[] = "$Id: dircode.c,v 2.10 2019/08/27 23:51:23 greg Exp $";
3 greg 2.1 #endif
4     /*
5 greg 2.7 * Compute a 4-byte direction code (externals defined in rtmath.h).
6 greg 2.3 *
7     * Mean accuracy is 0.0022 degrees, with a maximum error of 0.0058 degrees.
8 greg 2.1 */
9    
10 greg 2.7 #include "rtmath.h"
11 greg 2.1
12 greg 2.11 #define DCSCALE 11584.7 /* (1<<13)*sqrt(2) - .5 */
13 greg 2.1 #define FXNEG 01
14     #define FYNEG 02
15     #define FZNEG 04
16     #define F1X 010
17     #define F2Z 020
18     #define F1SFT 5
19     #define F2SFT 18
20     #define FMASK 0x1fff
21    
22 greg 2.5 int32
23 greg 2.8 encodedir(FVECT dv) /* encode a normalized direction vector */
24 greg 2.1 {
25 greg 2.8 int32 dc = 0;
26 greg 2.1 int cd[3], cm;
27 greg 2.8 int i;
28 greg 2.1
29     for (i = 0; i < 3; i++)
30     if (dv[i] < 0.) {
31 greg 2.11 cd[i] = (int)(dv[i] * -DCSCALE + .5);
32 greg 2.1 dc |= FXNEG<<i;
33     } else
34 greg 2.11 cd[i] = (int)(dv[i] * DCSCALE + .5);
35 greg 2.6 if (!(cd[0] | cd[1] | cd[2]))
36     return(0); /* zero normal */
37 greg 2.1 if (cd[0] <= cd[1]) {
38     dc |= F1X | cd[0] << F1SFT;
39     cm = cd[1];
40     } else {
41     dc |= cd[1] << F1SFT;
42     cm = cd[0];
43     }
44     if (cd[2] <= cm)
45     dc |= F2Z | cd[2] << F2SFT;
46     else
47     dc |= cm << F2SFT;
48 greg 2.6 if (!dc) /* don't generate 0 code normally */
49 greg 2.2 dc = F1X;
50 greg 2.1 return(dc);
51     }
52    
53 greg 2.9 #if 0 /* original version for reference */
54 greg 2.1
55     void
56 greg 2.8 decodedir(FVECT dv, int32 dc) /* decode a normalized direction vector */
57 greg 2.1 {
58     double d1, d2, der;
59    
60 greg 2.6 if (!dc) { /* special code for zero normal */
61     dv[0] = dv[1] = dv[2] = 0.;
62     return;
63     }
64 greg 2.1 d1 = ((dc>>F1SFT & FMASK)+.5)*(1./DCSCALE);
65     d2 = ((dc>>F2SFT & FMASK)+.5)*(1./DCSCALE);
66     der = sqrt(1. - d1*d1 - d2*d2);
67     if (dc & F1X) {
68     dv[0] = d1;
69     if (dc & F2Z) { dv[1] = der; dv[2] = d2; }
70     else { dv[1] = d2; dv[2] = der; }
71     } else {
72     dv[1] = d1;
73     if (dc & F2Z) { dv[0] = der; dv[2] = d2; }
74     else { dv[0] = d2; dv[2] = der; }
75     }
76     if (dc & FXNEG) dv[0] = -dv[0];
77     if (dc & FYNEG) dv[1] = -dv[1];
78     if (dc & FZNEG) dv[2] = -dv[2];
79     }
80    
81 greg 2.9 #else
82    
83     void
84     decodedir(FVECT dv, int32 dc) /* decode a normalized direction vector */
85     {
86     static const short itab[4][3] = {
87     {1,0,2},{0,1,2},{1,2,0},{0,2,1}
88     };
89     static const RREAL neg[2] = {1., -1.};
90     const int ndx = ((dc & F2Z) != 0)<<1 | ((dc & F1X) != 0);
91     double d1, d2, der;
92    
93     if (!dc) { /* special code for zero normal */
94     dv[0] = dv[1] = dv[2] = 0.;
95     return;
96     }
97 greg 2.11 d1 = (dc>>F1SFT & FMASK)*(1./DCSCALE);
98     d2 = (dc>>F2SFT & FMASK)*(1./DCSCALE);
99 greg 2.9 der = sqrt(1. - d1*d1 - d2*d2);
100     dv[itab[ndx][0]] = d1;
101     dv[itab[ndx][1]] = d2;
102     dv[itab[ndx][2]] = der;
103     dv[0] *= neg[(dc&FXNEG)!=0];
104     dv[1] *= neg[(dc&FYNEG)!=0];
105     dv[2] *= neg[(dc&FZNEG)!=0];
106     }
107    
108     #endif
109 greg 2.1
110     double
111 greg 2.8 dir2diff(int32 dc1, int32 dc2) /* approx. radians^2 between directions */
112 greg 2.1 {
113     FVECT v1, v2;
114    
115 greg 2.10 if (dc1 == dc2)
116     return 0.;
117    
118 greg 2.1 decodedir(v1, dc1);
119     decodedir(v2, dc2);
120    
121     return(2. - 2.*DOT(v1,v2));
122     }
123    
124     double
125 greg 2.8 fdir2diff(int32 dc1, FVECT v2) /* approx. radians^2 between directions */
126 greg 2.1 {
127     FVECT v1;
128    
129     decodedir(v1, dc1);
130    
131     return(2. - 2.*DOT(v1,v2));
132     }