/* :ts=8 */ /************************************************************************/ /* */ /* This file contains the hairy details of rotating and computing points*/ /* on splines. */ /* */ /* This file is a good candidate for optimization! */ /* */ /************************************************************************/ #include #include #include #include #include "general.h" #include "globals.h" #include "spl_math.h" /* The hermite basis matrix (Computer Graphics p. 484) */ static Matrix4_T M_Hermite = { 2, -2, 1, 1, -3, 3, -2, -1, 0, 0, 1, 0, 1, 0, 0, 0 }; void Matrix_Mult(Matrix_T MRes, Matrix_T M1, Matrix_T M2) /************************************************************************/ /* */ /* ==== == == */ /* MRes = M1 * M2 */ /* */ /************************************************************************/ { int i, j; for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) MRes[i][j] = M1[i][0] * M2[0][j] + M1[i][1] * M2[1][j] + M1[i][2] * M2[2][j]; } /* Matrix_Mult */ void Matrix4_Mult(Matrix4_T MRes, Matrix4_T M1, Matrix4_T M2) /************************************************************************/ /* */ /* ==== == == */ /* MRes = M1 * M2 */ /* */ /************************************************************************/ { int i, j; for (i = 0; i < 4; i++) for (j = 0; j < 4; j++) MRes[i][j] = M1[i][0] * M2[0][j] + M1[i][1] * M2[1][j] + M1[i][2] * M2[2][j] + M1[i][3] * M2[3][j]; } /* Matrix4_Mult */ void Matrix_VMult(Vector_T VRes, Vector_T V, Matrix_T M) /************************************************************************/ /* */ /* ---- - = */ /* VRes = V * M */ /* */ /************************************************************************/ { int i; for (i = 0; i < 3; i++) VRes[i] = V[0] * M[0][i] + V[1] * M[1][i] + V[2] * M[2][i]; } /* Matrix_VMult */ void Matrix4_VMult(Vector4_T VRes, Vector4_T V, Matrix4_T M) /************************************************************************/ /* */ /* ---- - = */ /* VRes = V * M */ /* */ /************************************************************************/ { int i; for (i = 0; i < 4; i++) VRes[i] = V[0] * M[0][i] + V[1] * M[1][i] + V[2] * M[2][i] + V[3] * M[3][i]; } /* Matrix4_VMult */ void Matrix_MultV(Vector_T VRes, Matrix_T M, Vector_T V) /************************************************************************/ /* */ /* ---- = - */ /* VRes = M * V */ /* */ /************************************************************************/ { int i; for (i = 0; i < 3; i++) VRes[i] = V[0] * M[i][0] + V[1] * M[i][1] + V[2] * M[i][2]; } /* Matrix_MultV */ void Compute_Rotation_Matrix(Matrix_T R_Matrix, Vector_T Rot) /* Degrees */ /************************************************************************/ /* */ /* Precompute 3d rotation matrix for a rotation of 'Rot[0]' degrees */ /* around the X-axis, 'Rot[1]' degrees around the Y-axis and 'Rot[2]' */ /* degrees around the Z-axis. */ /* */ /************************************************************************/ { Matrix_T M_X, M_Y, M_Z, M_Tmp; double Rot_X, Rot_Y, Rot_Z; double S, C; Rot_X = PI * Rot[0] / 180.0; Rot_Y = PI * Rot[1] / 180.0; Rot_Z = PI * Rot[2] / 180.0; S = sin(Rot_X); C = cos(Rot_X); M_X[0][0] = 1.0; M_X[1][0] = 0.0; M_X[2][0] = 0.0; M_X[0][1] = 0.0; M_X[1][1] = C; M_X[2][1] = S; M_X[0][2] = 0.0; M_X[1][2] = -S; M_X[2][2] = C; S = sin(Rot_Y); C = cos(Rot_Y); M_Y[0][0] = C; M_Y[1][0] = 0.0; M_Y[2][0] = -S; M_Y[0][1] = 0.0; M_Y[1][1] = 1.0; M_Y[2][1] = 0.0; M_Y[0][2] = S; M_Y[1][2] = 0.0; M_Y[2][2] = C; S = sin(Rot_Z); C = cos(Rot_Z); M_Z[0][0] = C; M_Z[1][0] = S; M_Z[2][0] = 0.0; M_Z[0][1] = -S; M_Z[1][1] = C; M_Z[2][1] = 0.0; M_Z[0][2] = 0.0; M_Z[1][2] = 0.0; M_Z[2][2] = 1.0; Matrix_Mult(M_Tmp, M_X, M_Y); Matrix_Mult(R_Matrix, M_Tmp, M_Z); } /* Compute_Rotation_Matrix */ void Compute_Inverse_Rotation_Matrix(Matrix_T R_Matrix, Vector_T Rot) /************************************************************************/ /* */ /* Precompute inverse 3d rotation matrix for a rotation of 'Rot[0]' */ /* degrees around the X-axis, 'Rot[1]' degrees around the Y-axis and */ /* 'Rot[2]' degrees around the Z-axis. */ /* */ /************************************************************************/ { Matrix_T M_X, M_Y, M_Z, M_Tmp; double Rot_X, Rot_Y, Rot_Z; double S, C; Rot_X = PI * Rot[0] / 180.0; Rot_Y = PI * Rot[1] / 180.0; Rot_Z = PI * Rot[2] / 180.0; S = sin(-Rot_X); C = cos(-Rot_X); M_X[0][0] = 1.0; M_X[1][0] = 0.0; M_X[2][0] = 0.0; M_X[0][1] = 0.0; M_X[1][1] = C; M_X[2][1] = S; M_X[0][2] = 0.0; M_X[1][2] = -S; M_X[2][2] = C; S = sin(-Rot_Y); C = cos(-Rot_Y); M_Y[0][0] = C; M_Y[1][0] = 0.0; M_Y[2][0] = -S; M_Y[0][1] = 0.0; M_Y[1][1] = 1.0; M_Y[2][1] = 0.0; M_Y[0][2] = S; M_Y[1][2] = 0.0; M_Y[2][2] = C; S = sin(-Rot_Z); C = cos(-Rot_Z); M_Z[0][0] = C; M_Z[1][0] = S; M_Z[2][0] = 0.0; M_Z[0][1] = -S; M_Z[1][1] = C; M_Z[2][1] = 0.0; M_Z[0][2] = 0.0; M_Z[1][2] = 0.0; M_Z[2][2] = 1.0; Matrix_Mult(M_Tmp, M_Z, M_Y); Matrix_Mult(R_Matrix, M_Tmp, M_X); } /* Compute_Inverse_Rotation_Matrix */ void Precompute_Kochanek_Bartels(Matrix4_T Mres, Spline_T *Spline, Knot_T *Knot) /************************************************************************/ /* */ /* Precompute matrix used for Kochanek-Bartels spline interpolation. */ /* See Image Synthesis p. 36 - 38. */ /* */ /************************************************************************/ { Matrix4_T M_Tmp; double t, b, c, tm, bm, cm, bp, cp; Knot_T *Knot0, *Knot1, *Knot2, *Knot3; short PId0, PId1, PId2, PId3; short Nbr_Knots; Nbr_Knots = Spline->Nbr_Knots; /* The knot before this knot */ Knot0 = Knot->Previous; if (Knot0 == NULL) Knot0 = Spline->Last; /* Loop */ /* This knot */ Knot1 = Knot; /* The knot after this knot */ Knot2 = Knot->Next; if (Knot2 == NULL) Knot2 = Spline->First; /* Loop */ /* The second knot after this knot */ Knot3 = Knot2->Next; if (Knot3 == NULL) Knot3 = Spline->First; /* Loop */ PId0 = Knot0->Point_Id; PId1 = Knot1->Point_Id; PId2 = Knot2->Point_Id; PId3 = Knot3->Point_Id; M_Tmp[0][0] = Points[PId1].Pos[0]; M_Tmp[0][1] = Points[PId1].Pos[1]; M_Tmp[0][2] = Points[PId1].Pos[2]; M_Tmp[0][3] = 0.0; M_Tmp[1][0] = Points[PId2].Pos[0]; M_Tmp[1][1] = Points[PId2].Pos[1]; M_Tmp[1][2] = Points[PId2].Pos[2]; M_Tmp[1][3] = 0.0; t = Knot1->Tension; tm = 1.0 - t; b = Knot1->Bias; bm = 1.0 - b; bp = 1.0 + b; c = Knot1->Continuity; cm = 1.0 - c; cp = 1.0 + c; if (Knot != Spline->First || Spline->Loop) { M_Tmp[2][0] = 0.5 * tm * ( bm * cm * (Points[PId2].Pos[0] - Points[PId1].Pos[0]) + bp * cp * (Points[PId1].Pos[0] - Points[PId0].Pos[0]) ); M_Tmp[2][1] = 0.5 * tm * ( bm * cm * (Points[PId2].Pos[1] - Points[PId1].Pos[1]) + bp * cp * (Points[PId1].Pos[1] - Points[PId0].Pos[1]) ); M_Tmp[2][2] = 0.5 * tm * ( bm * cm * (Points[PId2].Pos[2] - Points[PId1].Pos[2]) + bp * cp * (Points[PId1].Pos[2] - Points[PId0].Pos[2]) ); M_Tmp[2][3] = 0.0; } else { M_Tmp[2][0] = 0.0; M_Tmp[2][1] = 0.0; M_Tmp[2][2] = 0.0; M_Tmp[2][3] = 0.0; } /* if .. else .. */ t = Knot2->Tension; tm = 1.0 - t; b = Knot2->Bias; bm = 1.0 - b; bp = 1.0 + b; c = Knot2->Continuity; cm = 1.0 - c; cp = 1.0 + c; if (Knot->Next != NULL && Knot->Next->Next != NULL) { M_Tmp[3][0] = 0.5 * tm * ( bm * cp * (Points[PId3].Pos[0] - Points[PId2].Pos[0]) + bp * cm * (Points[PId2].Pos[0] - Points[PId1].Pos[0]) ); M_Tmp[3][1] = 0.5 * tm * ( bm * cp * (Points[PId3].Pos[1] - Points[PId2].Pos[1]) + bp * cm * (Points[PId2].Pos[1] - Points[PId1].Pos[1]) ); M_Tmp[3][2] = 0.5 * tm * ( bm * cp * (Points[PId3].Pos[2] - Points[PId2].Pos[2]) + bp * cm * (Points[PId2].Pos[2] - Points[PId1].Pos[2]) ); M_Tmp[3][3] = 0.0; } else { M_Tmp[3][0] = 0.0; M_Tmp[3][1] = 0.0; M_Tmp[3][2] = 0.0; M_Tmp[3][3] = 0.0; } /* if .. else .. */ Matrix4_Mult(Mres, M_Hermite, M_Tmp); } /* Precompute_Kochanek_Bartels */ void Compute_Position(Vector4_T Vres, double T, Matrix4_T M_Precomp) /************************************************************************/ /* */ /* Compute position on a Kochanek Bartels spline at parameter value T */ /* using the precomputed matrix M_Precomp. */ /* Return the result in Vres */ /* */ /************************************************************************/ { Vector4_T Vt; Vt[3] = 1.0; Vt[2] = T; Vt[1] = T * T; Vt[0] = Vt[1] * T; Matrix4_VMult(Vres, Vt, M_Precomp); } /* Compute_Position */