home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / modelers / geomview / source.lha / Geomview / src / lib / geometry / hpoint3 / hg4.c < prev    next >
Encoding:
Text File  |  1992-02-26  |  8.9 KB  |  486 lines
  1. /* Copyright (c) 1992 The Geometry Center; University of Minnesota
  2.    1300 South Second Street;  Minneapolis, MN  55454, USA;
  3.    
  4. This file is part of geomview/OOGL. geomview/OOGL is free software;
  5. you can redistribute it and/or modify it only under the terms given in
  6. the file COPYING, which you should have received along with this file.
  7. This and other related software may be obtained via anonymous ftp from
  8. geom.umn.edu; email: software@geom.umn.edu. */
  9.  
  10. /* Authors: Charlie Gunn, Pat Hanrahan, Stuart Levy, Tamara Munzner, Mark Phillips */
  11.  
  12. #
  13. /*
  14. **    hg4.h - procedural interface to homogeneous geometry
  15. **
  16. **    pat hanrahan
  17. **
  18. */
  19.  
  20. # include <math.h>
  21. # include <stdio.h>
  22. # include "tolerance.h"
  23. # include "transform3.h"
  24. # include "hg4.h"
  25.  
  26. char *
  27. Hg4Create()
  28. {
  29.     return (char *) malloc( sizeof(Hg4Tensor1) );
  30. }
  31.  
  32. void
  33. Hg4Delete( p )
  34.     Hg4Tensor1 p;
  35. {
  36.     free( (char *) p );
  37. }
  38.  
  39. void
  40. Hg4Print( p )
  41.   Hg4Tensor1 p;
  42. {
  43.     if( p ) 
  44.     printf( "%g %g %g %g\n", p[X], p[Y], p[Z], p[W] );
  45. }
  46.  
  47. void
  48. Hg4From( p, x, y, z, w )
  49.     Hg4Tensor1 p;
  50.     Hg4Coord x, y, z, w;
  51. {
  52.     p[X] = x;
  53.     p[Y] = y;
  54.     p[Z] = z;
  55.     p[W] = w;
  56. }
  57.  
  58. void
  59. Hg4Copy( a, b )
  60.     Hg4Tensor1 a, b;
  61. {
  62.     bcopy( (char *)a, (char *)b, sizeof(Hg4Tensor1) );
  63. }
  64.  
  65. void 
  66. Hg4Add( p1, p2, p3)
  67.     Hg4Tensor1 p1, p2, p3;
  68. {
  69.     register int i;
  70.     for (i=0; i<4; ++i)
  71.     p3[i] = p1[i] + p2[i];
  72. } 
  73.  
  74. int
  75. Hg4Compare( p1, p2 )
  76.     Hg4Tensor1 p1, p2;
  77. {  
  78.     Hg4Coord test;
  79.  
  80.     test = p1[X]*p2[Y] - p1[Y]*p2[X];
  81.     if( fneg(test) ) return -1;
  82.     if( fpos(test) ) return  1;
  83.  
  84.     test = p1[X]*p2[Z] - p1[Z]*p2[X];
  85.     if( fneg(test) ) return -1;
  86.     if( fpos(test) ) return  1;
  87.  
  88.     test = p1[Y]*p2[Z] - p1[Z]*p2[Y];
  89.     if( fneg(test) ) return -1;
  90.     if( fpos(test) ) return  1;
  91.  
  92.     test = p1[X]*p2[W] - p1[W]*p2[X];
  93.     if( fneg(test) ) return -1;
  94.     if( fpos(test) ) return  1;
  95.  
  96.     test = p1[Y]*p2[W] - p1[W]*p2[Y];
  97.     if( fneg(test) ) return -1;
  98.     if( fpos(test) ) return  1;
  99.  
  100.     test = p1[Z]*p2[W] - p1[W]*p2[Z];
  101.     if( fneg(test) ) return -1;
  102.     if( fpos(test) ) return  1;
  103.  
  104.     return 0;
  105. }
  106.  
  107.  
  108. int
  109. Hg4Coincident( p1, p2 )
  110.     Hg4Tensor1 p1;
  111.     Hg4Tensor1 p2;
  112. {
  113.     return Hg4Compare( p1, p2 ) == 0;
  114. }
  115.  
  116. int
  117. Hg4Undefined( a )
  118.     Hg4Tensor1 a;
  119. {
  120.     if( !fzero(a[X]) ) return 0;
  121.     if( !fzero(a[Y]) ) return 0;
  122.     if( !fzero(a[Z]) ) return 0;
  123.     if( !fzero(a[W]) ) return 0;
  124.     return 1;
  125. }
  126.  
  127.  
  128. int
  129. Hg4Infinity( p, dual )
  130.     Hg4Tensor1 p;
  131.     int dual;
  132. {
  133.     /* Assume not undefined */
  134.     if( dual ) { /* plane */
  135.     if( !fzero(p[X]) ) return 0;
  136.     if( !fzero(p[Y]) ) return 0;
  137.     if( !fzero(p[Z]) ) return 0;
  138.     return 1;
  139.     }
  140.     else { /* point */
  141.     if( !fzero(p[W]) ) return 0;
  142.     return 1;
  143.     }
  144. }
  145.  
  146. void
  147. Hg4Normalize( p, q )
  148.     Hg4Tensor1 p, q;
  149. {
  150.     Hg4Copy( p, q );
  151.     if( q[W] != 1. && q[W] != 0. ) {
  152.     q[X] /= q[W];
  153.     q[Y] /= q[W];
  154.     q[Z] /= q[W];
  155.     q[W]  = 1.;
  156.     }
  157. }
  158.  
  159. void
  160. Hg4Pencil( t1, p1, t2, p2, p )
  161.     Hg4Coord t1, t2;
  162.     Hg4Tensor1 p1, p2, p;
  163. {  
  164.     p[W] = t1 * p1[W] + t2 * p2[W];
  165.     /* Keep W positive */
  166.     if( p[W] < 0. ) {
  167.     p[W] = -p[W];
  168.     t1 = -t1;
  169.     t2 = -t2;
  170.     }
  171.     p[X] = t1 * p1[X] + t2 * p2[X];
  172.     p[Y] = t1 * p1[Y] + t2 * p2[Y];
  173.     p[Z] = t1 * p1[Z] + t2 * p2[Z];
  174. }
  175.  
  176. /*
  177.  * transform a 3d point
  178.  *
  179.  * pt2 = pt1 * [a]
  180.  *
  181.  */
  182. void
  183. Hg4Transform( T, p1, p2)
  184.     Transform3 T;
  185.     Hg4Tensor1 p1, p2;
  186. {
  187.     register Tm3Coord *aptr;
  188.     register Hg4Coord *pptr;
  189.     Hg4Coord x, y, z, w;
  190.     int register cnt;
  191.  
  192.     x = p1[X];
  193.     y = p1[Y];
  194.     z = p1[Z];
  195.     w = p1[W];
  196.     aptr= T[0];
  197.     pptr= p2;
  198.     cnt=4;
  199.     do{
  200.         *pptr++ = aptr[0]*x + aptr[4]*y + aptr[8]*z + aptr[12]*w;
  201.         ++aptr;
  202.     } while(--cnt);
  203. }
  204.  
  205. void
  206. Hg4Print2( L )
  207.     Hg4Tensor2 L;
  208. {
  209.     printf( "[%g %g %g %g\n",  L[X][X], L[X][Y], L[X][Z], L[X][W] );
  210.     printf( " %g %g %g %g\n",  L[Y][X], L[Y][Y], L[Y][Z], L[Y][W] );
  211.     printf( " %g %g %g %g\n",  L[Z][X], L[Z][Y], L[Z][Z], L[Z][W] );
  212.     printf( " %g %g %g %g]\n", L[W][X], L[W][Y], L[W][Z], L[W][W] );
  213. }
  214.  
  215. void
  216. Hg4Copy2( L, K )
  217.     Hg4Tensor2 L, K;
  218. {
  219.     bcopy( (char *)L, (char *)K, sizeof(Hg4Tensor2) );
  220. }
  221.  
  222. int
  223. Hg4Compare2( L, K )
  224.     Hg4Tensor2 L, K;
  225. {
  226.     Hg4Coord t;
  227.     Hg4Tensor2 N;
  228.  
  229.     Hg4ContractPijQjk( K, L, N );
  230.     t = Hg4ContractPii( N );
  231.  
  232.     if( fzero(t) ) return 0;
  233.     if( t < 0.   ) return -1;
  234.     if( t > 0.   ) return  1;
  235. }
  236.  
  237. int
  238. Hg4Undefined2( L )
  239.     Hg4Tensor2 L;
  240. {
  241.     if( !fzero(L[X][Y]) ) return 0;
  242.     if( !fzero(L[X][Z]) ) return 0;
  243.     if( !fzero(L[X][W]) ) return 0;
  244.     if( !fzero(L[Y][Z]) ) return 0;
  245.     if( !fzero(L[Y][W]) ) return 0;
  246.     if( !fzero(L[Z][W]) ) return 0;
  247.     return 1;
  248. }
  249.  
  250. int
  251. Hg4Infinity2( L, dual )
  252.     Hg4Tensor2 L;
  253.     int dual;
  254. {
  255.     /* plane form */
  256.     if( dual ) {
  257.     if( !fzero(L[X][Y]) ) return 0;
  258.     if( !fzero(L[X][Z]) ) return 0;
  259.     if( !fzero(L[Y][Z]) ) return 0;
  260.     return 1;
  261.     }
  262.     else {
  263.     if( !fzero(L[X][W]) ) return 0;
  264.     if( !fzero(L[Y][W]) ) return 0;
  265.     if( !fzero(L[Z][W]) ) return 0;
  266.     return 1;
  267.     }
  268. }
  269.  
  270.  
  271. void
  272. Hg4Transform2( T, p1, p2 )
  273.     Transform3 T;
  274.     Hg4Tensor2 p1, p2;
  275. {
  276.     Transform3 Tt;
  277.  
  278. fprintf(stderr,"\nWARNING: dubious procedure Hg4Transform2 being called.\n\
  279. This procedure may not have been correctly updated for new transform\n\
  280. library.  Ask me about this. --- mbp Mon Aug 19 10:38:19 1991.\n\n");
  281.  
  282. /*
  283.   In fact, this procedure has not been updated at all.  This is the old
  284.   version.   I don't know whether this depends on a notion of col vs row
  285.   vectors, or right vs left mult.  I think it does but I'm not sure how,
  286.   so I'll deal with it later.  -- mbp
  287. */
  288.  
  289.     /* Assume p1 is the plane-form */
  290.     Tm3Transpose( T, Tt );
  291.     Hg4ContractPijQjk( T,  p1, p2 );
  292.     Hg4ContractPijQjk( p2, Tt, p2 ); 
  293. }
  294.  
  295. void
  296. Hg4AntiProductPiQj( L, p1, p2 )
  297.     Hg4Tensor2 L;
  298.     Hg4Tensor1 p1, p2;
  299. {
  300.     L[X][X] = L[Y][Y] = L[Z][Z] = L[W][W] = 0.;
  301.  
  302.     L[X][Y] = p1[X]*p2[Y] - p1[Y]*p2[X];
  303.     L[X][Z] = p1[X]*p2[Z] - p1[Z]*p2[X];
  304.     L[X][W] = p1[X]*p2[W] - p1[W]*p2[X];
  305.     L[Y][Z] = p1[Y]*p2[Z] - p1[Z]*p2[Y];
  306.     L[Y][W] = p1[Y]*p2[W] - p1[W]*p2[Y];
  307.     L[Z][W] = p1[Z]*p2[W] - p1[W]*p2[Z];
  308.  
  309.     L[Y][X] = -L[X][Y];
  310.     L[Z][X] = -L[X][Z];
  311.     L[W][X] = -L[X][W];
  312.     L[Z][Y] = -L[Y][Z];
  313.     L[W][Y] = -L[Y][W];
  314.     L[W][Z] = -L[Z][W];
  315. }
  316.  
  317. Hg4Coord
  318. Hg4ContractPiQi( pl, pt )
  319.     Hg4Tensor1 pl, pt;
  320. {
  321.     Hg4Coord sum;
  322.  
  323.     sum = 0.;
  324.     sum += pl[X] * pt[X];
  325.     sum += pl[Y] * pt[Y];
  326.     sum += pl[Z] * pt[Z];
  327.     sum += pl[W] * pt[W];
  328.  
  329.     return sum;
  330. }
  331.  
  332. void
  333. Hg4AntiContractPijQj( L, p1, p2 )
  334.     Hg4Tensor2 L;
  335.     Hg4Tensor1 p1, p2;
  336. {
  337.     Hg4Coord x, y, z, w;
  338.     Hg4Coord xy, xz, xw, yz, yw, zw;
  339.  
  340.     x = p1[X];
  341.     y = p1[Y];
  342.     z = p1[Z];
  343.     w = p1[W];
  344.  
  345.     xy = L[X][Y];
  346.     xz = L[X][Z];
  347.     xw = L[X][W];
  348.     yz = L[Y][Z];
  349.     yw = L[Y][W];
  350.     zw = L[Z][W];
  351.  
  352.     p2[X] =  xy * y + xz * z + xw * w;
  353.     p2[Y] = -xy * x + yz * z + yw * w;
  354.     p2[Z] = -xz * x - yz * y + zw * w;
  355.     p2[W] = -xw * x - yw * y - zw * z;
  356. }
  357.  
  358. void
  359. Hg4ContractPijQjk( a, b, c )
  360.     Hg4Tensor2 a, b, c;
  361. {
  362.     Hg4Tensor2 d;
  363.     int i, j, k;
  364.  
  365.     /* This can be made more efficient */
  366.     for( i=0; i<4; i++ )
  367.     for( j=0; j<4; j++ ) {
  368.         d[i][j] = 0.;
  369.         for( k=0; k<4; k++ )
  370.         d[i][j] += a[i][k] * b[k][j];
  371.     }
  372.     Hg4Copy2( d, c );
  373. }
  374.  
  375. Hg4Coord
  376. Hg4ContractPii( L )
  377.     Hg4Tensor2 L;
  378. {
  379.     return L[X][X] + L[Y][Y] + L[Z][Z] + L[W][W];
  380. }
  381.  
  382. int
  383. Hg4Intersect2( L, a, b )
  384.     Hg4Tensor2 L;
  385.     Hg4Tensor1 a, b;
  386. {
  387.     Hg4AntiContractPijQj( L, a, b );
  388.     return Hg4Undefined( b );
  389. }
  390.  
  391. int
  392. Hg4Intersect3( a, b, c, p, dual )
  393.     Hg4Tensor1 a, b, c, p;
  394.     int dual;
  395. {
  396.     Hg4Tensor2 L;
  397.  
  398.     Hg4AntiProductPiQj( L, a, b );
  399.     if( dual )
  400.     Hg4Dual( L, L );
  401.  
  402.     Hg4AntiContractPijQj( L, c, p );
  403.  
  404.     return Hg4Undefined( p );
  405. }
  406.  
  407. /*
  408. **    Hg4Intersect4 - predicate which tests for 3d line intersection and
  409. **    if an intersection is found returns the point at which th
  410. **    two lines cross and the plane in which the two lines lie.
  411. **    
  412. **    Note: One of the lines should be in the "plane-form" and the
  413. **    other in the "point-form."
  414. **
  415. **    Assume K is plane-form, L is point-form
  416. */
  417. int
  418. Hg4Intersect4( K, L, pl, pt )
  419.     Hg4Tensor2 K, L;
  420.     Hg4Tensor1 pl;
  421.     Hg4Tensor1 pt;
  422. {
  423.     int flag;
  424.     int i, j;
  425.     Hg4Tensor2 N;
  426.     Hg4Coord t;
  427.  
  428.     Hg4ContractPijQjk( K, L, N );
  429.     Hg4From( pl, 0., 0., 0., 0. );
  430.     Hg4From( pt, 0., 0., 0., 0. );
  431.  
  432.     t = Hg4ContractPii( N );
  433.     if( fzero(t) ) {
  434.     /* Look for a non-zero row */
  435.     flag = 0;
  436.     for( i=0; i<4; i++ ) {
  437.         for( j=0; j<4; j++ ) {
  438.         pt[j] = N[i][j];
  439.         if( !fzero( pt[j] ) ) flag++;
  440.         }
  441.         if( flag ) break;
  442.     }
  443.  
  444.     /* Look for a non-zero col */
  445.     flag = 0;
  446.     for( i=0; i<4; i++ ) {
  447.         for( j=0; j<4; j++ ) {
  448.         pl[j] = N[j][i];
  449.         if( !fzero( pl[j] ) ) flag++;
  450.         }
  451.         if( flag ) break;
  452.     }
  453.     }
  454.  
  455.     return fzero(t);
  456. }
  457.  
  458. void
  459. Hg4Dual( L, K )
  460.     Hg4Tensor2 L, K;
  461. {
  462.     Hg4Coord p, q, r, u, t, s;
  463.  
  464.     /*
  465.     t = xy; xy = -zw; zw = -t;
  466.     t = xz; xz =  yw; yw =  t;
  467.     t = yz; yz = -xw; xw = -t;
  468.     */
  469.  
  470.     if( L != K )
  471.     Hg4Copy2( L, K );
  472.  
  473.     p = K[X][Y]; u = K[Z][W];
  474.     K[Z][W] = -p; K[W][Z] = p;
  475.     K[X][Y] = -u; K[Y][X] = u;
  476.  
  477.     q = K[Z][X]; t = K[W][Y];
  478.     K[Y][W] = -q; K[W][Y] = q;
  479.     K[X][Z] = -t; K[Z][X] = t;
  480.  
  481.     s = K[Y][Z]; r = K[X][W];
  482.     K[X][W] = -s; K[W][X] = s;
  483.     K[Y][Z] = -r; K[Z][Y] = r;
  484. }
  485.  
  486.