Skip to content
Permalink
Browse files

Fix dual implementation of division so that the value is consistent w…

…ith regular division (#1066)
  • Loading branch information...
fpsunflower committed Oct 3, 2019
1 parent 55bd37b commit 02d364e7cbf05409d23ef9c880be318184fc5b2f
@@ -386,8 +386,8 @@ operator* (const S &b, const Dual<T,P> &a) -> Dual<decltype(a.elem(0)*b),P>
template<class T, int P>
OSL_HOSTDEVICE inline OIIO_CONSTEXPR14 Dual<T,P> operator/ (const Dual<T,P> &a, const Dual<T,P> &b)
{
T bvalinv = 1.0f / b.val();
T aval_bval = a.val() * bvalinv;
T bvalinv = T(1) / b.val();
T aval_bval = a.val() / b.val();
Dual<T,P> result;
result.val() = aval_bval;
for (int i = 1; i <= P; ++i)
@@ -401,8 +401,13 @@ OSL_HOSTDEVICE inline OIIO_CONSTEXPR14 Dual<T,P> operator/ (const Dual<T,P> &a,
template<class T, int P>
OSL_HOSTDEVICE inline OIIO_CONSTEXPR14 Dual<T,P> operator/ (const Dual<T,P> &a, const T &b)
{
T binv = 1.0f / b;
return a * binv;
T bvalinv = T(1) / b;
T aval_bval = a.val() / b;
Dual<T,P> result;
result.val() = aval_bval;
for (int i = 1; i <= P; ++i)
result.elem(i) = bvalinv * a.elem(i);
return result;
}


@@ -411,8 +416,8 @@ OSL_HOSTDEVICE inline OIIO_CONSTEXPR14 Dual<T,P> operator/ (const Dual<T,P> &a,
template<class T, int P>
OSL_HOSTDEVICE inline OIIO_CONSTEXPR14 Dual<T,P> operator/ (const T &aval, const Dual<T,P> &b)
{
T bvalinv = 1.0f / b.val();
T aval_bval = aval * bvalinv;
T bvalinv = T(1) / b.val();
T aval_bval = aval / b.val();
Dual<T,P> result;
result.val() = aval_bval;
for (int i = 1; i <= P; ++i)
@@ -446,15 +446,17 @@ template<class T, int P>
OSL_HOSTDEVICE inline Dual<Imath::Vec3<T>,P>
normalize (const Dual<Imath::Vec3<T>,P> &a)
{
// NOTE: math must be consistent with osl_normalize_vv
// TODO: math for derivative elements could be further optimized ...
auto ax = comp (a, 0);
auto ay = comp (a, 1);
auto az = comp (a, 2);
auto length = sqrt(ax * ax + ay * ay + az * az);
if (length > 0.0f) {
// NOTE: do a full division here to match what OpenEXR Imath does in the non-dual case
ax = ax / length;
ay = ay / length;
az = az / length;
auto len = sqrt(ax * ax + ay * ay + az * az);
if (len > T(0)) {
auto invlen = T(1) / len;
ax = ax * invlen;
ay = ay * invlen;
az = az * invlen;
return make_Vec3 (ax, ay, az);
} else {
return Vec3(0,0,0);
@@ -715,11 +715,21 @@ osl_distance_dfvdv (void *result, void *a, void *b)
DFLOAT(result) = distance (VEC(a), DVEC(b));
}


OSL_SHADEOP void
osl_normalize_vv (void *result, void *a)
{
VEC(result) = VEC(a).normalized();
using std::sqrt;
// NOTE: must match with the Dual version of normalize used below
Vec3 v = VEC(a);
float len = sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
if (len > 0) {
float invlen = 1 / len;
v.x *= invlen;
v.y *= invlen;
v.z *= invlen;
} else
v.x = v.y = v.z = 0;
VEC(result) = v;
}

OSL_SHADEOP void
@@ -202,7 +202,7 @@ cross product: c = cross(1 -1 1, 1 1 1) = -2 0 2 Dx(c) = -4 0 4 Dy(c) = -4 0 4

comp assign: C[0]=u, C[1]=v: now C = 1 0 0, Dx(C) = 1 0 0, Dy(C) = 0 1 0

normalize: n = normalize(1 1 1) = 0.5774 0.5774 0.5774 Dx(n) = 3.441e-08 3.441e-08 3.441e-08 Dy(n) = 0.3849 -0.7698 0.3849
normalize: n = normalize(1 1 1) = 0.5774 0.5774 0.5774 Dx(n) = 5.96e-08 5.96e-08 5.96e-08 Dy(n) = 0.3849 -0.7698 0.3849

length: l = length(1 1 1) = 1.732 Dx(l) = 1.732 Dy(l) = 0.5774

@@ -309,7 +309,7 @@ cross product: c = cross(1 -1 1, 1 1 1) = -2 0 2 Dx(c) = -4 0 4 Dy(c) = -4 0 4

comp assign: C[0]=u, C[1]=v: now C = 0 1 0, Dx(C) = 1 0 0, Dy(C) = 0 1 0

normalize: n = normalize(1 -1 1) = 0.5774 -0.5774 0.5774 Dx(n) = 0.3849 0.7698 0.3849 Dy(n) = 3.441e-08 -3.441e-08 3.441e-08
normalize: n = normalize(1 -1 1) = 0.5774 -0.5774 0.5774 Dx(n) = 0.3849 0.7698 0.3849 Dy(n) = 5.96e-08 -5.96e-08 5.96e-08

length: l = length(1 -1 1) = 1.732 Dx(l) = 0.5774 Dy(l) = 1.732

@@ -416,7 +416,7 @@ cross product: c = cross(2 -2 2, 2 2 2) = -8 0 8 Dx(c) = -8 0 8 Dy(c) = -8 0 8

comp assign: C[0]=u, C[1]=v: now C = 1 1 0, Dx(C) = 1 0 0, Dy(C) = 0 1 0

normalize: n = normalize(2 0 2) = 0.7071 0 0.7071 Dx(n) = 2.107e-08 0.3536 2.107e-08 Dy(n) = 2.107e-08 -0.3536 2.107e-08
normalize: n = normalize(2 0 2) = 0.7071 0 0.7071 Dx(n) = 2.98e-08 0.3536 2.98e-08 Dy(n) = 2.98e-08 -0.3536 2.98e-08

length: l = length(2 0 2) = 2.828 Dx(l) = 1.414 Dy(l) = 1.414

Binary file not shown.
Binary file not shown.
Binary file not shown.

0 comments on commit 02d364e

Please sign in to comment.
You can’t perform that action at this time.