This is hg4.c in view mode; [Download] [Up]
/* Copyright (c) 1992 The Geometry Center; University of Minnesota
1300 South Second Street; Minneapolis, MN 55454, USA;
This file is part of geomview/OOGL. geomview/OOGL is free software;
you can redistribute it and/or modify it only under the terms given in
the file COPYING, which you should have received along with this file.
This and other related software may be obtained via anonymous ftp from
geom.umn.edu; email: software@geom.umn.edu. */
/* Authors: Charlie Gunn, Pat Hanrahan, Stuart Levy, Tamara Munzner, Mark Phillips */
#
/*
** hg4.h - procedural interface to homogeneous geometry
**
** pat hanrahan
**
*/
# include <math.h>
# include <stdio.h>
# include "tolerance.h"
# include "transform3.h"
# include "hg4.h"
char *
Hg4Create()
{
return (char *) malloc( sizeof(Hg4Tensor1) );
}
void
Hg4Delete( p )
Hg4Tensor1 p;
{
free( (char *) p );
}
void
Hg4Print( p )
Hg4Tensor1 p;
{
if( p )
printf( "%g %g %g %g\n", p[X], p[Y], p[Z], p[W] );
}
void
Hg4From( p, x, y, z, w )
Hg4Tensor1 p;
Hg4Coord x, y, z, w;
{
p[X] = x;
p[Y] = y;
p[Z] = z;
p[W] = w;
}
void
Hg4Copy( a, b )
Hg4Tensor1 a, b;
{
bcopy( (char *)a, (char *)b, sizeof(Hg4Tensor1) );
}
void
Hg4Add( p1, p2, p3)
Hg4Tensor1 p1, p2, p3;
{
register int i;
for (i=0; i<4; ++i)
p3[i] = p1[i] + p2[i];
}
int
Hg4Compare( p1, p2 )
Hg4Tensor1 p1, p2;
{
Hg4Coord test;
test = p1[X]*p2[Y] - p1[Y]*p2[X];
if( fneg(test) ) return -1;
if( fpos(test) ) return 1;
test = p1[X]*p2[Z] - p1[Z]*p2[X];
if( fneg(test) ) return -1;
if( fpos(test) ) return 1;
test = p1[Y]*p2[Z] - p1[Z]*p2[Y];
if( fneg(test) ) return -1;
if( fpos(test) ) return 1;
test = p1[X]*p2[W] - p1[W]*p2[X];
if( fneg(test) ) return -1;
if( fpos(test) ) return 1;
test = p1[Y]*p2[W] - p1[W]*p2[Y];
if( fneg(test) ) return -1;
if( fpos(test) ) return 1;
test = p1[Z]*p2[W] - p1[W]*p2[Z];
if( fneg(test) ) return -1;
if( fpos(test) ) return 1;
return 0;
}
int
Hg4Coincident( p1, p2 )
Hg4Tensor1 p1;
Hg4Tensor1 p2;
{
return Hg4Compare( p1, p2 ) == 0;
}
int
Hg4Undefined( a )
Hg4Tensor1 a;
{
if( !fzero(a[X]) ) return 0;
if( !fzero(a[Y]) ) return 0;
if( !fzero(a[Z]) ) return 0;
if( !fzero(a[W]) ) return 0;
return 1;
}
int
Hg4Infinity( p, dual )
Hg4Tensor1 p;
int dual;
{
/* Assume not undefined */
if( dual ) { /* plane */
if( !fzero(p[X]) ) return 0;
if( !fzero(p[Y]) ) return 0;
if( !fzero(p[Z]) ) return 0;
return 1;
}
else { /* point */
if( !fzero(p[W]) ) return 0;
return 1;
}
}
void
Hg4Normalize( p, q )
Hg4Tensor1 p, q;
{
Hg4Copy( p, q );
if( q[W] != 1. && q[W] != 0. ) {
q[X] /= q[W];
q[Y] /= q[W];
q[Z] /= q[W];
q[W] = 1.;
}
}
void
Hg4Pencil( t1, p1, t2, p2, p )
Hg4Coord t1, t2;
Hg4Tensor1 p1, p2, p;
{
p[W] = t1 * p1[W] + t2 * p2[W];
/* Keep W positive */
if( p[W] < 0. ) {
p[W] = -p[W];
t1 = -t1;
t2 = -t2;
}
p[X] = t1 * p1[X] + t2 * p2[X];
p[Y] = t1 * p1[Y] + t2 * p2[Y];
p[Z] = t1 * p1[Z] + t2 * p2[Z];
}
/*
* transform a 3d point
*
* pt2 = pt1 * [a]
*
*/
void
Hg4Transform( T, p1, p2)
Transform3 T;
Hg4Tensor1 p1, p2;
{
register Tm3Coord *aptr;
register Hg4Coord *pptr;
Hg4Coord x, y, z, w;
int register cnt;
x = p1[X];
y = p1[Y];
z = p1[Z];
w = p1[W];
aptr= T[0];
pptr= p2;
cnt=4;
do{
*pptr++ = aptr[0]*x + aptr[4]*y + aptr[8]*z + aptr[12]*w;
++aptr;
} while(--cnt);
}
void
Hg4Print2( L )
Hg4Tensor2 L;
{
printf( "[%g %g %g %g\n", L[X][X], L[X][Y], L[X][Z], L[X][W] );
printf( " %g %g %g %g\n", L[Y][X], L[Y][Y], L[Y][Z], L[Y][W] );
printf( " %g %g %g %g\n", L[Z][X], L[Z][Y], L[Z][Z], L[Z][W] );
printf( " %g %g %g %g]\n", L[W][X], L[W][Y], L[W][Z], L[W][W] );
}
void
Hg4Copy2( L, K )
Hg4Tensor2 L, K;
{
bcopy( (char *)L, (char *)K, sizeof(Hg4Tensor2) );
}
int
Hg4Compare2( L, K )
Hg4Tensor2 L, K;
{
Hg4Coord t;
Hg4Tensor2 N;
Hg4ContractPijQjk( K, L, N );
t = Hg4ContractPii( N );
if( fzero(t) ) return 0;
if( t < 0. ) return -1;
if( t > 0. ) return 1;
}
int
Hg4Undefined2( L )
Hg4Tensor2 L;
{
if( !fzero(L[X][Y]) ) return 0;
if( !fzero(L[X][Z]) ) return 0;
if( !fzero(L[X][W]) ) return 0;
if( !fzero(L[Y][Z]) ) return 0;
if( !fzero(L[Y][W]) ) return 0;
if( !fzero(L[Z][W]) ) return 0;
return 1;
}
int
Hg4Infinity2( L, dual )
Hg4Tensor2 L;
int dual;
{
/* plane form */
if( dual ) {
if( !fzero(L[X][Y]) ) return 0;
if( !fzero(L[X][Z]) ) return 0;
if( !fzero(L[Y][Z]) ) return 0;
return 1;
}
else {
if( !fzero(L[X][W]) ) return 0;
if( !fzero(L[Y][W]) ) return 0;
if( !fzero(L[Z][W]) ) return 0;
return 1;
}
}
void
Hg4Transform2( T, p1, p2 )
Transform3 T;
Hg4Tensor2 p1, p2;
{
Transform3 Tt;
fprintf(stderr,"\nWARNING: dubious procedure Hg4Transform2 being called.\n\
This procedure may not have been correctly updated for new transform\n\
library. Ask me about this. --- mbp Mon Aug 19 10:38:19 1991.\n\n");
/*
In fact, this procedure has not been updated at all. This is the old
version. I don't know whether this depends on a notion of col vs row
vectors, or right vs left mult. I think it does but I'm not sure how,
so I'll deal with it later. -- mbp
*/
/* Assume p1 is the plane-form */
Tm3Transpose( T, Tt );
Hg4ContractPijQjk( T, p1, p2 );
Hg4ContractPijQjk( p2, Tt, p2 );
}
void
Hg4AntiProductPiQj( L, p1, p2 )
Hg4Tensor2 L;
Hg4Tensor1 p1, p2;
{
L[X][X] = L[Y][Y] = L[Z][Z] = L[W][W] = 0.;
L[X][Y] = p1[X]*p2[Y] - p1[Y]*p2[X];
L[X][Z] = p1[X]*p2[Z] - p1[Z]*p2[X];
L[X][W] = p1[X]*p2[W] - p1[W]*p2[X];
L[Y][Z] = p1[Y]*p2[Z] - p1[Z]*p2[Y];
L[Y][W] = p1[Y]*p2[W] - p1[W]*p2[Y];
L[Z][W] = p1[Z]*p2[W] - p1[W]*p2[Z];
L[Y][X] = -L[X][Y];
L[Z][X] = -L[X][Z];
L[W][X] = -L[X][W];
L[Z][Y] = -L[Y][Z];
L[W][Y] = -L[Y][W];
L[W][Z] = -L[Z][W];
}
Hg4Coord
Hg4ContractPiQi( pl, pt )
Hg4Tensor1 pl, pt;
{
Hg4Coord sum;
sum = 0.;
sum += pl[X] * pt[X];
sum += pl[Y] * pt[Y];
sum += pl[Z] * pt[Z];
sum += pl[W] * pt[W];
return sum;
}
void
Hg4AntiContractPijQj( L, p1, p2 )
Hg4Tensor2 L;
Hg4Tensor1 p1, p2;
{
Hg4Coord x, y, z, w;
Hg4Coord xy, xz, xw, yz, yw, zw;
x = p1[X];
y = p1[Y];
z = p1[Z];
w = p1[W];
xy = L[X][Y];
xz = L[X][Z];
xw = L[X][W];
yz = L[Y][Z];
yw = L[Y][W];
zw = L[Z][W];
p2[X] = xy * y + xz * z + xw * w;
p2[Y] = -xy * x + yz * z + yw * w;
p2[Z] = -xz * x - yz * y + zw * w;
p2[W] = -xw * x - yw * y - zw * z;
}
void
Hg4ContractPijQjk( a, b, c )
Hg4Tensor2 a, b, c;
{
Hg4Tensor2 d;
int i, j, k;
/* This can be made more efficient */
for( i=0; i<4; i++ )
for( j=0; j<4; j++ ) {
d[i][j] = 0.;
for( k=0; k<4; k++ )
d[i][j] += a[i][k] * b[k][j];
}
Hg4Copy2( d, c );
}
Hg4Coord
Hg4ContractPii( L )
Hg4Tensor2 L;
{
return L[X][X] + L[Y][Y] + L[Z][Z] + L[W][W];
}
int
Hg4Intersect2( L, a, b )
Hg4Tensor2 L;
Hg4Tensor1 a, b;
{
Hg4AntiContractPijQj( L, a, b );
return Hg4Undefined( b );
}
int
Hg4Intersect3( a, b, c, p, dual )
Hg4Tensor1 a, b, c, p;
int dual;
{
Hg4Tensor2 L;
Hg4AntiProductPiQj( L, a, b );
if( dual )
Hg4Dual( L, L );
Hg4AntiContractPijQj( L, c, p );
return Hg4Undefined( p );
}
/*
** Hg4Intersect4 - predicate which tests for 3d line intersection and
** if an intersection is found returns the point at which th
** two lines cross and the plane in which the two lines lie.
**
** Note: One of the lines should be in the "plane-form" and the
** other in the "point-form."
**
** Assume K is plane-form, L is point-form
*/
int
Hg4Intersect4( K, L, pl, pt )
Hg4Tensor2 K, L;
Hg4Tensor1 pl;
Hg4Tensor1 pt;
{
int flag;
int i, j;
Hg4Tensor2 N;
Hg4Coord t;
Hg4ContractPijQjk( K, L, N );
Hg4From( pl, 0., 0., 0., 0. );
Hg4From( pt, 0., 0., 0., 0. );
t = Hg4ContractPii( N );
if( fzero(t) ) {
/* Look for a non-zero row */
flag = 0;
for( i=0; i<4; i++ ) {
for( j=0; j<4; j++ ) {
pt[j] = N[i][j];
if( !fzero( pt[j] ) ) flag++;
}
if( flag ) break;
}
/* Look for a non-zero col */
flag = 0;
for( i=0; i<4; i++ ) {
for( j=0; j<4; j++ ) {
pl[j] = N[j][i];
if( !fzero( pl[j] ) ) flag++;
}
if( flag ) break;
}
}
return fzero(t);
}
void
Hg4Dual( L, K )
Hg4Tensor2 L, K;
{
Hg4Coord p, q, r, u, t, s;
/*
t = xy; xy = -zw; zw = -t;
t = xz; xz = yw; yw = t;
t = yz; yz = -xw; xw = -t;
*/
if( L != K )
Hg4Copy2( L, K );
p = K[X][Y]; u = K[Z][W];
K[Z][W] = -p; K[W][Z] = p;
K[X][Y] = -u; K[Y][X] = u;
q = K[Z][X]; t = K[W][Y];
K[Y][W] = -q; K[W][Y] = q;
K[X][Z] = -t; K[Z][X] = t;
s = K[Y][Z]; r = K[X][W];
K[X][W] = -s; K[W][X] = s;
K[Y][Z] = -r; K[Z][Y] = r;
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.