Getting ray using gluUnProject or inverted MVP matrix

Question

I've read a lot of topics here, on SO, opengl.org etc. Example how gluUnProject should work (from NeHe tutorial):

winX = (float)x;
winY = (float)viewport[3] - (float)y;
glReadPixels( x, int(winY), 1, 1, GL_DEPTH_COMPONENT, GL_FLOAT, &winZ );
gluUnProject( winX, winY, winZ, modelview, projection, viewport, &posX, &posY, &posZ);

The problem is that gluUnProject returns mostly zeros, if not zeros then it returns Camera position in World Coords. Value winZ I'm getting from glReadPixels is always correct, if I scale it like below:

float z_distance = ProjectionMatrix[14] / (winZ * -2.0 + 1.0 - ProjectionMatrix[10]);

I get distance between near plane and parallel plane that contains mouse click point. Look at the sample below:

gluUnProject( winX, winY, -1, modelview, projection, viewport, &posX1, &posY1, &posZ1);
gluUnProject( winX, winY, 0, modelview, projection, viewport, &posX2, &posY2, &posZ2);
gluUnProject( winX, winY, 1, modelview, projection, viewport, &posX3, &posY3, &posZ3);
gluUnProject( winX, winY, winZ, modelview, projection, viewport, &posX4, &posY4, &posZ4);

The result is always the same (posX1,posY1,posZ1)=(posX2,posY2,posZ2)=(posX3,posY3,posZ3)=(posX4,posY4,posZ4) and is either (0,0,0) or (Camera-Position.x,Camera-Position.y,Camera-Position.z)

After gluUnProject I tried to do it manually:

float NDC_x = 2.0f * winX / Width - 1.0f;
float NDC_y = 1.0f - 2.0f * winY / Height;
float winZ = 0;
glReadPixels(winX, (Height - winY), 1, 1, GL_DEPTH_COMPONENT, GL_FLOAT, &winZ);
vec4 point1 = inverse(ViewProjectionMatrix) * vec4(NDC_x, NDC_y, winZ, 1.0f);
vec4 point2 = inverse(ViewProjectionMatrix) * vec4(NDC_x, NDC_y, -1.0f, 1.0f);
vec4 point3 = inverse(ViewProjectionMatrix) * vec4(NDC_x, NDC_y, 0.0f, 1.0f);
vec4 point4 = inverse(ViewProjectionMatrix) * vec4(NDC_x, NDC_y, 1.0f, 1.0f);

Now all points coordinates are different and never zeros but some values, seems like real ones but I don't know how to handle it or I missed something.

BTW I used different approach to get ray (using angles) and it works perfectly. Check here.

Q1: Any suggestions what's wrong with gluUnProject ?

Q2: What is my mistake in manual implementation?

I thought that I messed with row\column-major format, but transposition of matrices had no effect. I read\write matrices as column-major format so last column is (ViewMatrix[12], ViewMatrix[13], ViewMatrix[14], ViewMatrix[15])

Answer 1

Firstly, we should convert winX and winY values in pixels to Normalized Device Coordinates.

To get world coordinates from pixel position on screen we should multiply inverse ViewProjection matrix by clip space vec3(NDC_x, NDC_y, some_z, 1.0f) Actually every pixel on screen corresponds to line in world space. So the question is what some_z value to use. But no matter what values we will use - final points will always be laying on our ray.

Tricky part is to correctly use received values. Let's look at four possible situations I've mentioned in opening post:

float NDC_x = 2.0f * winX / Width - 1.0f;
float NDC_y = 1.0f - 2.0f * winY / Height;
float winZ = 0;
glReadPixels(winX, (Height - winY), 1, 1, GL_DEPTH_COMPONENT, GL_FLOAT, &winZ);
vec4 point1 = inverse(ViewProjectionMatrix) * vec4(NDC_x, NDC_y, winZ, 1.0f);
vec4 point2 = inverse(ViewProjectionMatrix) * vec4(NDC_x, NDC_y, -1.0f, 1.0f);
vec4 point3 = inverse(ViewProjectionMatrix) * vec4(NDC_x, NDC_y, 0.0f, 1.0f);
vec4 point4 = inverse(ViewProjectionMatrix) * vec4(NDC_x, NDC_y, 1.0f, 1.0f);

Now we need to do perspective division (this part that is missing very often).

Note: code below uses 'point' with no numeration because its the same for any case.

point /= point.w; //perspective division for chosen point
vec3 RayDirection = vec3(point.x, point.y, point.z) - Camera->Position; 
 //subtracting  camera position from 'point' coordinates gives as ray direction, length 
 //depends on `z-value` used: '1.0f' gives us ray that stops on far plane, '0.0f' - ray 
 //will stop at near plane, all negative values return point between camera and near plane,
 //using value from glReadPixels(from 0 to 1) returns point between near and far plane
ScaledRay = RayDirection * scale; //scale original ray in 'scale' times
NormalizedRay = normalize(ScaledRay); //normalizing ray length to 1, we can normalize first
                //and then scale to whatever length we want(for example range of shooting)
RayEndPoint = NormalizedRay + Camera->Position; //adding Camera position gives endpoint

Now when we know how to handle values - same goes for gluUnProject:

GLint vp[4];
glGetIntegerv(GL_VIEWPORT, vp);
GLdouble wx = 0, wy = 0, wz = 0;
gluUnProject(winX, (vp[3] - winY), 0, mvmatrix, projmatrix, vp, &wx, &wy, &wz);
vec3 RayEndpoint = vec3((float)wx, (float)wy, (float)wz);
vec3 RayDirection = RayEndpoint - Camera->Position; 
vec3 ScaledRay = RayDirection * Scale;
vec3 NormalizedRay = normalize(ScaledRay);
RayEndPoint = NormalizedRay + Camera->Position;

Note: gluUnProject is unstable for me - mostly it returns zeros or #INF, but once per 10 launches of exe it works flawlessly (I meant I never change exe but just starting it again and again).