What causes this distortion in my perspective projection at steep view angles?

Question

For various reasons – partly educational, partly technical – I'm writing a 3D game engine from scratch, using OpenGL 3.0; aiming for Core Profile, eventually, so of course I'm using a custom matrix stack and shaders.  I've got nearly all of that working, but I'm getting strange results and not sure where I'm making a mistake.

The problem:  using a perspective projection matrix (based on the formulas from the OpenGL 2.1 Reference), everything looks great… until I aim the camera up or down at a steep angle, when it gets horribly distorted, as shown in the screenshots below.

I've checked out the answers here (though the video is unavailable, so I'm not sure if it's the same sort of distortion), and double-checked my code against my references.  So far as I can tell the math is correct and I don't see any glaring typos, but I'm hoping a fresh pair of eyes can help me find the cause – thanks in advance.

Code

If there's anything else relevant you need to look at, just ask.

Vector3f Util::normalize(const Vector3f &v)
{
    float d = std::sqrt(dot(v, v));
    return d ? (v / d) : v;
}

Vector4f Util::buildQuaternion(const Angle &angle, const Vector3f &axis)
{
    Angle halfAngle = 0.5f * angle;
    float c = halfAngle.cos(), s = halfAngle.sin();
    return Vector4f(s * axis.x, s * axis.y, s * axis.z, c);
}

Vector3f Util::rotate(const Vector3f &p, const Vector4f &q)
{
    Vector4f r = qmul(qmul(q, Vector4f(p.x, p.y, p.z, 0.0f)), ~q);
    return Vector3f(r.x, r.y, r.z);
}

void Util::orbit(Vector3f &planet, const Vector3f &sun, const Vector3f &up, const Angle &yaw, const Angle &pitch)
{
    Vector3f f = planet - sun;
    Vector3f s = cross(up, f);
    Vector4f qx = buildQuaternion(pitch, normalize(s));
    Vector4f qy = buildQuaternion(yaw, normalize(up));
    Vector3f r = rotate(f, qmul(qx, qy));
    planet = r + sun;
}

Transform &Transform::makeRotate(const Vector4f &q)
{
    float xx = q.x * q.x;
    float xy = q.x * q.y;
    float xz = q.x * q.z;
    float xw = q.x * q.w;
    float yy = q.y * q.y;
    float yz = q.y * q.z;
    float yw = q.y * q.w;
    float zz = q.z * q.z;
    float zw = q.z * q.w;
    float ww = q.w * q.w;

    float *m = m_matrix;
    m[ 0] = 1-2*(yy+zz); m[ 4] =   2*(xy-zw); m[ 8] =   2*(xz+yw); m[12] = 0.0;
    m[ 1] =   2*(xy+zw); m[ 5] = 1-2*(xx+zz); m[ 9] =   2*(yz-xw); m[13] = 0.0;
    m[ 2] =   2*(xz-yw); m[ 6] =   2*(yz+xw); m[10] = 1-2*(xx+yy); m[14] = 0.0;
    m[ 3] = 0.0;         m[ 7] = 0.0;         m[11] = 0.0;         m[15] = 1.0;
    return *this;
}

Transform &Transform::makeRotate(const Angle &angle, const Vector3f &axis)
{
    return makeRotate(Util::buildQuaternion(angle, axis));
}

Transform &Transform::makePerspective(const Angle &fov, float aspect, float zNear, float zFar)
{
    float f = (fov * 0.5f).cot();
    float e = f / aspect;
    float g = (zFar + zNear) / (zNear - zFar);
    float h = 2.0f * zFar * zNear / (zNear - zFar);

    float *m = m_matrix;
    m[ 0] =  e ; m[ 4] = 0.0; m[ 8] =  0.0; m[12] = 0.0;
    m[ 1] = 0.0; m[ 5] =  f ; m[ 9] =  0.0; m[13] = 0.0;
    m[ 2] = 0.0; m[ 6] = 0.0; m[10] =   g ; m[14] =  h ;
    m[ 3] = 0.0; m[ 7] = 0.0; m[11] = -1.0; m[15] = 0.0;
    return *this;
}

Transform &Transform::makeLookAt(const Vector3f &eye, const Vector3f &center, const Vector3f &up)
{
    Vector3f f = Util::normalize(center - eye);
    Vector3f u = Util::normalize(up);
    Vector3f s = cross(f, u);
    u = cross(Util::normalize(s), f);

    float *m = m_matrix;
    m[ 0] =  s.x; m[ 4] =  s.y; m[ 8] =  s.z; m[12] = 0.0;
    m[ 1] =  u.x; m[ 5] =  u.y; m[ 9] =  u.z; m[13] = 0.0;
    m[ 2] = -f.x; m[ 6] = -f.y; m[10] = -f.z; m[14] = 0.0;
    m[ 3] =  0.0; m[ 7] =  0.0; m[11] =  0.0; m[15] = 1.0;
    return translate(-eye);
}

void Camera::orbit3P(const Angle &yaw, const Angle &pitch)
{
    /// 3rd-person orbit - move camera around subject
    Util::orbit(m_eye, m_center, m_up, yaw, pitch);
    m_dirty = true;
}

void Camera::orbit1P(const Angle &yaw, const Angle &pitch)
{
    /// 1st-person orbit - change camera direction
    Util::orbit(m_center, m_eye, m_up, -yaw, pitch);
    m_dirty = true;
}

Screenshots

Horizontal - all is well

About 45° - distortion begins to appear

Almost 90° - what is this I don't even

(The ortho overlay was to help visualize the camera eye and center locations.)

Answer 1

Inside makeLookAt contains your error:

Vector3f f = Util::normalize(center - eye);
Vector3f u = Util::normalize(up);
Vector3f s = cross(f, u);
u = cross(Util::normalize(s), f);

s is not normalized.

It should be:

Vector3f f = Util::normalize(center - eye);
Vector3f u = Util::normalize(up);
Vector3f s = Util::normalize(cross(f, u));
u = cross(s, f); // no norm needed because s,f both unit length and orthogonal

If up and f (forward) are perpendicular, then s is already (approximately) unit length. That is why it works for looking horizontally but not vertically.