Issue raycasting when direction vector nears zero in x and z axis

Question

I am working on a voxel game right now in OpenGL and working on implementing raycasting so I can pick blocks in my world. So far what I have created works quite well except for when the x or z axis that is created from my rayWOR vector nears 0 or ±1. I have tried to investigate this further by doing ceil when below zero and floor otherwise but still problematic. I attached what I have created - it is modified from this for C++.

I attempted to introduce an epsilon variable for calculating the d and max variables but that did not fix my problem and I noticed no changes. If anyone happens to notice something I am not doing - or shouldn't be doing, it would be greatly appreciated.

glm::vec3 chunkManager::mouseVoxel(Raycast &ray, Camera &camera) {
    using namespace glm;

    float distance = 10.0f;
    float cellSize = 1.0f;

    vec3 rayWOR = ray.getCurrentRay(); //normalized direction vector from ray (at enter of screen)
    vec3 startPoint = camera.getCameraWorldPosition();

    vec3 endPoint = startPoint + (rayWOR * distance);

    vec3 startCell (startPoint.x < 0.0f ? ceil(startPoint.x) : floor(startPoint.x), startPoint.y < 0.0f ? ceil(startPoint.y) : floor(startPoint.y) , startPoint.z < 0.0f ? ceil(startPoint.z) : floor(startPoint.z)); /*floor(startPoint)*/
    vec3 endCell ( (endPoint.x < 0.0f) ? ceil(endPoint.x) : floor(endPoint.x), (endPoint.y < 0.0f) ? ceil(endPoint.y) : floor(endPoint.y) , (endPoint.z < 0.0f) ? ceil(endPoint.z) : floor(endPoint.z)); /*floor(endPoint)*/

    vec3 direction = endPoint - startPoint;
    vec3 norm_direction = normalize(direction);

    // Determine the step direction for raycasting
    int stepX = sign(rayWOR.x);
    int stepY = sign(rayWOR.y);
    int stepZ = sign(rayWOR.z);

    // Determine the step direction for raycasting
    // Distance to nearest square side
    float near_x = (stepX >= 0) ? ((startCell.x + 1) * cellSize - startPoint.x) : (startPoint.x - (startCell.x * cellSize)) ;
    float near_y = (stepY >= 0) ? ((startCell.y + 1) * cellSize - startPoint.y) : (startPoint.y - (startCell.y * cellSize)) ;
    float near_z = (stepZ >= 0) ? ((startCell.z + 1) * cellSize - startPoint.z) : (startPoint.z - (startCell.z * cellSize)) ;

    const float EPSILON = 1e-6f;
    float maxX = (std::abs(norm_direction.x) > EPSILON) ? near_x / norm_direction.x : DBL_MAX;
    float maxY = (std::abs(norm_direction.y) > EPSILON) ? near_y / norm_direction.y : DBL_MAX;
    float masZ = (std::abs(norm_direction.z) > EPSILON) ? near_z / norm_direction.z : DBL_MAX;

    float dx = (std::abs(norm_direction.x) > EPSILON) ? cellSize / norm_direction.x : DBL_MAX;
    float dy = (std::abs(norm_direction.y) > EPSILON) ? cellSize / norm_direction.y : DBL_MAX;
    float dz = (std::abs(norm_direction.z) > EPSILON) ? cellSize / norm_direction.z : DBL_MAX;

    vec3 pos = startPoint; 

    int gridBoundX = std::abs(endCell.x - startCell.x);
    int gridBoundY = std::abs(endCell.y - startCell.y);
    int gridBoundZ = std::abs(endCell.z - startCell.z);

    int counter = 0;

    while(counter != (gridBoundX + gridBoundY + gridBoundZ)) {
      if(std::abs(maxX) < std::abs(maxY)) {
        if(std::abs(maxX) < std::abs(masZ)) {
          maxX += dx;
          pos.x += stepX;
        } else {
          masZ += dz;
          pos.z += stepZ;
        }
      } else {
        if(std::abs(maxY) < std::abs(masZ)) {
          maxY += dy;
          pos.y += stepY;
        } else {
          masZ += dz;
          pos.z += stepZ;
        }
      }
 
      if(pos.y < 0)
        break;

      if(blockExists(pos.x, pos.y, pos.z)) {
        return pos;
      }
      ++counter;
    }

    return vec3(0.0, -1.0, 0.0);
}