我帮你梳理了代码里的几个关键问题，这些问题导致旋转矩阵无法正确将[1,1,1]映射到目标向量，另外还要先明确一个核心前提：旋转矩阵是正交矩阵，会严格保持向量长度不变。所以如果你的目标向量[x,y,z]的长度不等于[1,1,1]的长度（即√3 ≈ 1.732），严格来说不存在这样的纯旋转矩阵——你要么调整目标向量使其长度匹配，要么改用包含缩放的刚体变换矩阵。

具体错误点与修正方式 #

1. 归一化逻辑误用
   你在函数开头直接对目标向量和[1,1,1]做了归一化，但需求是让旋转矩阵作用在原始的[1,1,1]上得到原始的目标向量。如果目标向量和[1,1,1]长度一致，归一化计算旋转矩阵是没问题的，但验证时应该用归一化后的[1,1,1]乘矩阵再缩放回原长度；如果长度不一致，纯旋转矩阵根本无法满足需求。

2. c的计算完全错误
   Rodrigues旋转公式里的c是两个归一化向量的点积（也就是cosθ，θ为两向量夹角），但你的代码写成了：

   float c = dot(vid,vector)*cos(alfa);
   

   因为dot(vid, vector)本身就是cos(alfa)（向量已归一化），这会让c变成cos²(alfa)，完全偏离公式要求。正确写法是：

   float c = dot(vid, vector);
   

3. s的计算错误
   公式里的s是叉积的长度（对应sinθ），但你的代码写成了：

   float s = len(v)*sin(alfa);
   

   而v是两个归一化向量的叉积，len(v)本身就是sin(alfa)，这会得到sin²(alfa)。正确写法是：

   float s = len(v);
   

   或者直接用sin(alfa)，因为alfa已经是两向量的夹角。

4. 向量旋转的矩阵乘法索引错误
   vector_rotation函数里的矩阵乘法逻辑完全错了，正确的旋转矩阵乘向量应该是对结果的每个元素i，计算矩阵第i行与输入向量的点积。修正后的函数：

   void vector_rotation(float rotation_matrix[3][3], float vector[3], float result[3]) {
       int i, k;
       for (i = 0; i < 3; i++) {
           result[i] = 0;
           for (k = 0; k < 3; k++) {
               result[i] += rotation_matrix[i][k] * vector[k];
           }
       }
   }
   

5. len与norm函数混用
   代码里定义了len函数计算向量长度，但打印时用了未定义的norm函数，这会导致编译错误，统一替换为len即可。

6. 边界情况未处理
   当两向量几乎平行时（s接近0），s*s会趋近于0，导致k=(1-c)/(s*s)出现除以0的问题。需要添加判断：如果s小于一个极小值（比如1e-6），说明两向量几乎同向或反向，此时旋转矩阵要么是单位矩阵，要么是绕垂直轴旋转180度的矩阵。

修正后的完整代码 #

#include <stdio.h>
#include <math.h>
#define pi 3.14159265359
typedef struct float3_ { float x,y,z; } float3;

float len(float* v) { 
    return sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]); 
}

float dot(float* v1,float* v2) { 
    return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]; 
}

void cross(float* v1,float* v2, float* v3) { 
    v3[0] = v1[1]*v2[2] - v2[1]*v1[2]; 
    v3[1] = v1[2]*v2[0] - v2[2]*v1[0]; 
    v3[2] = v1[0]*v2[1] - v2[0]*v1[1]; 
}

void normal(float* v) { 
    float temp = len(v); 
    if (temp < 1e-6) temp = 1; // 避免除以0
    v[0]/=temp; 
    v[1]/=temp; 
    v[2]/=temp; 
}

void matrix_product(float a[3][3], float b[3][3], float c[3][3]) {
    int i,j,k;
    for(i = 0; i < 3; i++ ) {
        for(j = 0; j < 3; j++ ) {
            c[i][j] = 0;
            for(k = 0; k < 3; k++ ) {
                c[i][j] += a[i][k] * b[k][j];
            }
        }
    }
}

void prod_matrix_scalar(float m[3][3], float k, float r[3][3]) {
    int i,j;
    for(i=0;i<3;i++) {
        for(j=0;j<3;j++) {
            r[i][j] = m[i][j] * k;
        }
    }
}

void sum_matrix(float m[3][3], float n[3][3], float r[3][3]) {
    int i,j;
    for(i=0;i<3;i++) {
        for(j=0;j<3;j++) {
            r[i][j] = m[i][j] + n[i][j];
        }
    }
}

void print_matrix(float matrix[3][3]) { 
    int i,j; 
    for(i=0; i<3; i++) { 
        printf("\n"); 
        for(j=0; j<3; j++) {
            printf("%f ", matrix[i][j]); 
        }
    } 
}

void print_vector(float v[3]) { 
    int i=0; 
    printf("["); 
    for(i=0; i<3; i++) {
        printf(" %.3f ",v[i]); 
    }
    printf("]\n"); 
}

void vector_rotation(float rotation_matrix[3][3], float vector[3], float result[3]) {
    int i, k;
    for(i = 0; i < 3; i++) {
        result[i] = 0;
        for(k = 0; k < 3; k++) {
            result[i] += rotation_matrix[i][k] * vector[k];
        }
    }
}

void vector2rotation_matrix (float target[3], float R[3][3]) {
    float vid[3] = {1,1,1};
    float vid_norm[3] = {1,1,1};
    float target_norm[3];
    // 保存原始长度
    float len_vid = len(vid);
    float len_target = len(target);
    
    // 归一化两个向量用于计算旋转矩阵
    normal(vid_norm);
    for(int i=0;i<3;i++) target_norm[i] = target[i];
    normal(target_norm);
    
    printf("Original target vector:"); print_vector(target);
    printf("Original [1,1,1] vector:"); print_vector(vid);
    printf("\nNormalized target:"); print_vector(target_norm);
    printf("Normalized [1,1,1]:"); print_vector(vid_norm);
    
    float v[3];
    cross(vid_norm, target_norm, v);
    float s = len(v);
    float c = dot(vid_norm, target_norm);
    
    // 处理边界情况：两向量几乎平行
    if(s < 1e-6) {
        if(c > 0) {
            // 同向，单位矩阵
            for(int i=0;i<3;i++) {
                for(int j=0;j<3;j++) {
                    R[i][j] = (i==j) ? 1 : 0;
                }
            }
        } else {
            // 反向，绕垂直轴旋转180度
            // 找一个垂直于vid_norm的向量
            float u[3];
            if(fabs(vid_norm[0]) > fabs(vid_norm[1])) {
                u[0] = -vid_norm[2];
                u[1] = 0;
                u[2] = vid_norm[0];
            } else {
                u[0] = 0;
                u[1] = -vid_norm[2];
                u[2] = vid_norm[1];
            }
            normal(u);
            float ux[3][3] = {0, -u[2], u[1], u[2], 0, -u[0], -u[1], u[0], 0};
            float ux2[3][3];
            matrix_product(ux, ux, ux2);
            // R = I + 2*ux² (Rodrigues公式当theta=pi时)
            float Id[3][3] = {1,0,0,0,1,0,0,0,1};
            float ux2_2[3][3];
            prod_matrix_scalar(ux2, 2, ux2_2);
            sum_matrix(Id, ux2_2, R);
        }
        printf("\nVectors are nearly parallel, using special case rotation matrix.\n");
        return;
    }
    
    float alfa = acos(c);
    printf("\nAngle between vectors: %.3f rad (%.3f degrees)\n", alfa, alfa*180/pi);
    
    float Id[3][3] = {1,0,0,0,1,0,0,0,1};
    float Vx[3][3] = {0, -v[2], v[1], v[2], 0, -v[0], -v[1], v[0], 0};
    float Vx2[3][3];
    float Vx2k[3][3];
    float Ris[3][3];
    
    matrix_product(Vx, Vx, Vx2);
    float k = (1 - c)/(s*s);
    prod_matrix_scalar(Vx2, k, Vx2k);
    
    // 应用Rodrigues公式：R = I + Vx + k*Vx²
    sum_matrix(Id, Vx, Ris);
    sum_matrix(Ris, Vx2k, R);
    
    // 注意：如果目标向量和[1,1,1]长度不同，需要缩放矩阵才能满足R*[1,1,1] = target
    // 这里如果需要严格满足等式，且长度不同，要乘以缩放因子 len_target / len_vid
    // 取消下面的注释即可：
    // float scale = len_target / len_vid;
    // prod_matrix_scalar(R, scale, R);
}

int main() {
    // 测试用目标向量，这里用和[1,1,1]等长的向量（长度√3）
    // float v[3] = {1,1,-1}; // 长度√3，纯旋转可满足
    float v[3] = {2,4,5}; // 长度≠√3，需取消scale注释才能匹配
    float matrix[3][3];
    vector2rotation_matrix(v, matrix);
    
    printf("\n\nFinal rotation matrix:\n");
    print_matrix(matrix);
    
    float test_vec[3] = {1,1,1};
    float result[3];
    vector_rotation(matrix, test_vec, result);
    
    printf("\n\nR * [1,1,1] = ");
    print_vector(result);
    printf("Expected target vector: ");
    print_vector(v);
    
    return 0;
}

补充说明 #

• 如果你的目标向量和[1,1,1]长度不同，且必须满足R*[1,1,1] = 目标向量，那你需要的不是纯旋转矩阵，而是旋转+缩放的矩阵。此时取消代码中scale相关的注释即可，缩放因子为目标向量长度除以[1,1,1]的长度。
• 修正后的代码处理了向量平行的边界情况，避免了除以0的错误。

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