Finding quaternion representing the rotation from one vector to another

Quaternion q;
vector a = crossproduct(v1, v2);
q.xyz = a;
q.w = sqrt((v1.Length ^ 2) * (v2.Length ^ 2)) + dotproduct(v1, v2);

Don't forget to normalize q.

Richard is right about there not being a unique rotation, but the above should give the "shortest arc," which is probably what you need.

The problem as stated is not well-defined: there is not a unique rotation for a given pair of vectors. Consider the case, for example, where u = <1, 0, 0> and v = <0, 1, 0>. One rotation from u to v would be a pi / 2 rotation around the z-axis. Another rotation from u to v would be a pi rotation around the vector <1, 1, 0>.

Some of the answers don't seem to consider possibility that cross product could be 0. Below snippet uses angle-axis representation:

//v1, v2 are assumed to be normalized
Vector3 axis = v1.cross(v2);
if (axis == Vector3::Zero())
    axis = up();
else
    axis = axis.normalized();

return toQuaternion(axis, ang);

The toQuaternion can be implemented as follows:

static Quaternion toQuaternion(const Vector3& axis, float angle)
{
    auto s = std::sin(angle / 2);
    auto u = axis.normalized();
    return Quaternion(std::cos(angle / 2), u.x() * s, u.y() * s, u.z() * s);
}

If you are using Eigen library, you can also just do:

Quaternion::FromTwoVectors(from, to)

Working just with normalized quaternions, we can express Joseph Thompson's answer in the follwing terms.

Let q_v = (0, u_x, v_y, v_z) and q_w = (0, v_x, v_y, v_z) and consider

q = q_v * q_w = (-u dot v, u x v).

So representing q as q(q_0, q_1, q_2, q_3) we have

q_r = (1 - q_0, q_1, q_2, q_3).normalize()

According to the derivation of the quaternion rotation between two angles, one can rotate a vector u to vector v with

function fromVectors(u, v) {

  d = dot(u, v)
  w = cross(u, v)

  return Quaternion(d + sqrt(d * d + dot(w, w)), w).normalize()
}

If it is known that the vectors u to vector v are unit vectors, the function reduces to

function fromUnitVectors(u, v) {
  return Quaternion(1 + dot(u, v), cross(u, v)).normalize()
}

Depending on your use-case, handling the cases when the dot product is 1 (parallel vectors) and -1 (vectors pointing in opposite directions) may be needed.

The Generalized Solution

function align(Q, u, v)
    U = quat(0, ux, uy, uz)
    V = quat(0, vx, vy, vz)
    return normalize(length(U*V)*Q - V*Q*U)

To find the quaternion of smallest rotation which rotate u to v, use

align(quat(1, 0, 0, 0), u, v)

Why This Generalization?

R is the quaternion closest to Q which will rotate u to v. More importantly, R is the quaternion closest to Q whose local u direction points in same direction as v.

This can be used to give you all possible rotations which rotate from u to v, depending on the choice of Q. If you want the minimal rotation from u to v, as the other solutions give, use Q = quat(1, 0, 0, 0).

Most commonly, I find that the real operation you want to do is a general alignment of one axis with another.

// If you find yourself often doing something like
quatFromTo(toWorldSpace(Q, localFrom), worldTo)*Q
// you should instead consider doing 
align(Q, localFrom, worldTo)

Example

Say you want the quaternion Y which only represents Q's yaw, the pure rotation about the y axis. We can compute Y with the following.

Y = align(quat(Qw, Qx, Qy, Qz), vec(0, 1, 0), vec(0, 1, 0))

// simplifies to
Y = normalize(quat(Qw, 0, Qy, 0))

Alignment as a 4x4 Projection Matrix

If you want to perform the same alignment operation repeatedly, because this operation is the same as the projection of a quaternion onto a 2D plane embedded in 4D space, we can represent this operation as the multiplication with 4x4 projection matrix, A*Q.

I = mat4(
    1, 0, 0, 0,
    0, 1, 0, 0,
    0, 0, 1, 0,
    0, 0, 0, 1)

A = I - leftQ(V)*rightQ(U)/length(U*V)

// which expands to
A = mat4(
    1 + ux*vx + uy*vy + uz*vz, uy*vz - uz*vy, uz*vx - ux*vz, ux*vy - uy*vx,
    uy*vz - uz*vy, 1 + ux*vx - uy*vy - uz*vz, uy*vx + ux*vy, uz*vx + ux*vz,
    uz*vx - ux*vz, uy*vx + ux*vy, 1 - ux*vx + uy*vy - uz*vz, uz*vy + uy*vz,
    ux*vy - uy*vx, uz*vx + ux*vz, uz*vy + uy*vz, 1 - ux*vx - uy*vy + uz*vz)

// A can be applied to Q with the usual matrix-vector multiplication
R = normalize(A*Q)

//LeftQ is a 4x4 matrix which represents the multiplication on the left
//RightQ is a 4x4 matrix which represents the multiplication on the Right
LeftQ(w, x, y, z) = mat4(
    w, -x, -y, -z,
    x,  w, -z,  y,
    y,  z,  w, -x,
    z, -y,  x,  w)

RightQ(w, x, y, z) = mat4(
    w, -x, -y, -z,
    x,  w,  z, -y,
    y, -z,  w,  x,
    z,  y, -x,  w)