Custom lookAt function goes wild

Question

I have written a custom lookAt function based on a lot of posts from all over the net, and it works very nice... except when a rotation (which is stored in a quaternion) crosses some 'threshold'. However, it works flawless when the camera's position is on the same horizontal plane as the model I'm lookingAt

On the below GIF, all I do is hold a key to lookAt the model together with a left arrow key to move the camera along the right vector:

My following understanding about how this function should look like (in a brief):

1 .Calculate a quaternion rotation needed to rotate an object from it's facing vector and a target vector.

1.a. Calculate an angle between those vectors by taking their dot product and a perpendicular axis by taking their cross product.

1.b Convert angle-axis info to a quaternion using this formula.

2 .Multiply the main quaternion rotation by the quaternion received above.

If I'm right, I believe that I may have a bug in my code, but I'm not sure. My code below:

lookAt: function(vec) {
    if (this.position.isEqual(vec)) {
        return;
    }
    this.rotation.multiply(new Quaternion().twoVecToQuat(this.rotation.getForwardVector(), this.position.sub(vec)));
    MVMatrix.setRotation(this.rotation.quaternionToMatrix());
},

// And the twoVecToQuat function based on Ogre3D quaternion library:

twoVecToQuat: function(v0, v1) {
    v0.normalize();
    v1.normalize();

    var dot = v0.dot(v1);
    // If dot === 1, vectors are the same
    if (dot >= 1) {
        return this.makeIdentity();
    }
    // If vectors are opposite
    if (dot < (1.0e-06 - 1)) {
        var fallbackAxis = new vec3().make(1, 0, 0);
        fallbackAxis.crossSelf(this);
        // If collinear, pick another vector
        if (Math.areScalarsSimilar(fallbackAxis.length(), 0)) {
            fallbackAxis.make(0, 1, 0);
            fallbackAxis.crossSelf(this);
        }
        fallbackAxis.normalize();
        this.axisToQuaternion(Math.PI, fallbackAxis);
    }
    else {
        var sDot = Math.sqrt((1 + dot) * 2);
        var invSDot = 1 / sDot;
        var tempVec = v0.cross(v1);

        this.x = tempVec.x * invSDot;
        this.y = tempVec.y * invSDot;
        this.z = tempVec.z * invSDot;
        this.w = sDot * 0.5;

        this.normalize();
    }

    return this;
},

UPDATE: Slin's solutions proved to be working (that one about calculating lookAt matrix), but I still don't know where is the bug in the quaternions library code. Working lookAt function below:

lookAt: function(at) {
    if (this.position.isEqual(at)) {
        return;
    }
    var zAxis = this.position.sub(at).normalize();
    var up = upVec(); // World's up vector instead of camera's up vector will unroll the camera
    var xAxis = up.cross(zAxis).normalize();
    var yAxis = zAxis.cross(xAxis);
    MVMatrix.make(
        xAxis.x,            xAxis.y,            xAxis.z,
        yAxis.x,            yAxis.y,            yAxis.z,
        zAxis.x,            zAxis.y,            zAxis.z,
        this.position.x,    this.position.y,    this.position.z);

    this.rotation.matrixToQuaternion(MVMatrix);
    this.copyEulerAngles();
},

Answer 1

I am not sure but this does not look like the formula you posted:

 var sDot = Math.sqrt((1 + dot) * 2);
 var invSDot = 1 / sDot;
 var tempVec = v0.cross(v1);

 this.x = tempVec.x * invSDot;
 this.y = tempVec.y * invSDot;
 this.z = tempVec.z * invSDot;
 this.w = sDot * 0.5;

 this.normalize();

Unless I do not understand some of the sine/cosine relationships or approximations, it should instead be something like this:

 var halfCosine = Math.sqrt((1 + dot) * 0.5);
 var halfSine = Math.sqrt(1 - halfCosine*halfCosine);
 var tempVec = v0.cross(v1);

 this.x = tempVec.x * halfSine;
 this.y = tempVec.y * halfSine;
 this.z = tempVec.z * halfSine;
 this.w = halfCosine;

 this.normalize();

But you could also just set your quaternion or rotation matrix directly instead of multiplying to it to get the same result. Calculating a look at rotation matrix is actually really simple as all you need is the direction to the target and two orthogonal vectors representing the up and right direction. Check this for more information: https://stackoverflow.com/questions/349050/calculating-a-lookat-matrix