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