src/rt/rayinit.cal

{ RCSid $Id: rayinit.cal,v 2.18 2018/11/30 18:43:57 greg Exp $ }
{
        Initialization file for Radiance.

        The following are predefined:

        Dx, Dy, Dz                      - ray direction
        Nx, Ny, Nz                      - surface normal
        Px, Py, Pz                      - intersection point
        T                               - distance from start
        Ts                              - single ray (shadow) distance
        Rdot                            - ray dot product
        S                               - world scale
        Tx, Ty, Tz                      - world origin
        Ix, Iy, Iz                      - world i unit vector
        Jx, Jy, Jz                      - world j unit vector
        Kx, Ky, Kz                      - world k unit vector
        arg(n)                          - real arguments, arg(0) is count

        For mesh objects, the following are available:

        Lu, Lv                          - local (u,v) coordinates

        For *func & *data materials, the following are also available:

        NxP, NyP, NzP                   - perturbed surface normal
        RdotP                           - perturbed ray dot product
        CrP, CgP, CbP                   - perturbed material color

        For prism1 and prism2 types, the following are available:

        DxA, DyA, DzA                   - direction to target light source

        Library functions:

        if(a, b, c)                     - if a positive, return b, else c

        select(N, a1, a2, ..)           - return aN

        min(a1, a2, ..)                 - return minimum argument
        max(a1, a2, ..)                 - return maximum argument

        sqrt(x)                         - square root function

        sin(x), cos(x), tan(x),
        asin(x), acos(x),
        atan(x), atan2(y,x)             - standard trig functions (radians)

        floor(x), ceil(x)               - g.l.b. & l.u.b.

        exp(x), log(x), log10(x)        - exponent and log functions

        erf(z), erfc(z)                 - error functions

        rand(x)                         - pseudo-random function (0 to 1)

        noise3(x,y,z), noise3x(x,y,z),
        noise3y(x,y,z), noise3z(x,y,z)  - noise function with gradient (-1 to 1)

        fnoise3(x,y,z)                  - fractal noise function (-1 to 1)
}

                        { Backward compatibility }
AC = arg(0);
A1 = arg(1); A2 = arg(2); A3 = arg(3); A4 = arg(4); A5 = arg(5);
A6 = arg(6); A7 = arg(7); A8 = arg(8); A9 = arg(9); A10 = arg(10);

noise3a(x,y,z) : noise3x(x,y,z);
noise3b(x,y,z) : noise3y(x,y,z);
noise3c(x,y,z) : noise3z(x,y,z);

                        { Forward compatibility (?) }
D(i) = select(i, Dx, Dy, Dz);
N(i) = select(i, Nx, Ny, Nz);
P(i) = select(i, Px, Py, Pz);
noise3d(i,x,y,z) : select(i, noise3x(x,y,z), noise3y(x,y,z), noise3z(x,y,z));

                        { More robust versions of library functions }
bound(a,x,b) : if(a-x, a, if(x-b, b, x));
Acos(x) : acos(bound(-1,x,1));
Asin(x) : asin(bound(-1,x,1));
Atan2(y,x) : if(x*x+y*y, atan2(y,x), 0);
Exp(x) : if(-x-100, 0, exp(x));
Sqrt(x) : if(x, sqrt(x), 0);

                        { Useful constants }
PI : 3.14159265358979323846;
DEGREE : PI/180;
FTINY : 1e-7;

                        { Useful functions }
and(a,b) : if( a, b, a );
or(a,b) : if( a, a, b );
not(a) : if( a, -1, 1 );
xor(a,b) : if( a, not(b), b );
abs(x) : if( x, x, -x );
sgn(x) : if( x, 1, if(-x, -1, 0) );
sq(x) : x*x;
inside(a,x,b) : and(x-a,b-x);
frac(x) : x - floor(x);
mod(n,d) : n - floor(n/d)*d;
tri(n,d) : abs( d - mod(n-d,2*d) );
linterp(t,p0,p1) : (1-t)*p0 + t*p1;

noop(v) : v;
clip(v) : bound(0,v,1);
noneg(v) : if(v,v,0);
red(r,g,b) : if(r,r,0);
green(r,g,b) : if(g,g,0);
blue(r,g,b) : if(b,b,0);
grey(r,g,b) : noneg(.265074126*r + .670114631*g + .064811243*b);
clip_r(r,g,b) : bound(0,r,1);
clip_g(r,g,b) : bound(0,g,1);
clip_b(r,g,b) : bound(0,b,1);
clipgrey(r,g,b) : min(grey(r,g,b),1);

dot(v1,v2) : v1(1)*v2(1) + v1(2)*v2(2) + v1(3)*v2(3);
cross(i,v1,v2) : select(i,      v1(2)*v2(3) - v1(3)*v2(2),
                                v1(3)*v2(1) - v1(1)*v2(3),
                                v1(1)*v2(2) - v1(2)*v2(1));

fade(near_val,far_val,dist) : far_val +
                if (16-dist, (near_val-far_val)/(1+dist*dist), 0);

hermite(p0,p1,r0,r1,t) :        p0 * ((2*t-3)*t*t+1) +
                                p1 * (-2*t+3)*t*t +
                                r0 * (((t-2)*t+1)*t) +
                                r1 * ((t-1)*t*t);

bezier(p1, p2, p3, p4, t) :     p1 * (1+t*(-3+t*(3-t))) +
                                p2 * 3*t*(1+t*(-2+t)) +
                                p3 * 3*t*t*(1-t) +
                                p4 * t*t*t ;

bspline(pp, p0, p1, pn, t) :    pp * (1/6+t*(-.5+t*(.5-1/6*t))) +
                                p0 * (2/3+t*t*(-1+.5*t)) +
                                p1 * (1/6+t*(.5+t*(.5-.5*t))) +
                                pn * (1/6*t*t*t) ;

turbulence(x,y,z,s) : if( s-1.01, 0, abs(noise3(x/s,y/s,z/s)*s) +
                                                turbulence(x,y,z,2*s) );
turbulencex(x,y,z,s) : if( s-1.01, 0,
                        sgn(noise3(x/s,y/s,z/s))*noise3x(x/s,y/s,z/s) +
                        turbulencex(x,y,z,2*s) );
turbulencey(x,y,z,s) : if( s-1.01, 0,
                        sgn(noise3(x/s,y/s,z/s))*noise3y(x/s,y/s,z/s) +
                        turbulencey(x,y,z,2*s) );
turbulencez(x,y,z,s) : if( s-1.01, 0,
                        sgn(noise3(x/s,y/s,z/s))*noise3z(x/s,y/s,z/s) +
                        turbulencez(x,y,z,2*s) );

                        { Normal distribution from uniform range (0,1) }

un2`P.(t) : t - (2.515517+t*(.802853+t*.010328))/
                (1+t*(1.432788+t*(.189269+t*.001308))) ;
un1`P.(p) : un2`P.(sqrt(-2*log(p))) ;

unif2norm(p) : if( .5-p, -un1`P.(p), un1`P.(1-p) ) ;

nrand(x) = unif2norm(rand(x));

                        { Local (u,v) coordinates for planar surfaces }
crosslen`P. = Nx*Nx + Ny*Ny;
                        { U is distance from projected Z-axis }
U = if( crosslen`P. - FTINY,
                (Py*Nx - Px*Ny)/crosslen`P.,
                Px);
                        { V is defined so that N = U x V }
V = if( crosslen`P. - FTINY,
                Pz - Nz*(Px*Nx + Py*Ny)/crosslen`P.,
                Py);

                        { Local hemisphere direction for *func & *data types }
                        { last 3 real args = unnormalized up-vector }
Vux`P. = arg(AC-1)*NzP - arg(AC)*NyP;
Vuy`P. = arg(AC)*NxP - arg(AC-2)*NzP;
Vuz`P. = arg(AC-2)*NyP - arg(AC-1)*NxP;
vnorm`P. = 1/sqrt(Vux`P.*Vux`P. + Vuy`P.*Vuy`P. + Vuz`P.*Vuz`P.);
Vnx`P. = Vux`P.*vnorm`P.;
Vny`P. = Vuy`P.*vnorm`P.;
Vnz`P. = Vuz`P.*vnorm`P.;
Unx`P. = NyP*Vnz`P. - NzP*Vny`P.;
Uny`P. = NzP*Vnx`P. - NxP*Vnz`P.;
Unz`P. = NxP*Vny`P. - NyP*Vnx`P.;
                        { Transform vectors, normalized (dx,dy,dz) away from surf }
Ldx(dx,dy,dz) = dx*Unx`P. + dy*Uny`P. + dz*Unz`P.;
Ldy(dx,dy,dz) = dx*Vnx`P. + dy*Vny`P. + dz*Vnz`P.;
Ldz(dx,dy,dz) = dx*NxP + dy*NyP + dz*NzP;
                        { Incident vector transformed to our coords }
Idx = Ldx(-Dx,-Dy,-Dz);
Idy = Ldy(-Dx,-Dy,-Dz);
Idz = RdotP;
Revision:	2.19
Committed:	Mon Jun 10 13:56:52 2019 UTC (4 years, 10 months ago) by greg
Branch:	MAIN
CVS Tags:	rad5R4, rad5R3, HEAD
Changes since 2.18:	+4 -3 lines
Log Message:	Added max() and min() as .cal library functions
#	User	Rev	Content
1	greg	2.19	{ RCSid $Id: rayinit.cal,v 2.18 2018/11/30 18:43:57 greg Exp $ }
2	greg	1.1	{
3			Initialization file for Radiance.
4
5			The following are predefined:
6
7			Dx, Dy, Dz - ray direction
8			Nx, Ny, Nz - surface normal
9			Px, Py, Pz - intersection point
10			T - distance from start
11	greg	2.3	Ts - single ray (shadow) distance
12	greg	1.1	Rdot - ray dot product
13			S - world scale
14			Tx, Ty, Tz - world origin
15			Ix, Iy, Iz - world i unit vector
16			Jx, Jy, Jz - world j unit vector
17			Kx, Ky, Kz - world k unit vector
18			arg(n) - real arguments, arg(0) is count
19	greg	2.14
20			For mesh objects, the following are available:
21
22			Lu, Lv - local (u,v) coordinates
23	greg	1.1
24	greg	2.15	For func & data materials, the following are also available:
25	greg	1.1
26			NxP, NyP, NzP - perturbed surface normal
27			RdotP - perturbed ray dot product
28			CrP, CgP, CbP - perturbed material color
29
30	greg	2.9	For prism1 and prism2 types, the following are available:
31
32			DxA, DyA, DzA - direction to target light source
33
34	greg	1.1	Library functions:
35
36			if(a, b, c) - if a positive, return b, else c
37
38			select(N, a1, a2, ..) - return aN
39
40	greg	2.19	min(a1, a2, ..) - return minimum argument
41			max(a1, a2, ..) - return maximum argument
42
43	greg	1.1	sqrt(x) - square root function
44
45			sin(x), cos(x), tan(x),
46			asin(x), acos(x),
47	greg	2.15	atan(x), atan2(y,x) - standard trig functions (radians)
48	greg	1.1
49			floor(x), ceil(x) - g.l.b. & l.u.b.
50
51			exp(x), log(x), log10(x) - exponent and log functions
52
53			erf(z), erfc(z) - error functions
54
55			rand(x) - pseudo-random function (0 to 1)
56
57	greg	2.11	noise3(x,y,z), noise3x(x,y,z),
58			noise3y(x,y,z), noise3z(x,y,z) - noise function with gradient (-1 to 1)
59	greg	1.1
60			fnoise3(x,y,z) - fractal noise function (-1 to 1)
61			}
62
63			{ Backward compatibility }
64			AC = arg(0);
65			A1 = arg(1); A2 = arg(2); A3 = arg(3); A4 = arg(4); A5 = arg(5);
66			A6 = arg(6); A7 = arg(7); A8 = arg(8); A9 = arg(9); A10 = arg(10);
67
68	greg	2.11	noise3a(x,y,z) : noise3x(x,y,z);
69			noise3b(x,y,z) : noise3y(x,y,z);
70			noise3c(x,y,z) : noise3z(x,y,z);
71
72	greg	1.1	{ Forward compatibility (?) }
73			D(i) = select(i, Dx, Dy, Dz);
74			N(i) = select(i, Nx, Ny, Nz);
75			P(i) = select(i, Px, Py, Pz);
76	greg	2.11	noise3d(i,x,y,z) : select(i, noise3x(x,y,z), noise3y(x,y,z), noise3z(x,y,z));
77	greg	1.1
78			{ More robust versions of library functions }
79			bound(a,x,b) : if(a-x, a, if(x-b, b, x));
80			Acos(x) : acos(bound(-1,x,1));
81			Asin(x) : asin(bound(-1,x,1));
82	greg	2.8	Atan2(y,x) : if(xx+yy, atan2(y,x), 0);
83	greg	2.2	Exp(x) : if(-x-100, 0, exp(x));
84	greg	1.1	Sqrt(x) : if(x, sqrt(x), 0);
85
86			{ Useful constants }
87			PI : 3.14159265358979323846;
88			DEGREE : PI/180;
89			FTINY : 1e-7;
90
91			{ Useful functions }
92			and(a,b) : if( a, b, a );
93			or(a,b) : if( a, a, b );
94			not(a) : if( a, -1, 1 );
95	greg	2.13	xor(a,b) : if( a, not(b), b );
96	greg	1.1	abs(x) : if( x, x, -x );
97			sgn(x) : if( x, 1, if(-x, -1, 0) );
98			sq(x) : x*x;
99			inside(a,x,b) : and(x-a,b-x);
100			frac(x) : x - floor(x);
101			mod(n,d) : n - floor(n/d)*d;
102			tri(n,d) : abs( d - mod(n-d,2*d) );
103			linterp(t,p0,p1) : (1-t)p0 + tp1;
104
105	gwlarson	2.12	noop(v) : v;
106			clip(v) : bound(0,v,1);
107			noneg(v) : if(v,v,0);
108			red(r,g,b) : if(r,r,0);
109			green(r,g,b) : if(g,g,0);
110			blue(r,g,b) : if(b,b,0);
111			grey(r,g,b) : noneg(.265074126r + .670114631g + .064811243*b);
112			clip_r(r,g,b) : bound(0,r,1);
113			clip_g(r,g,b) : bound(0,g,1);
114			clip_b(r,g,b) : bound(0,b,1);
115			clipgrey(r,g,b) : min(grey(r,g,b),1);
116	greg	1.1
117			dot(v1,v2) : v1(1)v2(1) + v1(2)v2(2) + v1(3)*v2(3);
118			cross(i,v1,v2) : select(i, v1(2)v2(3) - v1(3)v2(2),
119			v1(3)v2(1) - v1(1)v2(3),
120			v1(1)v2(2) - v1(2)v2(1));
121
122	gwlarson	2.12	fade(near_val,far_val,dist) : far_val +
123	greg	1.1	if (16-dist, (near_val-far_val)/(1+dist*dist), 0);
124
125	greg	2.16	hermite(p0,p1,r0,r1,t) : p0 * ((2t-3)t*t+1) +
126			p1 * (-2t+3)t*t +
127			r0 * (((t-2)t+1)t) +
128			r1 * ((t-1)tt);
129
130	gwlarson	2.12	bezier(p1, p2, p3, p4, t) : p1 * (1+t(-3+t(3-t))) +
131	greg	1.1	p2 * 3t(1+t*(-2+t)) +
132			p3 * 3tt*(1-t) +
133			p4 * ttt ;
134
135	gwlarson	2.12	bspline(pp, p0, p1, pn, t) : pp * (1/6+t(-.5+t(.5-1/6*t))) +
136	greg	1.1	p0 * (2/3+tt(-1+.5*t)) +
137			p1 * (1/6+t(.5+t(.5-.5*t))) +
138			pn * (1/6tt*t) ;
139
140	gwlarson	2.12	turbulence(x,y,z,s) : if( s-1.01, 0, abs(noise3(x/s,y/s,z/s)*s) +
141	greg	1.1	turbulence(x,y,z,2*s) );
142	gwlarson	2.12	turbulencex(x,y,z,s) : if( s-1.01, 0,
143	greg	2.11	sgn(noise3(x/s,y/s,z/s))*noise3x(x/s,y/s,z/s) +
144	gwlarson	2.12	turbulencex(x,y,z,2*s) );
145			turbulencey(x,y,z,s) : if( s-1.01, 0,
146	greg	2.11	sgn(noise3(x/s,y/s,z/s))*noise3y(x/s,y/s,z/s) +
147	gwlarson	2.12	turbulencey(x,y,z,2*s) );
148			turbulencez(x,y,z,s) : if( s-1.01, 0,
149	greg	2.11	sgn(noise3(x/s,y/s,z/s))*noise3z(x/s,y/s,z/s) +
150	gwlarson	2.12	turbulencez(x,y,z,2*s) );
151	greg	1.2
152	greg	1.3	{ Normal distribution from uniform range (0,1) }
153
154	greg	2.17	un2`P.(t) : t - (2.515517+t(.802853+t.010328))/
155	greg	1.3	(1+t(1.432788+t(.189269+t*.001308))) ;
156	greg	2.17	un1`P.(p) : un2`P.(sqrt(-2*log(p))) ;
157	greg	1.3
158	greg	2.17	unif2norm(p) : if( .5-p, -un1`P.(p), un1`P.(1-p) ) ;
159	greg	1.3
160			nrand(x) = unif2norm(rand(x));
161
162	greg	1.2	{ Local (u,v) coordinates for planar surfaces }
163	greg	2.17	crosslen`P. = NxNx + NyNy;
164	greg	2.7	{ U is distance from projected Z-axis }
165	greg	2.17	U = if( crosslen`P. - FTINY,
166			(PyNx - PxNy)/crosslen`P.,
167	greg	1.2	Px);
168			{ V is defined so that N = U x V }
169	greg	2.17	V = if( crosslen`P. - FTINY,
170			Pz - Nz(PxNx + Py*Ny)/crosslen`P.,
171	greg	1.2	Py);
172	greg	2.18
173			{ Local hemisphere direction for func & data types }
174			{ last 3 real args = unnormalized up-vector }
175			Vux`P. = arg(AC-1)NzP - arg(AC)NyP;
176			Vuy`P. = arg(AC)NxP - arg(AC-2)NzP;
177			Vuz`P. = arg(AC-2)NyP - arg(AC-1)NxP;
178			vnorm`P. = 1/sqrt(Vux`P.Vux`P. + Vuy`P.Vuy`P. + Vuz`P.*Vuz`P.);
179			Vnx`P. = Vux`P.*vnorm`P.;
180			Vny`P. = Vuy`P.*vnorm`P.;
181			Vnz`P. = Vuz`P.*vnorm`P.;
182			Unx`P. = NyPVnz`P. - NzPVny`P.;
183			Uny`P. = NzPVnx`P. - NxPVnz`P.;
184			Unz`P. = NxPVny`P. - NyPVnx`P.;
185			{ Transform vectors, normalized (dx,dy,dz) away from surf }
186			Ldx(dx,dy,dz) = dxUnx`P. + dyUny`P. + dz*Unz`P.;
187			Ldy(dx,dy,dz) = dxVnx`P. + dyVny`P. + dz*Vnz`P.;
188			Ldz(dx,dy,dz) = dxNxP + dyNyP + dz*NzP;
189			{ Incident vector transformed to our coords }
190			Idx = Ldx(-Dx,-Dy,-Dz);
191			Idy = Ldy(-Dx,-Dy,-Dz);
192			Idz = RdotP;