transformationmatrix.cpp

linux下开源浏览器WebKit的源码,市面上的很多商用浏览器都是移植自WebKit

    v3Scale(row[0], 1.0);

    // Compute XY shear factor and make 2nd row orthogonal to 1st.
    result.skewXY = v3Dot(row[0], row[1]);
    v3Combine(row[1], row[0], row[1], 1.0, -result.skewXY);

    // Now, compute Y scale and normalize 2nd row.
    result.scaleY = v3Length(row[1]);
    v3Scale(row[1], 1.0);
    result.skewXY /= result.scaleY;

    // Compute XZ and YZ shears, orthogonalize 3rd row.
    result.skewXZ = v3Dot(row[0], row[2]);
    v3Combine(row[2], row[0], row[2], 1.0, -result.skewXZ);
    result.skewYZ = v3Dot(row[1], row[2]);
    v3Combine(row[2], row[1], row[2], 1.0, -result.skewYZ);

    // Next, get Z scale and normalize 3rd row.
    result.scaleZ = v3Length(row[2]);
    v3Scale(row[2], 1.0);
    result.skewXZ /= result.scaleZ;
    result.skewYZ /= result.scaleZ;
 
    // At this point, the matrix (in rows[]) is orthonormal.
    // Check for a coordinate system flip.  If the determinant
    // is -1, then negate the matrix and the scaling factors.
    v3Cross(row[1], row[2], pdum3);
    if (v3Dot(row[0], pdum3) < 0) {
        for (i = 0; i < 3; i++) {
            result.scaleX *= -1;
            row[i][0] *= -1;
            row[i][1] *= -1;
            row[i][2] *= -1;
        }
    }
 
    // Now, get the rotations out, as described in the gem.
    
    // FIXME - Add the ability to return either quaternions (which are
    // easier to recompose with) or Euler angles (rx, ry, rz), which
    // are easier for authors to deal with. The latter will only be useful
    // when we fix https://bugs.webkit.org/show_bug.cgi?id=23799, so I
    // will leave the Euler angle code here for now.

    // ret.rotateY = asin(-row[0][2]);
    // if (cos(ret.rotateY) != 0) {
    //     ret.rotateX = atan2(row[1][2], row[2][2]);
    //     ret.rotateZ = atan2(row[0][1], row[0][0]);
    // } else {
    //     ret.rotateX = atan2(-row[2][0], row[1][1]);
    //     ret.rotateZ = 0;
    // }
    
    double s, t, x, y, z, w;

    t = row[0][0] + row[1][1] + row[2][2] + 1.0;

    if (t > 1e-4) {
        s = 0.5 / sqrt(t);
        w = 0.25 / s;
        x = (row[2][1] - row[1][2]) * s;
        y = (row[0][2] - row[2][0]) * s;
        z = (row[1][0] - row[0][1]) * s;
    } else if (row[0][0] > row[1][1] && row[0][0] > row[2][2]) { 
        s = sqrt (1.0 + row[0][0] - row[1][1] - row[2][2]) * 2.0; // S=4*qx 
        x = 0.25 * s;
        y = (row[0][1] + row[1][0]) / s; 
        z = (row[0][2] + row[2][0]) / s; 
        w = (row[2][1] - row[1][2]) / s;
    } else if (row[1][1] > row[2][2]) { 
        s = sqrt (1.0 + row[1][1] - row[0][0] - row[2][2]) * 2.0; // S=4*qy
        x = (row[0][1] + row[1][0]) / s; 
        y = 0.25 * s;
        z = (row[1][2] + row[2][1]) / s; 
        w = (row[0][2] - row[2][0]) / s;
    } else { 
        s = sqrt(1.0 + row[2][2] - row[0][0] - row[1][1]) * 2.0; // S=4*qz
        x = (row[0][2] + row[2][0]) / s;
        y = (row[1][2] + row[2][1]) / s; 
        z = 0.25 * s;
        w = (row[1][0] - row[0][1]) / s;
    }

    result.quaternionX = x;
    result.quaternionY = y;
    result.quaternionZ = z;
    result.quaternionW = w;
    
    return true;
}

// Perform a spherical linear interpolation between the two
// passed quaternions with 0 <= t <= 1
static void slerp(double qa[4], const double qb[4], double t)
{
    double ax, ay, az, aw;
    double bx, by, bz, bw;
    double cx, cy, cz, cw;
    double angle;
    double th, invth, scale, invscale;

    ax = qa[0]; ay = qa[1]; az = qa[2]; aw = qa[3];
    bx = qb[0]; by = qb[1]; bz = qb[2]; bw = qb[3];

    angle = ax * bx + ay * by + az * bz + aw * bw;

    if (angle < 0.0) {
        ax = -ax; ay = -ay;
        az = -az; aw = -aw;
        angle = -angle;
    }

    if (angle + 1.0 > .05) {
        if (1.0 - angle >= .05) {
            th = acos (angle);
            invth = 1.0 / sin (th);
            scale = sin (th * (1.0 - t)) * invth;
            invscale = sin (th * t) * invth;
        } else {
            scale = 1.0 - t;
            invscale = t;
        }
    } else {
        bx = -ay;
        by = ax;
        bz = -aw;
        bw = az;
        scale = sin(piDouble * (.5 - t));
        invscale = sin (piDouble * t);
    }

    cx = ax * scale + bx * invscale;
    cy = ay * scale + by * invscale;
    cz = az * scale + bz * invscale;
    cw = aw * scale + bw * invscale;

    qa[0] = cx; qa[1] = cy; qa[2] = cz; qa[3] = cw;
}

// End of Supporting Math Functions

TransformationMatrix& TransformationMatrix::scale(double s)
{
    return scaleNonUniform(s, s);
}

TransformationMatrix& TransformationMatrix::rotateFromVector(double x, double y)
{
    return rotate(rad2deg(atan2(y, x)));
}

TransformationMatrix& TransformationMatrix::flipX()
{
    return scaleNonUniform(-1.0f, 1.0f);
}

TransformationMatrix& TransformationMatrix::flipY()
{
    return scaleNonUniform(1.0f, -1.0f);
}

FloatPoint TransformationMatrix::projectPoint(const FloatPoint& p) const
{
    // This is basically raytracing. We have a point in the destination
    // plane with z=0, and we cast a ray parallel to the z-axis from that
    // point to find the z-position at which it intersects the z=0 plane
    // with the transform applied. Once we have that point we apply the
    // inverse transform to find the corresponding point in the source
    // space.
    // 
    // Given a plane with normal Pn, and a ray starting at point R0 and
    // with direction defined by the vector Rd, we can find the
    // intersection point as a distance d from R0 in units of Rd by:
    // 
    // d = -dot (Pn', R0) / dot (Pn', Rd)
    
    double x = p.x();
    double y = p.y();
    double z = -(m13() * x + m23() * y + m43()) / m33();

    double outX = x * m11() + y * m21() + z * m31() + m41();
    double outY = x * m12() + y * m22() + z * m32() + m42();

    double w = x * m14() + y * m24() + z * m34() + m44();
    if (w != 1 && w != 0) {
        outX /= w;
        outY /= w;
    }

    return FloatPoint(static_cast<float>(outX), static_cast<float>(outY));
}

FloatQuad TransformationMatrix::projectQuad(const FloatQuad& q) const
{
    FloatQuad projectedQuad;
    projectedQuad.setP1(projectPoint(q.p1()));
    projectedQuad.setP2(projectPoint(q.p2()));
    projectedQuad.setP3(projectPoint(q.p3()));
    projectedQuad.setP4(projectPoint(q.p4()));
    return projectedQuad;
}

FloatPoint TransformationMatrix::mapPoint(const FloatPoint& p) const
{
    double x, y;
    multVecMatrix(p.x(), p.y(), x, y);
    return FloatPoint(static_cast<float>(x), static_cast<float>(y));
}

FloatPoint3D TransformationMatrix::mapPoint(const FloatPoint3D& p) const
{
    double x, y, z;
    multVecMatrix(p.x(), p.y(), p.z(), x, y, z);
    return FloatPoint3D(static_cast<float>(x), static_cast<float>(y), static_cast<float>(z));
}

IntPoint TransformationMatrix::mapPoint(const IntPoint& point) const
{
    double x, y;
    multVecMatrix(point.x(), point.y(), x, y);

    // Round the point.
    return IntPoint(lround(x), lround(y));
}

IntRect TransformationMatrix::mapRect(const IntRect &rect) const
{
    return enclosingIntRect(mapRect(FloatRect(rect)));
}

FloatRect TransformationMatrix::mapRect(const FloatRect& r) const
{
    FloatQuad resultQuad = mapQuad(FloatQuad(r));
    return resultQuad.boundingBox();
}

FloatQuad TransformationMatrix::mapQuad(const FloatQuad& q) const
{
    FloatQuad result;
    result.setP1(mapPoint(q.p1()));
    result.setP2(mapPoint(q.p2()));
    result.setP3(mapPoint(q.p3()));
    result.setP4(mapPoint(q.p4()));
    return result;
}

TransformationMatrix& TransformationMatrix::scaleNonUniform(double sx, double sy)
{
    TransformationMatrix mat;
    mat.m_matrix[0][0] = sx;
    mat.m_matrix[1][1] = sy;

    multLeft(mat);
    return *this;
}

TransformationMatrix& TransformationMatrix::scale3d(double sx, double sy, double sz)
{
    TransformationMatrix mat;
    mat.m_matrix[0][0] = sx;
    mat.m_matrix[1][1] = sy;
    mat.m_matrix[2][2] = sz;

    multLeft(mat);
    return *this;
}

TransformationMatrix& TransformationMatrix::rotate3d(double x, double y, double z, double angle)
{
    // angles are in degrees. Switch to radians
    angle = deg2rad(angle);
    
    angle /= 2.0f;
    double sinA = sin(angle);
    double cosA = cos(angle);
    double sinA2 = sinA * sinA;
    
    // normalize
    double length = sqrt(x * x + y * y + z * z);
    if (length == 0) {
        // bad vector, just use something reasonable
        x = 0;
        y = 0;
        z = 1;
    } else if (length != 1) {
        x /= length;
        y /= length;
        z /= length;
    }
    
    TransformationMatrix mat;

    // optimize case where axis is along major axis
    if (x == 1.0f && y == 0.0f && z == 0.0f) {
        mat.m_matrix[0][0] = 1.0f;
        mat.m_matrix[0][1] = 0.0f;
        mat.m_matrix[0][2] = 0.0f;
        mat.m_matrix[1][0] = 0.0f;
        mat.m_matrix[1][1] = 1.0f - 2.0f * sinA2;
        mat.m_matrix[1][2] = 2.0f * sinA * cosA;
        mat.m_matrix[2][0] = 0.0f;
        mat.m_matrix[2][1] = -2.0f * sinA * cosA;
        mat.m_matrix[2][2] = 1.0f - 2.0f * sinA2;
        mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
        mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
        mat.m_matrix[3][3] = 1.0f;
    } else if (x == 0.0f && y == 1.0f && z == 0.0f) {
        mat.m_matrix[0][0] = 1.0f - 2.0f * sinA2;
        mat.m_matrix[0][1] = 0.0f;
        mat.m_matrix[0][2] = -2.0f * sinA * cosA;
        mat.m_matrix[1][0] = 0.0f;
        mat.m_matrix[1][1] = 1.0f;
        mat.m_matrix[1][2] = 0.0f;
        mat.m_matrix[2][0] = 2.0f * sinA * cosA;
        mat.m_matrix[2][1] = 0.0f;
        mat.m_matrix[2][2] = 1.0f - 2.0f * sinA2;
        mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
        mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
        mat.m_matrix[3][3] = 1.0f;
    } else if (x == 0.0f && y == 0.0f && z == 1.0f) {
        mat.m_matrix[0][0] = 1.0f - 2.0f * sinA2;
        mat.m_matrix[0][1] = 2.0f * sinA * cosA;
        mat.m_matrix[0][2] = 0.0f;
        mat.m_matrix[1][0] = -2.0f * sinA * cosA;
        mat.m_matrix[1][1] = 1.0f - 2.0f * sinA2;
        mat.m_matrix[1][2] = 0.0f;
        mat.m_matrix[2][0] = 0.0f;
        mat.m_matrix[2][1] = 0.0f;
        mat.m_matrix[2][2] = 1.0f;
        mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
        mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
        mat.m_matrix[3][3] = 1.0f;
    } else {
        double x2 = x*x;
        double y2 = y*y;
        double z2 = z*z;
    
        mat.m_matrix[0][0] = 1.0f - 2.0f * (y2 + z2) * sinA2;
        mat.m_matrix[0][1] = 2.0f * (x * y * sinA2 + z * sinA * cosA);
        mat.m_matrix[0][2] = 2.0f * (x * z * sinA2 - y * sinA * cosA);
        mat.m_matrix[1][0] = 2.0f * (y * x * sinA2 - z * sinA * cosA);
        mat.m_matrix[1][1] = 1.0f - 2.0f * (z2 + x2) * sinA2;
        mat.m_matrix[1][2] = 2.0f * (y * z * sinA2 + x * sinA * cosA);
        mat.m_matrix[2][0] = 2.0f * (z * x * sinA2 + y * sinA * cosA);
        mat.m_matrix[2][1] = 2.0f * (z * y * sinA2 - x * sinA * cosA);
        mat.m_matrix[2][2] = 1.0f - 2.0f * (x2 + y2) * sinA2;
        mat.m_matrix[0][3] = mat.m_matrix[1][3] = mat.m_matrix[2][3] = 0.0f;
        mat.m_matrix[3][0] = mat.m_matrix[3][1] = mat.m_matrix[3][2] = 0.0f;
        mat.m_matrix[3][3] = 1.0f;
    }
    multLeft(mat);
    return *this;
}