CUDA Ray-Sphere intersection random walk spooky values

Results

The above results are the |X|Y|Z|AbsDistance of each sphere intersection, random spooky values appear probably because of a newbie mistake, but I really can't get it.

To be as specific as I can:

The following snippet is supposed to calculate the intersection point between a ray and a spherical boundary with a predefined radius and the origin as the center.

To give more context:

1- The RandomWalk starts from the origin and moves with a randomly generated _step and _direction.

2- After each step, the ray is checked for hitting possibility by comparing the absolute distance to the radius of the boundary.

3- getIntersectionPoint() returns the point of intersection, but as the number of points or number of steps increases, the probability of outcasts increases, messing up the whole thing.

Here's what I've done:

main.cu

__global__ void finalPosition(unsigned int seed, curandState_t* states, Point* _gpuPoints,Boundary boundary,RNG rng) {
    int idx = blockIdx.x*blockDim.x+threadIdx.x;
    curand_init(seed, idx, 0, &states[idx]);
    Point finalPos;
    finalPos = randomWalk(states, idx, boundary, rng);
    _gpuPoints[idx] = finalPos;
.
.
.
int main(){
    int nBlocks = N/THREADS_PER_BLOCK + 1;

    curandState_t* states;
    cudaMalloc((void**) &states, N * sizeof(curandState_t));

// Allocate host memory for final positions
    Point * _cpuPoints= (Point*)malloc(sizeof(Point) * N);

// Allocate device  memory for final positions
    Point* _gpuPoints = nullptr;
    cudaMalloc((void**) &_gpuPoints, N * sizeof(Point));

// Initializing the Boundary and the Random Number Generator
    Boundary boundary = Boundary(BOUNDARY_RADIUS, Point(0.f, 0.f, 0.f));
    RNG rng;

// Call Kernel
    finalPosition<<<nBlocks,THREADS_PER_BLOCK>>>(time(0), states , _gpuPoints, boundary, rng);

// Copy device data to host memory to stream them out
    cudaMemcpy(_cpuPoints, _gpuPoints, N* sizeof( Point), cudaMemcpyDeviceToHost);


    streamOut (&_cpuPoints[0]);

    free(_cpuPoints);
    cudaFree(_gpuPoints);


    return 0;
}
}

RandomWalk.h

/**
 * @brief randomWalk
 * keeps wandering around with the photon in the 3D space
 * @return The Point where the Photon hits the Boundary
 */
__device__ Point randomWalk(curandState_t *states, int idx, Boundary boundary, RNG rng)
{
    Ray ray = Ray(Point(0.f, 0.f, 0.f), Point(0.f, 0.f, 0.f));

    while (!boundary.isCrossed(ray))
    {
        ray.move(rng.getRandomPoint(states, idx), rng.getRandomStep(states, idx));
    }
    return boundary.getIntersectionPoint(ray);
}

Boundary.cu

__device__ bool Boundary::isCrossed(Ray ray){
    float absDistance = (float) sqrtf((float) powf(ray.getCurrentPos().getX(),2)
                        + (float) powf(ray.getCurrentPos().getY(),2) 
                        + (float) powf(ray.getCurrentPos().getZ(),2));
    if(absDistance >= _radius){
        return true;
    } else {
        return false;
    }
}


__device__ Point Boundary::getIntersectionPoint(Ray ray){
        /**
            P(t) = A + tB
            P(t) is a point on the ray 
            A is the ray origin
            B is the ray direction
            t is a parameter used to move away from ray origin
            S = P - Center
            ||S||^2 = r^2
            Sphere: dot(S,S) = r^2
            Ray: P(t) = A + tB
            Combined: dot((A + tB - Center),(A + tB - Center)) = r^2
            in Quadratic form: t^2.dot(B,B) + 2t.dot(B, A - C) + dot(A - C, A - C) - r^2 = 0
            let a = dot(B,B)
                b = 2.dot(B, A - C)
                c = dot(A - C, A - C) - r^2
            t1, t2 = (-b (+/-) sqrt(b^2 - 4ac) / 2a)
        */
        Point A = ray.getPrevPos();
        Point B = ray.getDirection();
        Point S = A.add(_center);
        Point A_C = A.subtract(_center);
        float a = dot(B, B);
        float b = 2.0 * dot(B, A_C);
        float c = dot(A_C, A_C) - _radius*_radius;
        float discriminant = b*b - 4*a*c;
        float t1 = (-b + sqrtf(discriminant)) / (2.0*a);
        float t2 = (-b - sqrtf(discriminant)) / (2.0*a);
        float t;

        if(t1 < 0){
            t = t2;
        } else {
            t = t1;
        }

        return Point((A.getX()+B.getX()*t),(A.getY()+B.getY()*t),(A.getZ()+B.getZ()*t));
}

RNG.cu

__device__  float RNG::generate( curandState* globalState, int i) 
{
    curandState localState = globalState[i];
    float random = curand_uniform( &localState );
    globalState[i] = localState;
    return random;
}

__device__   float RNG::getRandomStep( curandState* globalState , int i) { 
    float step = 0.f;       // Intialize for step value
    step = generate (globalState, i);
    return step;
 } 

__device__  Point RNG::getRandomPoint( curandState* globalState , int i)
{

    float u = generate (globalState , i);
    float v = generate (globalState, i);

    float theta = 2 * M_PI * u;
    float phi = acos(1 - 2 * v);

    // Transforming into the cartesian space
    float x = sin(phi) * cos(theta);
    float y = sin(phi) * sin(theta);
    float z = cos(phi);

    return Point(x,y,z);
}