Ray Cast / Ray Tracing

Hello, I would like to ask some theoretical Questions here , combined with some pseudo code to explain what i need or mean. So recently i implemented traditional ray casting from the camera’s position into the scene. I will devide my post in 2 parts:

PART ONE:
Ray Unprojecting from far to near plane


        float mouse_x = (float)InputState::getMOUSE_X();
	float mouse_y = WINDOW_HEIGHT - (float)InputState::getMOUSE_Y();
	glm::vec4 viewport = glm::vec4(0.0f, 0.0f, WINDOW_WIDTH, WINDOW_HEIGHT);
	this->ray_start = glm::unProject(glm::vec3(mouse_x, mouse_y, 0.0f), camera->getViewMatrix(), Projection, viewport);
	this->ray_end = glm::unProject(glm::vec3(mouse_x, mouse_y, 1.0f), camera->getViewMatrix(), Projection, viewport);


Then i implemented a simple function to detect if the ray has collided with a bounding sphere.


       glm::vec3 vSphereCenter = Sphere->getCenter();
	glm::vec3 vA = this->ray_start;
	glm::vec3 vB = this->ray_end;
	float fSphereRadius = Sphere->getRadius();

	glm::vec3 vDirToSphere = vSphereCenter - vA;
	glm::vec3 vLineDir = glm::normalize(vB - vA);
	float fLineLength = glm::distance(vA, vB);
	float t = glm::dot(vDirToSphere, vLineDir);
	glm::vec3 vClosestPoint;
	
	if (t <= 0.0f)
		vClosestPoint = vA;
	else if (t >= fLineLength)
		vClosestPoint = vB;
	else
		vClosestPoint = vA + vLineDir*t;

	return glm::distance(vSphereCenter, vClosestPoint) <= fSphereRadius;

So far so good. The code above works. Here starts my broken “theory” / “questions”. Let’s say that we have a scene that consists of models. All models have some auto-generated bounding spheres. So we want to check whitch object the player has “clicked” on , or the object the ray intersects with.First when i was thinking about this. I didnt think that several models can be in a line (in each axis), thats why you should loop over all the models that have been hit and find the closest to the camera. Look at the pseudo code mix below, added some comments as well.


//Here i present you the Pseudo C++ mix :D
std::vector<models> scene_models;   //All Models in the scene
std::vector<models> selected_models; //All models that have been hit by the ray

//Loop over every model find where the ray intersects
for(i,scene_models.size())
{
  if(RayIntersectsWith(model[i]->getBoundSphere())
      {
         selected_models.pushback(model[i])
      }
}
//Loop over all intersected model and find the one closest to the camera
for(i,selected_models.size())
{
  //get the distance(camera->pos,selected_models[i]->getBoundSphere()->getCenter());
  //Find the one closest to the camera = closest_model_to_the_camera;
}
//Do something with closest_model_to_the_camera

My Question is as follows: Is this the right way to go ? Would it work ? Is there a more efficient way to apply this ?

PART TWO
This part is not connected with the part one. Here i wanted to ask how would one implement terrain responsive ray. A ray that can return the world coords on a 3D Terrain. Currently i have something that works , but not that well. At the moment i am using the ray and performing a binary search to find the point on the terrain


//What i mean by binary search:
Take the starting point of the ray. Set Ray Lenght. 
1. Take the middle point of the ray.
2. Determine if that point is above of below the terrain
3. if the point is below the terrain take the upper half / else take the lower half of the ray
4. Repeat N times steps 1-4. (Recursivly)
//This is not very accurate nor efficient way but does find the exact world position on the terrain, the more times(the bigger the N) you repeat the process the more accurate the result is.

. And it does work , the only problem is that i have to use big camera angle (~near 90 degrees). The terrain i am using is generated from a height map.
My question is: Is there some stable, reliable way to return the world coordinates on terrain from the ray.

[QUOTE=Asmodeus;1271984]


//Here i present you the Pseudo C++ mix :D
std::vector<models> scene_models;   //All Models in the scene
std::vector<models> selected_models; //All models that have been hit by the ray

//Loop over every model find where the ray intersects
for(i,scene_models.size())
{
  if(RayIntersectsWith(model[i]->getBoundSphere())
      {
         selected_models.pushback(model[i])
      }
}
//Loop over all intersected model and find the one closest to the camera
for(i,selected_models.size())
{
  //get the distance(camera->pos,selected_models[i]->getBoundSphere()->getCenter());
  //Find the one closest to the camera = closest_model_to_the_camera;
}
//Do something with closest_model_to_the_camera

My Question is as follows: Is this the right way to go ? Would it work ? Is there a more efficient way to apply this ?[/QUOTE]

I’m not entirely sure if it would be more efficient in your case, but couldn’t you create a 2D array (full of 0’s) representing what the camera is looking at and put a unique index at the locations where there’s supposed to be a model? The index could then be used in a hash/lookup table to find a pointer directing to the model(s). That way, all you’d have to do is to find the 2D position of the mouse, transform that position into indices, look up if there’s something in the array at those indices and run your ‘closest to camera’ function only on the models represented by the index. In cases where you have loads of model, you’d at least save the time you’d waste iterating on non-useful models.

Yea its possible. Also keep in mind that the vector<ofmodels> wont represent all models in the world , only those closest to the camer / or those that will be visible by it. So they wont be that much. Also alot of the scene is usually static. I was just asking if my approach is correct ? Am i on the right path ?

If you only need to know if it works, then yes, I do believe that your method would work.

One question tho’, i am having not very good time calculating the bounding sphere. Here is the sample snippet, the sphere does seem correct but i dont think that the center is calculated correclty. For example with a simple tree the sphere tends to be around the stem of the tree, same goes for a human-type mesh the sphere’s center is around the feet/knees i would expect it to be above the waist


void BoundingSphere::calculateBoundingSphere(std::vector<glm::vec3> vertices)
{
	vec3 _center = vec3(0.0f, 0.0f, 0.0f);
	for (int i = 0; i<vertices.size(); i++) {
		_center += vertices[i];
	}
	_center /= (float)vertices.size();

	float _radius = 0.0f;
	for (int i = 0; i<vertices.size(); i++) {
		vec3 v = vertices[i] - _center;
		float distSq = glm::length(v);
		if (distSq > _radius)
			_radius = distSq;
	}
	_radius = sqrtf(_radius);
	this->center = _center;
	this->radius = _radius;
}

EDIT: Also having mostly uneven meshes i should probably move to OBB or AABB

Unless I’m misunderstanding your approach, it’s incorrrect. If the ray intersects the terrain in multiple locations, it won’t necessarily find the intersection which is closest to the viewpoint.

You can’t do this by divide-and-conquer unless you have quadtrees (i.e. mipmaps) containing the minimum and maximum heights for each node. Even then, you may need to recurse into the nearest node first before determining that there’s no collision there then recursing into the farther nodes. IOW, in the worst case (for a ray tangential to the surface) it’s O(n) not O(log(n)).

You know (or presumably can easily find out) the start and end points of the ray in world coordinates, which can be interpolated to find any point on the ray.

Your calculation will find a point which is biased toward whichever side of the model has more vertices. There are various algorithms for finding the minimum bounding sphere or a close approximation to it (Link), but they’re significantly more complex. A simple alternative is to use the midpoint of the bounding box.

From your descriptions of the behaviour with a human and a tree, it sounds as if the vertical component is inverted. Are you neglecting to apply any relevant transformations to the centroid?

Well i apply scaling to the sphere and transform its position if thats what you mean. While i was waiting for new responses here i managed to create simple AABB generation for meshes and then found a very useful (custom) algorithm that finds intersection between ray and obb. The cool thing is that you can pass a AABB + ModelMatrix and the function that checks the intersection will convert it into an OBB , then check for intersection and return the value. Its handy and pretty accurate. One thing i am missing is that i found out that if i pass only TranfRotation (its all working smooth and fine, pixel perfect detection) but if i pass as model matrix TransfRotation*Scaling , the intersection is detected as if the obb is downscaled instead of the opposite.
Also forgot to mention that i scale the aabb max and min with the scaling factor of the mesh. (i assume this is correct since if i dont do this the box’s size is with the mesh’s base one)

If you’re transforming the ray into the coordinate system of the AABB, you need to use the inverse matrix.

If you’re composing the matrix from primitive transformations, you can avoid a generalised matrix inverse by using (AB)-1=B-1A-1 and inverting the individual transformations. Rotation can be inverted by negating the angle, translation by negating the vector, scaling by using the reciprocal of the scale factors.

As far as i can see the function does not take in account the scaling. Just a sample snippet from the code. As far as i can see the snippet below shows extrapolation of position and rotation from the matrix.


        glm::vec3 OBBposition_worldspace(ModelMatrix[3].x, ModelMatrix[3].y, ModelMatrix[3].z);
	glm::vec3 delta = OBBposition_worldspace - ray_origin;
	{
              
		glm::vec3 xaxis(ModelMatrix[0].x, ModelMatrix[0].y, ModelMatrix[0].z); 
		float e = glm::dot(xaxis, delta);
		float f = glm::dot(ray_direction, xaxis);

                //Where below calculations are processed for axis y,z

About my Part 2 question for the terrain , could you elaborate ?

“You know (or presumably can easily find out) the start and end points of the ray in world coordinates, which can be interpolated to find any point on the ray.”
I have the start and end points of the ray. But what do you mean by interpolating ?

Any scale factor will be reflected in the magnitudes of e and f (the dot product of two vectors is equal to the product of their magnitudes multiplied by the cosine of the angle between them).

If the code is written on the assumption that the axes are unit-length vectors, scaling will interfere with that. Or the scale factors may end up cancelling out.

[QUOTE=Asmodeus;1272004]
About my Part 2 question for the terrain , could you elaborate ?

“You know (or presumably can easily find out) the start and end points of the ray in world coordinates, which can be interpolated to find any point on the ray.”
I have the start and end points of the ray. But what do you mean by interpolating ?[/QUOTE]

For instance, when you find the midpoint, if you have the start and end points in world coordinates, then the average of those two values will be the midpoint in world coordinates.

Well is there any nice terrain-ray intersection algorithm available out there ? Since mine is just for testing and not accurate enought.
I will probably use both sphere and box intersections. So i will need a bit “loose” sphere around the mesh and if the ray intersects i will pass to the more accurate box intersection. That should save some significant time when testing more objects.

EDIT: I was wondering since i have generated a AABB , how could i convert that to OBB (this is done in the intersection check function , as i pass the matrices , but i was wondering how would it look like outside of that)

Project the ray onto the ground plane and trace the resulting line through the height map (as if you were drawing a line on a bitmap surface, e.g. Bresenham’s algorithm or DDA) starting at the edge nearest the viewpoint.

If you know the minimum and maximum height values, find the intersection between the ray and the corresponding planes. Then you only need to trace a line between those two points rather than across the entire height map.

If you want better average-case performance at the expense of complexity, generate two mipmaps for the minimum and maximum heights. You effectively then have a series of successively finer grids of AABBs bounding the terrain. You only need to examine the higher-resolution data when the lower-resolution data indicates a potential intersection.

I will leave out the terrain, for now cause i have something thats working to some extend. But i was looking over and coudn’t find any meaningful tutorial on how to generate OBBs

As for my last post i am still searching a nice OBB explanation. How its done more in-depth information. Can it be generated from a AABB. If possible some samples . Thanks

Well, an OBB (oriented bounding box) is just a bounding box which isn’t necessarily axis-aligned.

An OBB is typically an AABB plus a transformation. Testing for intersection with a ray can be performed by transforming the ray by the inverse transformation then using an AABB test.

Both of these are special cases of a convex hull. The specialisation just allows some of the calculations to be simplified. E.g. for an AABB, the calculation of the dot product with a face normal involves extracting one of the three coordinates (in essence, you’re multiplying by values which are zero or one, then adding three values of which two are known to be zero).

The convex hull of a set of points can be found deterministically, using e.g. Quickhull or other methods. This can be used as the basis for finding an OBB or bounding sphere (you only need to consider the vertices of the convex hull rather than all vertices).

Aham , but wait. So this is what i’ve got so far (from reading around). AABB is axil aligned , you can apply transformations or scaling, same applies for OBB but you can add Rotations ? My question was more of how to actually calculate one. I found various examples . Simple(on first glance) Code snippet i found. Once i started tracing the methods included in this constructor sh*t got real. The Covariance calc here invokes server other complex methods


OBB::Set( const Vector3* points, unsigned int nPoints )
{
    // //ASSERT( points );

    Vector3 centroid;

    // compute covariance matrix
    Matrix33 C;
    ComputeCovarianceMatrix( C, centroid, points, nPoints );

    // get basis vectors
    Vector3 basis[3];
    GetRealSymmetricEigenvectors( basis[0], basis[1], basis[2], C );
    mRotation.SetColumns( basis[0], basis[1], basis[2] );

    Vector3 min(FLT_MAX, FLT_MAX, FLT_MAX);
    Vector3 max(FLT_MIN, FLT_MIN, FLT_MIN);

    // compute min, max projections of box on axes
    // for each point do 
    unsigned int i;
    for ( i = 0; i < nPoints; ++i )
    {
        Vector3 diff = points[i]-centroid;
        for (int j = 0; j < 3; ++j)
        {
            float length = diff.Dot(basis[j]);
            if (length > max[j])
            {
                max[j] = length;
            }
            else if (length < min[j])
            {
                min[j] = length;
            }
        }
    }

    // compute center, extents
    mCenter = centroid;
    for ( i = 0; i < 3; ++i )
    {
        mCenter += 0.5f*(min[i]+max[i])*basis[i];
        mExtents[i] = 0.5f*(max[i]-min[i]);
    }

}   // End of OBB::Set()


Then i find another example on stackoverflow : c++ - obb rotation calculation - Stack Overflow. Here the OP uses completely different approach to construct the box.
Just the problem is i cant find ONE definitive way of doing this.

The reason for that is that computing an optimal OBB (for some definition of “optimal”; there’s more than one) would be insanely slow for any non-trivial mesh, so you’re left with a choice between various algorithms, some of which are closer to optimal, some are faster, some are simpler.

OTOH, finding the minimum AABB is trivial. While finding the convex hull is more complex, you at least have the advantage that there’s only one correct answer; the different algorithms just strike a different balance between performance and simplicity.

BTW, the algorithm you quote appears to be principal component analysis (PCA), which could equally be applied to finding an oriented bounding ellipsoid. Roughly speaking, It’s a multi-dimensional equivalent of least-squares line fitting.

Which method would you recommend. And if possible to apply some examples. I myself am new to this, not entirely sure if i understand the complete math involved but i started to clear some things up. Maybe i will better off with a simpler method, just for starters. Thanks :slight_smile:

Whichever one works best for your data.

E.g. PCA is normally used with data which follows something approximating a Gaussian distribution. The vertex coordinates of a typical mesh are probably a long way from that, so I’m not sure how well it will work in practice.

A fairly simple option for OBBs is to just test multiple orientations and use whichever one produces the smallest box. You don’t need to try very many to find an orientation that’s not too far off (or at least avoid the worst case, where all of your vertices are in a line which happens to be the box’s diagonal).