Math/MAT3geom.c

   1 /* #include "HEADERS.h" */
   2 /* Copyright 1988, Brown Computer Graphics Group.  All Rights Reserved. */
   3
   4 /* --------------------------------------------------------------------------
   5  * This file contains routines that perform geometry-related operations
   6  * on matrices.
   7  * -------------------------------------------------------------------------*/
   8
   9 #include "mat3defs.h"
  10
  11 /* --------------------------  Static Routines  ---------------------------- */
  12
  13 /* -------------------------  Internal Routines  --------------------------- */
  14
  15 /* --------------------------  Public Routines  ---------------------------- */
  16
  17 /*
  18  * This takes a matrix used to transform points, and returns a corresponding
  19  * matrix that can be used to transform direction vectors (between points).
  20  */
  21
  22 void
  23 MAT3direction_matrix(result_mat, mat)
  24 register MAT3mat result_mat, mat;
  25 {
  26    register int i;
  27
  28    MAT3copy(result_mat, mat);
  29
  30    for (i = 0; i < 4; i++) result_mat[i][3] = result_mat[3][i] = 0.0;
  31
  32    result_mat[3][3] = 1.0;
  33 }
  34
  35 /*
  36  * This takes a matrix used to transform points, and returns a corresponding
  37  * matrix that can be used to transform vectors that must remain perpendicular
  38  * to planes defined by the points.  It is useful when you are transforming
  39  * some object that has both points and normals in its definition, and you
  40  * only have the transformation matrix for the points.  This routine returns
  41  * FALSE if the normal matrix is uncomputable.  Otherwise, it returns TRUE.
  42  *
  43  * Spike sez: "This is the adjoint for the non-homogeneous part of the
  44  *             transformation."
  45  */
  46
  47 int
  48 MAT3normal_matrix(result_mat, mat)
  49 register MAT3mat result_mat, mat;
  50 {
  51    register int ret;
  52    MAT3mat      tmp_mat;
  53
  54    MAT3direction_matrix(result_mat, mat);
  55
  56    if ( (ret = MAT3invert(tmp_mat, tmp_mat)) ) {
  57        MAT3transpose(result_mat, tmp_mat);
  58    }
  59
  60    return(ret);
  61 }
  62
  63 /*
  64  * Sets the given matrix to be a scale matrix for the given vector of
  65  * scale values.
  66  */
  67
  68 void
  69 MAT3scale(result_mat, scale)
  70 MAT3mat result_mat;
  71 MAT3vec scale;
  72 {
  73    MAT3identity(result_mat);
  74
  75    result_mat[0][0] = scale[0];
  76    result_mat[1][1] = scale[1];
  77    result_mat[2][2] = scale[2];
  78 }
  79
  80 /*
  81  * Sets up a matrix for a rotation about an axis given by the line from
  82  * (0,0,0) to axis, through an angle (in radians).
  83  * Looking along the axis toward the origin, the rotation is counter-clockwise.
  84  */
  85
  86 #define SELECT  .7071   /* selection constant (roughly .5*sqrt(2) */
  87
  88 void
  89 MAT3rotate(result_mat, axis, angle_in_radians)
  90 MAT3mat result_mat;
  91 MAT3vec axis;
  92 double  angle_in_radians;
  93 {
  94    MAT3vec      naxis,  /* Axis of rotation, normalized                 */
  95                 base2,  /* 2nd unit basis vec, perp to axis             */
  96                 base3;  /* 3rd unit basis vec, perp to axis & base2     */
  97    double       dot;
  98    MAT3mat      base_mat,       /* Change-of-basis matrix               */
  99                 base_mat_trans; /* Inverse of c-o-b matrix              */
 100    register int i;
 101
 102    /* Step 1: extend { axis } to a basis for 3-space: { axis, base2, base3 }
 103     * which is orthonormal (all three have unit length, and all three are
 104     * mutually orthogonal). Also should be oriented, i.e. axis cross base2 =
 105     * base3, rather than -base3.
 106     *
 107     * Method: Find a vector linearly independent from axis. For this we
 108     * either use the y-axis, or, if that is too close to axis, the
 109     * z-axis. 'Too close' means that the dot product is too near to 1.
 110     */
 111
 112    MAT3_COPY_VEC(naxis, axis);
 113    MAT3_NORMALIZE_VEC(naxis, dot);
 114
 115    if (dot == 0.0) {
 116        /* ERR_ERROR(MAT3_errid, ERR_SEVERE,
 117                    (ERR_S, "Zero-length axis vector given to MAT3rotate")); */
 118       return;
 119    }
 120
 121    MAT3perp_vec(base2, naxis, TRUE);
 122    MAT3cross_product(base3, naxis, base2);
 123
 124    /* Set up the change-of-basis matrix, and its inverse */
 125    MAT3identity(base_mat);
 126    MAT3identity(base_mat_trans);
 127    MAT3identity(result_mat);
 128
 129    for (i = 0; i < 3; i++){
 130       base_mat_trans[i][0] = base_mat[0][i] = naxis[i];
 131       base_mat_trans[i][1] = base_mat[1][i] = base2[i];
 132       base_mat_trans[i][2] = base_mat[2][i] = base3[i];
 133    }
 134
 135    /* If T(u) = uR, where R is base_mat, then T(x-axis) = naxis,
 136     * T(y-axis) = base2, and T(z-axis) = base3. The inverse of base_mat is
 137     * its transpose.  OK?
 138     */
 139
 140    result_mat[1][1] =   result_mat[2][2] = cos(angle_in_radians);
 141    result_mat[2][1] = -(result_mat[1][2] = sin(angle_in_radians));
 142
 143    MAT3mult(result_mat, base_mat_trans, result_mat);
 144    MAT3mult(result_mat, result_mat,     base_mat);
 145 }
 146
 147 /*
 148  * Sets the given matrix to be a translation matrix for the given vector of
 149  * translation values.
 150  */
 151
 152 void
 153 MAT3translate(result_mat, trans)
 154 MAT3mat result_mat;
 155 MAT3vec trans;
 156 {
 157    MAT3identity(result_mat);
 158
 159    result_mat[3][0] = trans[0];
 160    result_mat[3][1] = trans[1];
 161    result_mat[3][2] = trans[2];
 162 }
 163
 164 /*
 165  * Sets the given matrix to be a shear matrix for the given x and y shear
 166  * values.
 167  */
 168
 169 void
 170 MAT3shear(result_mat, xshear, yshear)
 171 MAT3mat result_mat;
 172 double  xshear, yshear;
 173 {
 174    MAT3identity(result_mat);
 175
 176    result_mat[2][0] = xshear;
 177    result_mat[2][1] = yshear;
 178 }
 179