home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / gems / g_gemsi.lha / GraphicsGems / Forms.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-04-10  |  5.2 KB  |  210 lines
  1. /*
  2. Forms, Vectors,and Transforms
  3. by Bob Wallis
  4. from "Graphics Gems", Academic Press, 1990
  5. */
  6.  
  7. /*----------------------------------------------------------------------
  8. The main program below is set up to solve the Bezier subdivision problem
  9. in "Forms, Vectors, and Transforms".  The subroutines are useful in
  10. solving general problems which require manipulating matrices via exact
  11. integer arithmetic.  The intended application is validating or avoiding
  12. tedious algebraic calculations.  As such, no thought was given to
  13. efficiency.  
  14. ----------------------------------------------------------------------*/
  15. #define ABS(x) ((x)>(0)? (x):(-x))
  16. #define N 4            /* size of matrices to deal with */
  17. int     M[N][N] =        /* Bezier weights */
  18. {
  19.       1,   0,   0,   0,
  20.      -3,   3,   0,   0,
  21.       3,  -6,   3,   0,
  22.      -1,   3,  -3,   1,
  23. };
  24. int     T[N][N] =    /* re-parameterization xform for top half */
  25. {
  26.       1,  -1,   1,  -1,
  27.       0,   2,  -4,   6,
  28.       0,   0,   4, -12,
  29.       0,   0,   0,   8
  30. };
  31. main ()
  32. {
  33.     int     i,
  34.             j,
  35.             scale,
  36.             gcd,
  37.             C[N][N],
  38.             S[N][N],
  39.             Madj[N][N],
  40.             Tadj[N][N],
  41.             Mdet,
  42.             Tdet;
  43.  
  44.  
  45.     Tdet = adjoint (T, Tadj);   /* inverse without division by */
  46.     Mdet = adjoint (M, Madj);   /* determinant of T and M */
  47.     matmult (Madj, Tadj, C);
  48.     matmult (C, M, S);        /* Madj*Tadj*M -> S */
  49.     scale = gcd = Mdet * Tdet;  /* scale factors of both determinants */
  50.     for (i = 0; i < N; i++)    /* find the greatest common */
  51.     {                /* demoninator of S and determinants */
  52.         for (j = 0; j < N; j++)
  53.             gcd = Gcd (gcd, S[i][j]);
  54.     }
  55.     scale /= gcd;        /* divide everything by gcd to get */
  56.     for (i = 0; i < N; i++)    /* matrix and scale factor in lowest */
  57.     {                /* integer terms possible */
  58.         for (j = 0; j < N; j++)
  59.             S[i][j] /= gcd;
  60.     }
  61.     printf ("scale factor = 1/%d  ", scale);
  62.     print_mat ("M=", M, N);     /* display the results */
  63.     print_mat ("T=", T, N);
  64.     print_mat ("S=", S, N);     /* subdivision matrix */
  65.     exit (0);
  66. }
  67. Gcd (a, b)            /*returns greatest common demoninator */
  68. int     a,            /* of (a,b) */
  69.            b;
  70. {
  71.     int        i,
  72.                 r;
  73.     a = ABS (a);        /* force positive */
  74.     b = ABS (b);
  75.     if (a < b)            /* exchange so that a >= b */
  76.     {
  77.         i = b;
  78.         b = a;
  79.         a = i;
  80.     }
  81.     if     (b == 0)
  82.         return (a);    /* finished */
  83.     r = a % b;            /* remainder */
  84.     if (r == 0)
  85.         return (b);    /* finished */
  86.     else            /* recursive call */
  87.         return (Gcd (b, r));
  88. }
  89.  
  90. adjoint (A, Aadj)            /* returns determinant of A */
  91. int     A[N][N],            /* input matrix */
  92.         Aadj[N][N];            /* output = adjoint of A */
  93. {                    /* must have N >= 3 */
  94.     int       i,
  95.             j,
  96.             I[N],            /* arrays of row and column indices */
  97.             J[N],
  98.             Isub[N],            /* sub-arrays of the above */
  99.             Jsub[N],
  100.             cofactor,
  101.             det;
  102.     if (N < 3)
  103.     {
  104.         printf ("must have N >= 3\n");
  105.         exit (1);
  106.     }
  107.     for (i = 0; i < N; i++)
  108.     {                    /* lookup tables to select a */
  109.         I[i] = i;        /* particular subset of */
  110.         J[i] = i;        /* rows and columns */
  111.     }
  112.     det = 0;
  113.     for (i = 0; i < N; i++)
  114.     {                    /* delete ith row */
  115.         subarray (I, Isub, N, i);
  116.         for (j = 0; j < N; j++)
  117.         {                /* delete jth column */
  118.                 subarray (J, Jsub, N, j);
  119.                 cofactor = determinant (A, Isub, Jsub, N - 1,(i + j) & 1);
  120.                 if (j == 0)        /* use 0th column for det */
  121.                 det += cofactor * A[i][0];
  122.             Aadj[j][i] = cofactor;
  123.         }
  124.        }
  125.     return (det);
  126. }
  127. determinant (A, I, J, n, parity)/* actually gets a sub-determinant */
  128. int     A[N][N],            /* input = entire matrix */
  129.            I[N],                /* row sub-array we want */
  130.         J[N],                /* col sub-array we want */
  131.         parity,                /* 1-> flip polarity */
  132.         n;                /* # elements in subarrays */
  133.  
  134. {
  135.     int     i,
  136.             j,
  137.             det,
  138.             j_,
  139.             Jsub[N];
  140.     if (n <= 2)            /* call ourselves till we get down to */
  141.     {                /* a 2x2 matrix */
  142.         det =
  143.             (A[I[0]][J[0]] * A[I[1]][J[1]]) -
  144.             (A[I[1]][J[0]] * A[I[0]][J[1]]);
  145.         if (parity)
  146.         det = -det;
  147.         return (det);
  148.     }                    /* if (n <= 2) */
  149.     det = 0;                /* n > 2; call recursively */
  150.     i = I[0];                /* strike out 0th row */
  151.     for (j_ = 0; j_ < n; j_++)
  152.     {                    /* strike out jth column */
  153.         subarray (J, Jsub, n, j_);
  154.         j = J[j_];        /* I + 1 => struck out 0th row */
  155.         det += A[i][j] * determinant (A, I + 1, Jsub, n - 1, j_ & 1);
  156.     }
  157.     if (parity)
  158.         det = -det;
  159.     return (det);
  160. }
  161. subarray (src, dest, n, k)        /* strike out kth row/column */
  162. int      *src,                /*     source array of n indices */
  163.          *dest,                /* dest array formed by deleting k  */
  164.             n,
  165.             k;
  166. {
  167.     int     i;
  168.     for (i = 0; i < n; i++, src++)
  169.         if (i != k)            /* skip over k */
  170.             *dest++ = *src;
  171. }
  172.  
  173. matmult (A, B, C)                /* C = A*B */
  174. int     A[N][N],
  175.         B[N][N],
  176.         C[N][N];
  177. {
  178.     int        i,
  179.         j,
  180.         k,
  181.         sum;
  182.     for (i = 0; i < N; i++)
  183.     {
  184.         for (k = 0; k < N; k++)
  185.         {
  186.             sum = 0;
  187.             for (j = 0; j < N; j++)
  188.                 sum += A[i][j] * B[j][k];
  189.             C[i][k] = sum;
  190.         }
  191.     }
  192. }
  193. print_mat (string, mat, n)
  194. char   *string;
  195. int     mat[N][N],
  196.          n;
  197. {
  198.     int     i,
  199.              j;
  200.     printf ("%s\n", string);
  201.     for (i = 0; i < n; i++)
  202.     {
  203.         for (j = 0; j < n; j++)
  204.         printf (" %8ld", mat[i][j]);
  205.     printf ("\n");
  206.     }
  207. }
  208.  
  209.  
  210.