How to calculate look at point to move the camera

This will be confusing for me to explain so please bear with me.

I've already implemented most type of movements and rotations in my camera class, everything is working with the keyboard, now I want to implement the mouse. I capture the mouse movement like this:

#define SENSITIVITY 25.0f

void main(void) {
    (...)
    glutPassiveMotionFunc(processPassiveMotion);    
    glutWarpPointer(WINDOW_WIDTH / 2, WINDOW_HEIGHT / 2);
    glutSetCursor(GLUT_CURSOR_NONE);
    (...)
}

void processPassiveMotion(int x, int y) {
    int centerX = WINDOW_WIDTH / 2;
    int centerY = WINDOW_HEIGHT / 2;

    int deltaX = -1 * (x - centerX);
    int deltaY = -1 * (y - centerY);

    if(deltaX != 0 || deltaY != 0) {
        mainCamera.Rotate(deltaX / SENSITIVITY, deltaY / SENSITIVITY);

        glutWarpPointer(centerX, centerY);
    }
}

After everything I've read, I believe this is enough in my situation. However I must state that first I tried to call the Pitch() and Yaw() camera functions but it was a no go, I had to create an extra function to rotate both axis "at the same time".

That rotate function goes something like this:

#define DEG2RAD(a) (a * (M_PI / 180.0f))
#define SINDEG(a)  sin(DEG2RAD(a))
#define COSDEG(a)  cos(DEG2RAD(a))

void Camera::Rotate(GLfloat angleX, GLfloat angleY) {
    Reference = NormalizeVector(
        Reference * COSDEG(angleY) + UpVector * SINDEG(angleY)
    );

    Reference = NormalizeVector(
        Reference * COSDEG(angleX) - RightVector * SINDEG(angleX)
    );

    UpVector = CrossProduct(&Reference, &RightVector) * (-1);
    RightVector = CrossProduct(&Reference, &UpVector);
}

The Reference is the viewing direction, the point the camera is looking at. And since it's a normalized vector, it goes from -1.0 to 1.0. This vector, or point, is later used together with another vector (Position, which is the camera location) to calculate the real look at point to use in gluLookAt, like this:

void Camera::LookAt(void) {
    Vector3D viewPoint = Position + Reference;

    gluLookAt(
        Position.x, Position.y, Position.z,
        viewPoint.x, viewPoint.y, viewPoint.z,
        UpVector.x, UpVector.y, UpVector.z
    );
}

All vector operations above like +, - and * are overloaded of course.

Now I'm going to try to describe my problem...

The rotate function above works just fine in the sense that it correctly performs a pitch and yaw by using the mouse. However, those rotations don't look like the ones in First Person Shooter games. In those games, when one looks at sky and then looks left/right, one expects to keep looking at the sky. Imagining we are inside a sphere, a movement like that should "draw" a circle in the top part of the sphere.

But that's not what happens because that's not what a yaw does. A yaw movement will rotate around an arbitrary axis, which I think is the up vector in this situation. So, the problem is in the yaw movement because the pitch seems to work fine.

In other words, my code above can't keep the horizon leveled and that's what must happen cause that's happens in games when one looks at the sky and then look left/right, the horizon is always leveled. The same will not happen with my code, I look up and then left/right, and the horizon will be all twisted.

Did I make myself clear enough? I'm not sure how can I explain this any better. :( Hopefully it's enough for anyone to understand.

I'm not sure how can I fix this problem... How can I look left/right correctly after looking up/down, keeping the horizon leveled?

EDIT:

My rotate function code is taken from both the Yaw and Pitch functions which also exist so I can call those rotations independently. For reference purposes I'll add them below along with the Roll function too (which I'll probably never use, but in case I need it, it's there):

void Camera::Pitch(GLfloat angle) {
    Reference = NormalizeVector(
        Reference * COSDEG(angle) + UpVector * SINDEG(angle)
    );

    UpVector = CrossProduct(&Reference, &RightVector) * (-1);
}

void Camera::Yaw(GLfloat angle) {
    Reference = NormalizeVector(
        Reference * COSDEG(angle) - RightVector * SINDEG(angle)
    );

    RightVector = CrossProduct(&Reference, &UpVector);
}

void Camera::Roll(GLfloat angle) {
    RightVector = NormalizeVector(
        RightVector * COSDEG(angle) - UpVector * SINDEG(angle)
    );

    UpVector = CrossProduct(&Reference, &RightVector) * (-1);
}

Your problem appears to be in the statement:

UpVector = CrossProduct(&Reference, &RightVector) * (-1);

Once you've rotated Reference toward RightVector in the previous statement, their cross product will no longer result in an UpVector that gives you a horizontal horizon. Try it with your arms. Furthermore, Reference and RightVector are not separated by 90 degrees, so UpVector won't even be a unit vector either. (Finally, you should really just switch the order of the cross product for clarity, rather than multiplying by (-1).)

Honestly, if I were doing it, I would take a different approach. I don't see any logical reason why the two rotations have to be in one function. I also avoid explicit sines and cosines at all costs when working with vectors. I think what you really need is a function to Rotate About an Arbitrary Axis. If nothing else, it's very useful. Fortunately all the details are taken care of by Mr. Murray! If you implement this function, then it becomes very simple. Define a constant SkyVector that always points upward. Then in pseudocode,

AxisRotation( Vector vec, Vector axis, float angle ) {
    Vector result;

    // The axis is assumed to be normalized:  
    //    (just make sure you're not modifying the original)
    axis = NormalizeVector( &axis );

    // expanded for clarity:
    float u = axis.x;
    float v = axis.y;
    float w = axis.z;
    float x = vec.x;
    float y = vec.y;
    float z = vec.z;
    float c = cos(angle);
    float s = sin(angle);

    // Apply the formula verbatim from the linked page:
    result.x = u*(u*x + v*y + w*z)*(1.-c) + x*c + (-w*y + v*z)*s;
    result.y = v*(u*x + v*y + w*z)*(1.-c) + y*c + ( w*x - u*z)*s;
    result.z = w*(u*x + v*y + w*z)*(1.-c) + z*c + (-v*x + u*y)*s;

    return result;
}

Yaw(angleX) {
    Reference = AxisRotation( &Reference, &SkyVector, angleX );
    RightVector = NormalizeVector( CrossProduct( &Reference, &SkyVector ) );
    UpVector = CrossProduct( &RightVector, &Reference );
}

Pitch(angleY) {
    Reference = AxisRotation( &Reference, &RightVector, angleY );
    //RightVector doesn't change!
    UpVector = CrossProduct( &RightVector, &Reference );
}

If you go through that operation by operation, it should hopefully make some sense. Finally, I'll add that quaternions are really the 'correct' way to do this stuff and avoid gimbal lock, but I usually do pretty much exactly what you've done. You might have to check every now and then to make sure your vectors stay nice and perpendicular. Quaternions are more stable.

Edit: If the axis rotation function is overkill, you can still implement this with simple vectors and rotation matrices. The only thing is you'll have to start projecting things into the horizontal plane so that you can do the two rotations independently And it'll still take some sines and cosines. Your time is probably better spent implementing the axis rotation function!