如何实现绕任意轴的缩放矩阵？射线追踪厚射线应用需求

核心原理 #

绕任意单位轴$\vec{v}$进行缩放，缩放因子为$s$，本质是对垂直于$\vec{v}$的所有方向做$s$倍缩放，$\vec{v}$自身方向保持不变。标准实现思路有两种：

• 直接通过向量运算构造矩阵元素
• 先将轴对齐到标准坐标轴（如Z轴），完成缩放后再转回到原轴方向

你的代码问题分析 #

你的思路方向正确，但问题出在矩阵乘法顺序和旋转矩阵的逆处理上：

• 多数矩阵库中，Multiply(B)是将当前矩阵左乘B（即this = this * B），而变换的应用顺序是从右到左，你当前的矩阵组合顺序完全颠倒
• 正交旋转矩阵的逆等于其转置，无需重新构造反向旋转矩阵，直接转置即可，否则可能引入额外误差

正确实现方案 #

方案1：直接构造矩阵（推荐，无旋转依赖） #

基于线性代数公式直接生成缩放矩阵，无需依赖旋转矩阵实现，效率更高：

#include <cmath>

// 基础向量结构
struct Vector {
    float x, y, z;
    Vector(float x=0, float y=0, float z=0) : x(x), y(y), z(z) {}
    
    float Dot(const Vector& other) const {
        return x*other.x + y*other.y + z*other.z;
    }
    
    Vector Normalized() const {
        float len = sqrt(x*x + y*y + z*z);
        return Vector(x/len, y/len, z/len);
    }
};

// 4x4矩阵结构
struct Matrix {
    float m[4][4];
    
    Matrix() {
        // 初始化单位矩阵
        for (int i = 0; i < 4; i++) {
            for (int j = 0; j < 4; j++) {
                m[i][j] = (i == j) ? 1.0f : 0.0f;
            }
        }
    }

    // 矩阵乘法：this = this * other
    Matrix& Multiply(const Matrix& other) {
        Matrix temp;
        for (int i = 0; i < 4; i++) {
            for (int j = 0; j < 4; j++) {
                temp.m[i][j] = 0.0f;
                for (int k = 0; k < 4; k++) {
                    temp.m[i][j] += m[i][k] * other.m[k][j];
                }
            }
        }
        *this = temp;
        return *this;
    }

    // 绕任意轴构造缩放矩阵
    static Matrix ScaleAroundAxis(const Vector& axis, float scale) {
        Vector v = axis.Normalized();
        float vx = v.x, vy = v.y, vz = v.z;
        float s = scale;
        float oneMinusS = 1.0f - s;

        Matrix mat;
        mat.m[0][0] = oneMinusS * vx*vx + s;
        mat.m[0][1] = oneMinusS * vx*vy;
        mat.m[0][2] = oneMinusS * vx*vz;

        mat.m[1][0] = oneMinusS * vy*vx;
        mat.m[1][1] = oneMinusS * vy*vy + s;
        mat.m[1][2] = oneMinusS * vy*vz;

        mat.m[2][0] = oneMinusS * vz*vx;
        mat.m[2][1] = oneMinusS * vz*vy;
        mat.m[2][2] = oneMinusS * vz*vz + s;

        // 第四行第四列保持1，其余为0
        return mat;
    }
};

方案2：旋转+缩放+逆旋转（修正你的思路） #

如果需要基于现有旋转矩阵实现，需严格控制乘法顺序，并使用转置作为逆矩阵：

// 补充Matrix类的旋转矩阵构造和转置方法
Matrix Matrix::BuildRotateFromTo(const Vector& from, const Vector& to) {
    Vector f = from.Normalized();
    Vector t = to.Normalized();
    Vector axis = Vector(
        f.y*t.z - f.z*t.y,
        f.z*t.x - f.x*t.z,
        f.x*t.y - f.y*t.x
    );
    float dot = f.Dot(t);
    float sinTheta = sqrt(axis.Dot(axis));
    
    // 处理共线情况
    if (sinTheta < 1e-6) {
        return dot > 0 ? Matrix() : Matrix(); // 共线同向返回单位矩阵，反向返回负单位矩阵（简化处理）
    }

    axis = axis.Normalized();
    float tx = axis.x, ty = axis.y, tz = axis.z;
    float cosTheta = dot;
    float oneMinusCos = 1.0f - cosTheta;

    Matrix mat;
    mat.m[0][0] = cosTheta + tx*tx*oneMinusCos;
    mat.m[0][1] = tx*ty*oneMinusCos - tz*sinTheta;
    mat.m[0][2] = tx*tz*oneMinusCos + ty*sinTheta;

    mat.m[1][0] = ty*tx*oneMinusCos + tz*sinTheta;
    mat.m[1][1] = cosTheta + ty*ty*oneMinusCos;
    mat.m[1][2] = ty*tz*oneMinusCos - tx*sinTheta;

    mat.m[2][0] = tz*tx*oneMinusCos - ty*sinTheta;
    mat.m[2][1] = tz*ty*oneMinusCos + tx*sinTheta;
    mat.m[2][2] = cosTheta + tz*tz*oneMinusCos;
    return mat;
}

Matrix Matrix::Transposed() const {
    Matrix mat;
    for (int i = 0; i < 4; i++) {
        for (int j = 0; j < 4; j++) {
            mat.m[i][j] = m[j][i];
        }
    }
    return mat;
}

// 修正后的绕轴缩放方法
Matrix& Matrix::ScaleAroundAxis3D(const Vector& theAxis, float theScale) {
    Vector axis = theAxis.Normalized();
    Matrix R = BuildRotateFromTo(axis, Vector(0, 0, 1)); // 将轴对齐到Z轴
    Matrix S;
    S.m[0][0] = theScale;
    S.m[1][1] = theScale;
    S.m[2][2] = 1.0f; // Z轴保持原长度

    Matrix R_inv = R.Transposed(); // 正交矩阵的逆等于转置
    Matrix result = R_inv;
    result.Multiply(S);
    result.Multiply(R); // 矩阵顺序：R_inv * S * R，对应变换顺序：R → S → R_inv

    return Multiply(result);
}

验证说明 #

• 对于你测试的Vector(0.707, 0, 0.707)轴，两种方案生成的矩阵都会保持所有点在原平面内，仅沿垂直于轴的方向缩放，符合预期
• 所有输入轴必须先归一化，否则矩阵计算会出现偏差

内容的提问来源于stack exchange，提问作者KiraHoneybee