rayCast() in cuda_gvdb_raycast.cuh may need a condition check when DDA.mask is populated with zero values.

[digbeta]

I've been chasing a weird bug for some time that I think I may have tracked down into the master raycast function in cuda_gvdb_raycast.cuh.

I have a fix, although it's not as elegant as what Rama and/or Neil may identify, so I'll explain the issue and my fix.

I've noticed that when rays are fired that have zero values in the direction vector (may also be a requirement that the rays hit the box boundary itself), the DDA state will be populated with nan and/or inf values that appear (I believe) when the bounding box intersection tests are performed in rayBoxIntersect(). The results are cascaded into the mask value for the DDA function which will end up being all zeroes:

    // Compute the next time step and direction from DDA, but don't step yet.
    __device__ void Next() {
        mask.x = int((tSide.x < tSide.y) & (tSide.x <= tSide.z));
        mask.y = int((tSide.y < tSide.z) & (tSide.y <= tSide.x));
        mask.z = int((tSide.z < tSide.x) & (tSide.z <= tSide.y));
        t.y = mask.x ? tSide.x : (mask.y ? tSide.y : tSide.z);
    }

When this happens, dda.Step() is basically a NOOP, because it doesn't actually step:

    // Step the DDA to the next point.
    __device__ void Step() {
        t.x = t.y;
        tSide += make_float3(mask) * tDel;
        p += mask * pStep;
    }

The result is the for loop basically loops, doing no real work, until the MAX_ITER check is hit and breaks the loop. This may have been the purpose of MAX_ITER, I am not sure, but instead I added a check after dda.Next() which simply breaks the loop if the mask is zeroed out.

        if (dda.mask.x == 0 && dda.mask.y == 0 && dda.mask.z == 0)
        {
            break;
        }

There may be other consequences to this that I am missing, but this was causing problems for me when I set MAX_ITER very high or removed it altogether. If I am missing something else here, please let me know! Thanks in advance.

[NBickford-NV]

That makes sense! I think the NaNs and Infs are coming from the division by dir in dda.Prepare():

    // Set up variables for a node to prepare for stepping.
    __device__ void Prepare(float3 vmin, float3 vdel) {
        tDel = fabs3(vdel / dir);
        float3 pFlt = (pos + t.x * dir - vmin) / vdel;
        tSide = ((floor3(pFlt) - pFlt + 0.5f) * make_float3(pStep) + 0.5) * tDel + t.x;
        p = make_int3(floor3(pFlt));
    }

[digbeta]

Right... I didn't dive into the DDA algorithm, which is why I went with the very basic value check so I didn't have to touch any of your code. I wasn't sure whether this code took into account those boundary cases or not... :-)