How to create a User Geometry in Embree 4? (example of a sphere)

Question

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.