OpenGL: gluUnProject returns wrong coordinates

I set up my scene (Android app, OpenGL ES) like this:

GLU.gluPerspective(gl, 60, viewRatio, 0.1f, 1000.0f);
// ...
GLU.gluLookAt(gl, cameraX, cameraY, cameraZ, cameraX, cameraY, cameraZ - 1f, 0f, 1f, 0f);
// Scene scaled down because object positions are in [-10000..10000] for x & y
gl.glScalef(0.001f, 0.001f, 1.0f);

The scene is rendered well, it contains quads at z=-10 and one giant background quad at z=-30. I am now trying to implement ray picking, like this (taken from this thread:

public float[] unproject(float rx, float ry, float rz) {
    float[] xyzw = {0, 0, 0, 0};
    int[] viewport = {0, 0, mDisplayWidth, mDisplayHeight};
    android.opengl.GLU.gluUnProject(rx, ry, rz, mMatrixGrabber.mModelView, 0, mMatrixGrabber.mProjection, 0, viewport, 0, xyzw, 0);
    xyzw[0] /= xyzw[3];
    xyzw[1] /= xyzw[3];
    xyzw[2] /= xyzw[3];
    xyzw[3] = 1;
    return xyzw;
}

Taps onto the screen are then tested this way:

unproject(tapX, mDisplayHeight - tapY, BACKGROUND_Z);

Expected behaviour: Returned X,Y coordinates resemble the tapped point at depth = -30 (BACKGROUND_Z)

What it actually does: Returns correct X,Y coordinates when tapped at the screen center; But the distance between screen center and tapped-point is not interpreted correctly - it seems to be scaled down by some factor. E.g. tap on screen center gives [0,0], tap on left screen edge should give ~[-3600,0], but instead gives [-1.13,0]. Panning the view such that the previously touched point is in the screen center gives ~[-3600,0], as it should.

Can this be due to the scaling of the scene by 0.001? I have tested various configurations (moving scaling before gluLookAt, or entirely removing it), but the problem persists :(

Sample data output for touch at left screen edge:

Touch point passed to unproject():

x=3, y = 554 (0,0 at bottom left; w=720, h=1280)

Projection matrix:

[2.8482616, 0.0,       0.0,         0.0, 
 0.0,       1.7320509, 0.0,         0.0, 
 0.0,       0.0,      -1.0001999,  -1.0, 
 0.0,       0.0,      -0.20001999,  0.0]

Model matrix:

[0.0010, 0.0,    0.0, 0.0, 
 0.0,    0.0010, 0.0, 0.0, 
 0.0,    0.0,    1.0, 0.0, 
 0.0,    0.0,    1.0, 1.0]

unproject() output:

[-1.1232367, -0.11801138, -1.0032262, 1.0]

What comes to my mind here is that the returned z value seems to be wrong - the z value passed into unproject() is -30 (as mentioned above).

Solution

Let's see -- I wasn't able to reproduce the values you got, however I think this does not matter.

You must not call unproject with z-values not within [-1;+1]. The first thing unproject must do is reverse the windowing transform, and after that, the coordinates must be in NDC-space. In OpenGL, this is the cube [-1;+1]^3 (songho contains some nice graphics). Calling unproject with -30 means that the point must have been before the near plane. Calling it with z=-1.0f results in the near plane, and z=+1.0f in the far plane.

How to choose z s.t. after all matrix inversions, z=-30 holds -- I don't know.