Can't find an analytical sphere (using tutorial)

[AdriMerle]

(I'm french so I apologize in advance if I make any english mistake)

Hello everyone,

I am having trouble creating a User Geometry in Embree 4. I need to create a set of them (sphere, cylinder...) but let's focus on the sphere for now.

For some reason I cannot find an intersection with the sphere I created. I basically copy pasted this tutorial (just modified a few things so that it works in my setup).

I know the code is quite long and I am sorry about it.

void sphereBoundsFunc(const struct RTCBoundsFunctionArguments* args)
{
    const Sphere* spheres = (const Sphere*)args->geometryUserPtr;
    RTCBounds* bounds = args->bounds_o;
    const Sphere& sphere = spheres[args->primID];
    bounds->lower_x = sphere.p.x - sphere.r;
    bounds->lower_y = sphere.p.y - sphere.r;
    bounds->lower_z = sphere.p.z - sphere.r;
    bounds->upper_x = sphere.p.x + sphere.r;
    bounds->upper_y = sphere.p.y + sphere.r;
    bounds->upper_z = sphere.p.z + sphere.r;
}

RTC_SYCL_INDIRECTLY_CALLABLE void sphereIntersectFunc(const RTCIntersectFunctionNArguments* args)
{
    int* valid = args->valid;
    void* ptr = args->geometryUserPtr;
    RTCRayHit *rayhit = (RTCRayHit*)args->rayhit;
    // RTCHit* hit = (RTCHit*)&ray->Ng.x; // ??????????
    unsigned int primID = args->primID;

    assert(args->N == 1);
    const Sphere* spheres = (const Sphere*)ptr;
    const Sphere& sphere = spheres[primID];

    // valid=0  -> le rayon est invalide
    // valid=-1 -> le rayon est valide
    if (!valid[0]) return;
    valid[0] = 0;

    // A*X^2 + B*X + C
    const Point dir = (rayhit->ray.dir_x, rayhit->ray.dir_y, rayhit->ray.dir_z);
    const Point org = (rayhit->ray.org_x, rayhit->ray.org_y, rayhit->ray.org_z);
    const Point v = org - sphere.p;
    const float A = dir.dot(dir); 
    const float B = 2.0f * v.dot(dir);
    const float C = v.dot(v) - sphere.r * sphere.r;
    const float D = B * B - 4.0f * A * C; //delta = b^2 -4*a*c
    if (D < 0.0f) return;
    const float Q = sqrt(D);
    const float t0 = 0.5f * (-B - Q) / A;
    const float t1 = 0.5f * (-B + Q) / A;

    RTCHit potentialHit;
    potentialHit.u = 0.0f;
    potentialHit.v = 0.0f;
    copyInstanceIdStack(args->context, potentialHit.instID);
    potentialHit.geomID = sphere.geomID;
    potentialHit.primID = primID;

    if ((rayhit->ray.tnear < t0) && (t0 < rayhit->ray.tfar))
    {
        int imask;
        bool mask = 1;
        {
            imask = mask ? -1 : 0;
        }

        const Point Ng = org + dir * t0 - sphere.p;
        potentialHit.Ng_x = Ng.x;
        potentialHit.Ng_y = Ng.y;
        potentialHit.Ng_z = Ng.z;

        RTCFilterFunctionNArguments fargs;
        fargs.valid = (int*)&imask;
        fargs.geometryUserPtr = ptr;
        fargs.context = args->context;
        fargs.ray = (RTCRayN*)args->rayhit;
        fargs.hit = (RTCHitN*)&potentialHit;
        fargs.N = 1;

        const float old_t = rayhit->ray.tfar;
        rayhit->ray.tfar = t0;
#if USE_ARGUMENT_CALLBACKS
        contextFilterFunction(&fargs);
#else
        rtcInvokeIntersectFilterFromGeometry(args, &fargs);
#endif

        if (imask == -1) {
            rayhit->hit = potentialHit;
            valid[0] = -1;
        }
        else
            rayhit->ray.tfar = old_t;
    }

    if ((rayhit->ray.tnear < t1) && (t1 < rayhit->ray.tfar))
    {
        int imask;
        bool mask = 1;
        {
            imask = mask ? -1 : 0;
        }

        const Point Ng = org + dir * t1 - sphere.p;
        potentialHit.Ng_x = Ng.x;
        potentialHit.Ng_y = Ng.y;
        potentialHit.Ng_z = Ng.z;

        RTCFilterFunctionNArguments fargs;
        fargs.valid = (int*)&imask;
        fargs.geometryUserPtr = ptr;
        fargs.context = args->context;
        fargs.ray = (RTCRayN*)args->rayhit;
        fargs.hit = (RTCHitN*)&potentialHit;
        fargs.N = 1;

        const float old_t = rayhit->ray.tfar;
        rayhit->ray.tfar = t1;
#if USE_ARGUMENT_CALLBACKS
        contextFilterFunction(&fargs);
#else
        rtcInvokeIntersectFilterFromGeometry(args, &fargs);
#endif

        if (imask == -1) {
            rayhit->hit = potentialHit;
            valid[0] = -1;
        }
        else
            rayhit->ray.tfar = old_t;
    }
}

RTC_SYCL_INDIRECTLY_CALLABLE void sphereOccludedFunc(const RTCOccludedFunctionNArguments* args)
{
    int* valid = args->valid;
    void* ptr = args->geometryUserPtr;
    RTCRayHit* rayhit = (RTCRayHit*)args->ray;
    unsigned int primID = args->primID;

    assert(args->N == 1);
    const Sphere* spheres = (const Sphere*)ptr;
    const Sphere& sphere = spheres[primID];

    if (!valid[0]) return;
    valid[0] = 0;

    // A*X^2 + B*X + C
    const Point dir = (rayhit->ray.dir_x, rayhit->ray.dir_y, rayhit->ray.dir_z);
    const Point org = (rayhit->ray.org_x, rayhit->ray.org_y, rayhit->ray.org_z);
    const Point v = org - sphere.p;
    const float A = dir.dot(dir);
    const float B = 2.0f * v.dot(dir);
    const float C = v.dot(v) - sphere.r * sphere.r;
    const float D = B * B - 4.0f * A * C; //delta
    if (D < 0.0f) return;
    const float Q = sqrt(D);
    const float t0 = 0.5f * (-B - Q) / A;
    const float t1 = 0.5f * (-B + Q) / A;

    RTCHit potentialHit;
    potentialHit.u = 0.0f;
    potentialHit.v = 0.0f;
    copyInstanceIdStack(args->context, potentialHit.instID);
    potentialHit.geomID = sphere.geomID;
    potentialHit.primID = primID;

    if ((rayhit->ray.tnear < t0) && (t0 < rayhit->ray.tfar))
    {
        int imask;
        bool mask = 1;
        {
            imask = mask ? -1 : 0;
        }

        const Point Ng = org + dir * t0 - sphere.p;
        potentialHit.Ng_x = Ng.x;
        potentialHit.Ng_y = Ng.y;
        potentialHit.Ng_z = Ng.z;

        RTCFilterFunctionNArguments fargs;
        fargs.valid = (int*)&imask;
        fargs.geometryUserPtr = ptr;
        fargs.context = args->context;
        fargs.ray = (RTCRayN*)args->ray;
        fargs.hit = (RTCHitN*)&potentialHit;
        fargs.N = 1;

        const float old_t = rayhit->ray.tfar;
        rayhit->ray.tfar = t0;
#if USE_ARGUMENT_CALLBACKS
        contextFilterFunction(&fargs);
#else
        rtcInvokeOccludedFilterFromGeometry(args, &fargs);
#endif

        if (imask == -1) {
            rayhit->hit = potentialHit;
            valid[0] = -1;
        }
        else
            rayhit->ray.tfar = old_t;
    }

    if ((rayhit->ray.tnear < t1) && (t1 < rayhit->ray.tfar))
    {
        int imask;
        bool mask = 1;
        {
            imask = mask ? -1 : 0;
        }

        const Point Ng = org + dir * t1 - sphere.p;
        potentialHit.Ng_x = Ng.x;
        potentialHit.Ng_y = Ng.y;
        potentialHit.Ng_z = Ng.z;

        RTCFilterFunctionNArguments fargs;
        fargs.valid = (int*)&imask;
        fargs.geometryUserPtr = ptr;
        fargs.context = args->context;
        fargs.ray = (RTCRayN*)args->ray;
        fargs.hit = (RTCHitN*)&potentialHit;
        fargs.N = 1;

        const float old_t = rayhit->ray.tfar;
        rayhit->ray.tfar = t1;
#if USE_ARGUMENT_CALLBACKS
        contextFilterFunction(&fargs);
#else
        rtcInvokeOccludedFilterFromGeometry(args, &fargs);
#endif

        if (imask == -1) {
            rayhit->hit = potentialHit;
            valid[0] = -1;
        }
        else
            rayhit->ray.tfar = old_t;
    }
}

RTC_SYCL_INDIRECTLY_CALLABLE void contextIntersectFunc(const RTCIntersectFunctionNArguments* args)
{
    UserGeometryType* type = (UserGeometryType*)args->geometryUserPtr;
    sphereIntersectFunc(args);
    //if (*type == USER_GEOMETRY_SPHERE) sphereIntersectFunc(args);
    // Lister ici les types de géométries pour appeler la bonne fonction à chaque fois
}

RTC_SYCL_INDIRECTLY_CALLABLE void sphereFilterFunction(const RTCFilterFunctionNArguments* args)
{
    int* valid = args->valid;
    const RayQueryContext* context = (const RayQueryContext*)args->context;
    RTCRay* ray = (RTCRay*)args->ray;
    RTCHit* hit = (RTCHit*)args->hit;
    const unsigned int N = args->N;
    assert(N == 1);
    // _unused(N); // ???????

    /* avoid crashing when debug visualizations are used */
    if (context == nullptr)
        return;

    /* ignore inactive rays */
    if (valid[0] != -1) return;

    /* carve out parts of the sphere */
    const Point h = (ray->org_x + ray->dir_x * ray->tfar,
        ray->org_y + ray->dir_y * ray->tfar,
        ray->org_z + ray->dir_z * ray->tfar);
    float v = abs(sin(10.0f * h.x) * cos(10.0f * h.y) * sin(10.0f * h.z));
    float T = clamp((v - 0.1f) * 3.0f, 0.0f, 1.0f);

    /* reject some hits */
    if (T < 0.5f) valid[0] = 0;
}

RTC_SYCL_INDIRECTLY_CALLABLE void contextFilterFunction(const RTCFilterFunctionNArguments* args)
{
    int* valid = args->valid;
    if (!valid[0]) return;

    RTCHit* potential_hit = (RTCHit*)args->hit;
    if (potential_hit->instID[0] == 0)
        sphereFilterFunction(args);
}

RTCGeometry UDG::analyticalSphere(RTCDevice device, const Point o, float r)
{
    RTCGeometry geom = rtcNewGeometry(device, RTC_GEOMETRY_TYPE_USER);
    Sphere sphere = Sphere();
    sphere.type = USER_GEOMETRY_SPHERE;
    sphere.p = o;
    sphere.r = r;
    sphere.geometry = geom;
    //sphere.geomID = rtcAttachGeometry(scene, geom);
    rtcSetGeometryUserPrimitiveCount(geom, 1);
    rtcSetGeometryUserData(geom, &sphere);
    rtcSetGeometryBoundsFunction(geom, sphereBoundsFunc, nullptr);
#if !USE_ARGUMENT_CALLBACKS
    rtcSetGeometryIntersectFunction(geom, sphereIntersectFunc);
    rtcSetGeometryOccludedFunction(geom, sphereOccludedFunc);
#endif
    rtcCommitGeometry(geom);
    //rtcReleaseGeometry(geom);
    return geom;
}

I tried the code above using this setup :

    Point center = (0.f, 0.f, 0.f);
    RTCGeometry geom = UDG::analyticalSphere(api.device, center, 2);
    string scene = "UDG";

    api.newScene(scene);
    unsigned int geomID = api.attachGeometry(geom, scene); // custom rtcAttachGeometry that works (tested with triangles etc.)

    Point o;    o.x = 0.0f;    o.y = 0.0f;     o.z = -10.0f;
    Point d;    d.x = 0;       d.y = 0;        d.z = 1;
    RTCRayHit ray = api.newRay(o, d); // Creates a RayHit of origin o and direction d

    api.getIntersection(ray, scene);

And embree doesn't update the tfar component (hence no intersection).

Please let me know if you find any mistake in that code.

Thanks in advance.

Adrian

[svenwoop]

Does your attachGeometry call rtcCommit on the scene? I would just start with a single triangle to get Embree integration working, and then go to more complicated primitives. BTW, there is also direct support for spheres using the RTC_GEOMETRY_TYPE_SPHERE_POINT geometry type.

[AdriMerle]

Thanks for your reply.

Yes it does. Here is its code :

unsigned int LibEmbree::attachGeometry(RTCGeometry geom, const string sceneName)
{
    rtcCommitGeometry(geom);
    unsigned int geomID = rtcAttachGeometry(mapScenes[sceneName], geom);
    rtcCommitScene(mapScenes[sceneName]);
    return geomID;
}

Yes, I found that, but I will need to make cylinders and cones in the future and I believe that they are not supported natively. That's why I chose to learn how to do this with analytical method.

[svenwoop]

That code looks correct. The mapScenes is likely an std:map and looking up scenes there will significantly reduce performance.
Anyway, I would start with getting triangles working first. Also note that the sphere intersection you copied has some filter function attached, which will reject some hits (that is its purpose).

[AdriMerle]

Triangles are definitely working in my setup. I know that the map is reducing performances, yet it is simpler for me (for the time that I am constructing my program) to keep it.

I narrowed the error to these lines :

RTC_SYCL_INDIRECTLY_CALLABLE void sphereIntersectFunc(const RTCIntersectFunctionNArguments* args)
{
    cout << __FUNCTION__ << endl;
    int* valid = args->valid;
    void* ptr = args->geometryUserPtr;
    RTCRayHit *rayhit = (RTCRayHit*)args->rayhit;
    // RTCHit* hit = (RTCHit*)&ray->Ng.x; // ??????????
    unsigned int primID = args->primID;

    assert(args->N == 1);
    const Sphere* spheres = (const Sphere*)ptr; // PROBLEME ICI
    const Sphere& sphere = spheres[primID];

The ptr seems to be wrong, since the sphere is completely wrong here (hence messing with my computations).

Any idea ?

Thanks for your answer !

[svenwoop]

This line here "rtcSetGeometryUserData(geom, &sphere);" sets a pointer to a local variable. You have to allocate the sphere object on the heap.

[AdriMerle]

Oops you're right, I already used an allocation trying to debug (not in the first code I sent though) but forgot to remove the & before sphere. Hence, the pointer was numb.

Well, I believe that my question is now answered.

Thanks a lot for your help ! Have a good day.

[AdriMerle]

Hi, I'm back again :)

These two functions together work well with triangular shapes (tested with cubes and triangles) but doesn't seem to work with my analytical sphere.

What should be added/modified to these functions in order to gather all hits ? (whether they are from analytical or tessellated primitives, I don't care, I only want to be able to sort them later by geomID).

Filter :

RTC_SYCL_INDIRECTLY_CALLABLE void gather_all_hits(const RTCFilterFunctionNArguments* args)
{
    assert(*args->valid == -1);
    RayQueryContext* context = (RayQueryContext*)args->context;
    HitList& hits = context->hits;
    RTCRay* ray = (RTCRay*)args->ray;
    RTCHit* hit = (RTCHit*)args->hit;
    assert(args->N == 1);
    args->valid[0] = 0; // ignore all hits

    /* avoid overflow of hits array */
    if (hits.end >= MAX_TOTAL_HITS) return;

    hits.hits[hits.end++] = HitList::Hit(false, ray->tfar, hit->primID, hit->geomID, hit->instID[0]);
}

Function that gathers hits :

HitList* LibEmbree::getHits(RTCRayHit rayhit, const string sceneName)
{
    HitList * hits_o = new HitList();
    /* trace ray to gather all hits */
    RayQueryContext context(MAX_TOTAL_HITS, *hits_o);
    rtcInitRayQueryContext(&context.context);

    RTCIntersectArguments args;
    rtcInitIntersectArguments(&args);
    args.context = &context.context;
    args.filter = gather_all_hits;
    args.flags = RTC_RAY_QUERY_FLAG_INVOKE_ARGUMENT_FILTER; // invoke filter for each geometry
    rtcIntersect1(mapScenes[sceneName], &rayhit, &args);

    /* sort hits by extended order */
    //std::sort(&context.hits.hits[context.hits.begin], &context.hits.hits[context.hits.end]);

    /* drop hits in case we found too many */
    hits_o->end = std::min(hits_o->end, (unsigned int)MAX_TOTAL_HITS);

    // Pour debug : Afficher les membres du vecteur
    //for (auto d : distances)
    //    cout << d << " ";    hits_o = HitList();
    return hits_o;
}

Thanks again !

[svenwoop]

Make sure the ray.tfar is keeping old value, thus this line should get executed in sphere intersector: rayhit->ray.tfar = old_t;

[AdriMerle]

Yes it is, at the end of the ifs (that check which solution should be kept between t0 and t1), we put back old_t if the test with mask fails.

I have absolutely no idea how to pursue to make the "multi-hit gathering" with an alaytical sphere. Thanks again for your help.

Here is the last version of my code :

RTC_SYCL_INDIRECTLY_CALLABLE void sphereIntersectFunc(const RTCIntersectFunctionNArguments* args)
{
#ifdef _DEBUG
    cout << __FUNCTION__ << endl;
#endif // _DEBUG
    int* valid = args->valid;
    void* ptr = args->geometryUserPtr;
    RTCRayHit *rayhit = (RTCRayHit*)args->rayhit;
    // RTCHit* hit = (RTCHit*)&ray->Ng.x; // ??????????
    unsigned int primID = args->primID;

    assert(args->N == 1);
    const Sphere* spheres = (const Sphere*)ptr; // PROBLEME ICI
    const Sphere& sphere = spheres[primID];

    // valid=0  -> le rayon est invalide
    // valid=-1 -> le rayon est valide
    if (!valid[0]) return;
    valid[0] = 0;

    // A*X^2 + B*X + C
    const Point dir = Point(rayhit->ray.dir_x, rayhit->ray.dir_y, rayhit->ray.dir_z);
    const Point org = Point(rayhit->ray.org_x, rayhit->ray.org_y, rayhit->ray.org_z);
    const Point v = org - sphere.p;
    const float A = dir.dot(dir); 
    const float B = 2.0f * v.dot(dir);
    const float C = v.dot(v) - sphere.r * sphere.r;
    const float D = B * B - 4.0f * A * C; //delta = b^2 -4*a*c
    if (D < 0.0f) return;
    const float Q = sqrt(D);
    const float t0 = 0.5f * (-B - Q) / A;
    const float t1 = 0.5f * (-B + Q) / A;

    RTCHit potentialHit;
    potentialHit.u = 0.0f;
    potentialHit.v = 0.0f;
    copyInstanceIdStack(args->context, potentialHit.instID);
    potentialHit.geomID = sphere.geomID;
    potentialHit.primID = primID;

    if ((rayhit->ray.tnear < t0) && (t0 < rayhit->ray.tfar))
    {
        int imask;
        bool mask = 1;
        {
            imask = mask ? -1 : 0;
        }

        const Point Ng = org + dir * t0 - sphere.p;
        potentialHit.Ng_x = Ng.x;
        potentialHit.Ng_y = Ng.y;
        potentialHit.Ng_z = Ng.z;

        RTCFilterFunctionNArguments fargs;
        fargs.valid = (int*)&imask;
        fargs.geometryUserPtr = ptr;
        fargs.context = args->context;
        fargs.ray = (RTCRayN*)args->rayhit;
        fargs.hit = (RTCHitN*)&potentialHit;
        fargs.N = 1;

        const float old_t = rayhit->ray.tfar;
        rayhit->ray.tfar = t0;
#if USE_ARGUMENT_CALLBACKS
        contextFilterFunction(&fargs);
#else
        rtcInvokeIntersectFilterFromGeometry(args, &fargs);
#endif

        if (imask == -1) {
            rayhit->hit = potentialHit;
            valid[0] = -1;
        }
        else
            rayhit->ray.tfar = old_t;
    }

    if ((rayhit->ray.tnear < t1) && (t1 < rayhit->ray.tfar))
    {
        int imask;
        bool mask = true;
        {
            imask = mask ? -1 : 0;
        }

        const Point Ng = org + dir * t1 - sphere.p;
        potentialHit.Ng_x = Ng.x;
        potentialHit.Ng_y = Ng.y;
        potentialHit.Ng_z = Ng.z;

        RTCFilterFunctionNArguments fargs;
        fargs.valid = (int*)&imask;
        fargs.geometryUserPtr = ptr;
        fargs.context = args->context;
        fargs.ray = (RTCRayN*)args->rayhit;
        fargs.hit = (RTCHitN*)&potentialHit;
        fargs.N = 1;

        const float old_t = rayhit->ray.tfar;
        rayhit->ray.tfar = t1;
#if USE_ARGUMENT_CALLBACKS
        contextFilterFunction(&fargs);
#else
        rtcInvokeIntersectFilterFromGeometry(args, &fargs);
#endif

        if (imask == -1) {
            rayhit->hit = potentialHit;
            valid[0] = -1;
        }
        else
            rayhit->ray.tfar = old_t;
    }
}

RTC_SYCL_INDIRECTLY_CALLABLE void sphereFilterFunc(const RTCFilterFunctionNArguments* args)
{
#ifdef _DEBUG
    cout << __FUNCTION__ << endl;
#endif // _DEBUG
    int* valid = args->valid;
    const RayQueryContext* context = (const RayQueryContext*)args->context;
    RTCRay* ray = (RTCRay*)args->ray;
    RTCHit* hit = (RTCHit*)args->hit;
    const unsigned int N = args->N;
    assert(N == 1);
    // _unused(N); // ???????

    /* avoid crashing when debug visualizations are used */
    if (context == nullptr)
        return;

    /* ignore inactive rays */
    if (valid[0] != -1) return;

    /* carve out parts of the sphere */
    const Point h = (ray->org_x + ray->dir_x * ray->tfar,
        ray->org_y + ray->dir_y * ray->tfar,
        ray->org_z + ray->dir_z * ray->tfar);
    float v = abs(sin(10.0f * h.x) * cos(10.0f * h.y) * sin(10.0f * h.z));
    float T = clamp((v - 0.1f) * 3.0f, 0.0f, 1.0f);

    /* reject some hits */
    if (T < 0.5f) valid[0] = 0;
}

Sphere* UDG::analytical::sphere(RTCDevice device, const Point o, float r)
{
#ifdef _DEBUG
    cout << __FUNCTION__ << endl;
#endif // _DEBUG
    RTCGeometry geom = rtcNewGeometry(device, RTC_GEOMETRY_TYPE_USER);
    Sphere * sphere = new Sphere();
    sphere->type = USER_GEOMETRY_SPHERE;
    sphere->p = o;
    sphere->r = r;
    sphere->geometry = geom;
    //sphere.geomID = rtcAttachGeometry(scene, geom);
    rtcSetGeometryUserPrimitiveCount(geom, 1);
    rtcSetGeometryUserData(geom, sphere);
    rtcSetGeometryBoundsFunction(geom, sphereBoundsFunc, nullptr);
    rtcSetGeometryIntersectFunction(geom, sphereIntersectFunc);
    rtcSetGeometryIntersectFilterFunction(geom, sphereFilterFunc);
    rtcCommitGeometry(geom);
    //rtcReleaseGeometry(geom);
    return sphere;
}

And I have a question :
Is the intersect function called only one time because the ray enters only one time the bounding box ?
If yes, does it mean that I have to make manually the multihit using the sphere filter function (once I have one intersection I tell the filter to compute the other point) ?

Thanks in advance.

[svenwoop]

Is the gather_all_hits callback invoked at all? The rtcInvokeIntersectFilterFromGeometry call just invokes the filter function assigned to the geometry, not the one in the context. The context callback has to get invoked manually for user defined geometries, thus instead of rtcInvokeIntersectFilterFromGeometry inoke your gather_all_hits directly (either by just calling that function directly or by also adding the RTCIntersectArguments to your RayQueryContext and invoking it through the function pointer).

[AdriMerle]

Thanks a lot for your help, I guess that it's all working now