Cycles: svn merge -r41225:41232 ^/trunk/blender
[blender.git] / intern / cycles / kernel / kernel_differential.h
1 /*
2  * Copyright 2011, Blender Foundation.
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18
19 CCL_NAMESPACE_BEGIN
20
21 /* See "Tracing Ray Differentials", Homan Igehy, 1999. */
22
23 __device void differential_transfer(differential3 *dP_, const differential3 dP, float3 D, const differential3 dD, float3 Ng, float t)
24 {
25         /* ray differential transfer through homogenous medium, to
26          * compute dPdx/dy at a shading point from the incoming ray */
27
28         float3 tmp = D/dot(D, Ng);
29         float3 tmpx = dP.dx + t*dD.dx;
30         float3 tmpy = dP.dy + t*dD.dy;
31
32         dP_->dx = tmpx - dot(tmpx, Ng)*tmp;
33         dP_->dy = tmpy - dot(tmpy, Ng)*tmp;
34 }
35
36 __device void differential_incoming(differential3 *dI, const differential3 dD)
37 {
38         /* compute dIdx/dy at a shading point, we just need to negate the
39          * differential of the ray direction */
40
41         dI->dx = -dD.dx;
42         dI->dy = -dD.dy;
43 }
44
45 __device void differential_dudv(differential *du, differential *dv, float3 dPdu, float3 dPdv, differential3 dP, float3 Ng)
46 {
47         /* now we have dPdx/dy from the ray differential transfer, and dPdu/dv
48          * from the primitive, we can compute dudx/dy and dvdx/dy. these are
49          * mainly used for differentials of arbitrary mesh attributes. */
50
51         /* find most stable axis to project to 2D */
52         float xn= fabsf(Ng.x);
53         float yn= fabsf(Ng.y);
54         float zn= fabsf(Ng.z);
55
56         if(zn < xn || zn < yn) {
57                 if(yn < xn || yn < zn) {
58                         dPdu.x = dPdu.y;
59                         dPdv.x = dPdv.y;
60                         dP.dx.x = dP.dx.y;
61                         dP.dy.x = dP.dy.y;
62                 }
63
64                 dPdu.y = dPdu.z;
65                 dPdv.y = dPdv.z;
66                 dP.dx.y = dP.dx.z;
67                 dP.dy.y = dP.dy.z;
68         }
69
70         /* using Cramer's rule, we solve for dudx and dvdx in a 2x2 linear system,
71          * and the same for dudy and dvdy. the denominator is the same for both
72          * solutions, so we compute it only once.
73          *
74      * dP.dx = dPdu * dudx + dPdv * dvdx;
75      * dP.dy = dPdu * dudy + dPdv * dvdy; */
76
77         float det = (dPdu.x*dPdv.y - dPdv.x*dPdu.y);
78
79         if(det != 0.0f)
80                 det = 1.0f/det;
81
82         du->dx = (dP.dx.x*dPdv.y - dP.dx.y*dPdv.x)*det;
83         dv->dx = (dP.dx.y*dPdu.x - dP.dx.x*dPdu.y)*det;
84
85         du->dy = (dP.dy.x*dPdv.y - dP.dy.y*dPdv.x)*det;
86         dv->dy = (dP.dy.y*dPdu.x - dP.dy.x*dPdu.y)*det;
87 }
88
89 CCL_NAMESPACE_END
90