ftp.nice.ch/pub/next/graphics/3d/geomview.1.4.1.s.tar.gz#/Geomview/src/lib/geometry/transformn/transformn.c

This is transformn.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, Olaf Holt */

#include "hpointn.h"	/* Defines HPointN and TransformN */
#include "transformn.h"
#include <ooglutil.h>
#include <transform3.h>
#include <math.h>

	/* Construct a transform.  NULL a => identity */
TransformN *
TmNCreate( int idim, int odim, HPtNCoord *a )
{
	TransformN *T = OOGLNewE(TransformN, "new TransformN");
	RefInit((Ref *)T, TMNMAGIC);
	if( idim <= 0 ) idim = 1;
	if( odim <= 0 ) odim = 1;
	T->idim = idim;	 T->odim = odim;
	T->a = OOGLNewNE(HPtNCoord, idim*odim, "new TransformN data");
	if( a == NULL )
		bzero(T->a,idim*odim*sizeof(HPtNCoord));
	else
		memcpy(T->a, a, idim*odim*sizeof(HPtNCoord));
	return(T);
}

	/* Destroy */
void 
TmNDelete( TransformN *T )
{
	if(T && RefDecr((Ref*)T) == 0) {
		if(T->a) OOGLFree(T->a);
		OOGLFree(T);
	}
}

	/* Get and set space */
int 
TmNSpace( const TransformN *T );

TransformN *
TmNSetSpace( TransformN *T, int space );

	/* Invert */
TransformN * 
TmNInvert( const TransformN *T, TransformN *Tinv )
{
	register int i, j, k;
	HPtNCoord x;
	HPtNCoord f;
	short dim = T->idim; 
	TransformN *t = TmNCreate(dim,dim,T->a);

	if(T->odim != dim) {
		OOGLError(1,"Matrix for inversion is not square");
		return(0);
	}
	if(Tinv == NULL)
		Tinv = TmNCreate(dim,dim,NULL);
	else if (Tinv->idim != dim || Tinv->odim != dim) {
		Tinv->a = OOGLRenewNE(HPtNCoord,Tinv->a,dim*dim, "renew TransformN");
		Tinv-> idim = dim; Tinv->odim = dim;
	}

	TmNIdentity(Tinv);

	/* Components of unrolled inner loops: */
#define	SUB(v,k)  v[j*dim+k] -= f*v[i*dim+k]
#define	SWAP(v,k) x = v[i*dim+k], v[i*dim+k] = v[largest*dim+k], v[largest*dim+k] = x


	for (i = 0; i < dim; i++) {
		int largest = i;
		HPtNCoord largesq = t->a[i*dim+i]*t->a[i*dim+i];
		for (j = i+1; j < dim; j++)
			if ((x = t->a[j*dim+i]*t->a[j*dim+i]) > largesq)
				largest = j,  largesq = x;

		/* swap t->a[i][] with t->a[largest][] */
		for(k = 0; k < dim; k++) {
		    SWAP(t->a,k);
		    SWAP(Tinv->a,k);
		}

		for (j = i+1; j < dim; j++) {
			f = t->a[j*dim+i] / t->a[i*dim+i];
			/* subtract f*t->a[i][] from t->a[j][] */
			for(k = 0; k < dim; k++) {
			    SUB(t->a,k);
			    SUB(Tinv->a,k);
			}
		}
	}
	for (i = 0; i < dim; i++) {
		f = t->a[i*dim+i];
		for (k = 0; k < dim; k++) {
			t->a[i*dim+k] /= f;
			Tinv->a[i*dim+k] /= f;
		}
	}
	for (i = dim-1; i >= 0; i--)
		for (j = i-1; j >= 0; j--) {
			f = t->a[j*dim+i];
			for(k = 0; k < dim; k++) {
			    SUB(t->a,k);
			    SUB(Tinv->a,k);
			}
		}
#undef SUB
#undef SWAP
				
	return Tinv;
}

	/* Transpose */
TransformN *
TmNTranspose( const TransformN *from, TransformN *to )
{
	register int i,j;
	short idim = from->idim, odim = from->odim;
	HPtNCoord t;

	if( from != to) {
		if(to == NULL)
			to = TmNCreate(odim,idim,NULL);
		else if (to->idim != odim || to->odim != idim) {
			to->a = OOGLRenewNE(HPtNCoord,to->a,idim*odim, "renew TransformN");
			to-> idim = odim; to->odim = idim;
		}
		for( i=0; i<idim; i++)
			for( j=0; j<odim; j++)
				to->a[j*idim+i] = from->a[i*odim+j];
	} else {

#define SWITCH(a,b)		 t = a; a = b; b = t;						

		if (idim == odim) {
			for( i=0; i<idim; i++)
				for( j= 0; j<i; j++) {
					SWITCH(to->a[i*odim+j], to->a[j*odim+i]);
				}
		} else {
			int remainder, dividend;
			to->idim = odim; to->odim = idim;
			for( i=0; i<idim; i++)
				for( j=0; j<odim; j++) {
					remainder = (i*idim +j)%odim;
					dividend = (i*idim + j -remainder)/odim;
					SWITCH(to->a[i*odim+j],to->a[dividend*odim+remainder]);
				}
		}
#undef SWITCH
	}

	return(to);
}

	/* Multiply transforms */
TransformN *
TmNConcat( const TransformN *A, const TransformN *B, TransformN *result )
{
	register int i,j,k;
	short dim1 = A->idim, dim2 = A->odim, dim3 = B->odim;

#define MAKEPRODUCT(T,A,B)											\
	T->idim = dim1; T->odim = dim3;									\
	for( i=0; i<dim1; i++)											\
		for( j=0; j<dim3; j++) {									\
			T->a[i*dim3+j] = 0;										\
			for( k=0; k<dim2; k++ )									\
				T->a[i*dim3+j] += A->a[i*dim2+k] * B->a[k*dim3+j];	\
		}

	if(B->idim != dim2) {
	  if ( B->idim > dim2) {
		TransformN *Atmp = TmNCreate(dim1,dim2,NULL);
		dim2 = B->idim;
		TmNPadZero((TransformN *)A,dim1,dim2,Atmp);
		if( B == result ) {
			TransformN *T = TmNCreate( dim1, dim3, NULL);

			MAKEPRODUCT(T,Atmp,B);
			TmNCopy(T,result);
		} else {
			if(result == NULL)
				result = TmNCreate(dim1,dim3,NULL);
			else if (result->idim != dim1 || result->odim != dim3) {
				result->a = OOGLRenewNE(HPtNCoord,result->a,dim1*dim3, "renew TransformN");
				result-> idim = dim1; result->odim = dim3;
			}
			MAKEPRODUCT(result,Atmp,B);
		}
	  } else {  /* B->idim < dim2 */
		TransformN *Btmp = TmNCreate(dim2,dim3,NULL);
		TmNPadZero((TransformN *)B,dim2,dim3,Btmp);
		if( A == result ) {
			TransformN *T = TmNCreate( dim1, dim3, NULL);

			MAKEPRODUCT(T,A,Btmp);
			TmNCopy(T,result);
		} else {
			if(result == NULL)
				result = TmNCreate(dim1,dim3,NULL);
			else if (result->idim != dim1 || result->odim != dim3) {
				result->a = OOGLRenewNE(HPtNCoord,result->a,dim1*dim3, "renew TransformN");
				result-> idim = dim1; result->odim = dim3;
			}
			MAKEPRODUCT(result,A,Btmp);
		}
	  }
	} else {

		if( A == result || B == result ) {
			TransformN *T = TmNCreate( dim1, dim3, NULL);

			MAKEPRODUCT(T,A,B);
			TmNCopy(T,result);
		} else {
			if(result == NULL)
				result = TmNCreate(dim1,dim3,NULL);
			else if (result->idim != dim1 || result->odim != dim3) {
				result->a = OOGLRenewNE(HPtNCoord,result->a,dim1*dim3, "renew TransformN");
				result-> idim = dim1; result->odim = dim3;
			}
			MAKEPRODUCT(result,A,B);
		}
	}

#undef MAKEPRODUCT

	return(result);
}

	/* Copy */
TransformN *
TmNCopy( const TransformN *Tsrc, TransformN *Tdst )
{
	if(Tdst == NULL) {
		Tdst = TmNCreate(Tsrc->idim, Tsrc->odim, Tsrc->a);
	} else {
		if(Tdst->idim != Tsrc->idim || Tdst->odim != Tsrc->odim) {
			Tdst->a = OOGLRenewNE(HPtNCoord,Tdst->a,Tsrc->idim*Tsrc->odim, "renew TransformN");
			Tdst->idim = Tsrc->idim; Tdst->odim = Tsrc->odim;
		}
		memcpy(Tdst->a, Tsrc->a, Tsrc->idim*Tsrc->odim*sizeof(HPtNCoord));
	}
	return(Tdst);
}

	/* Set to identity */
TransformN * 
TmNIdentity( TransformN *T )
{


	if(T == NULL) {
		T = TmNCreate(1,1,NULL);
		T->a[0] = 1;
	} else {
  		int i;
	  short idim = T->idim, odim = T->odim;
		bzero(T->a, idim*odim*sizeof(HPtNCoord));
		if (idim == odim || idim < odim) {
			for(i = 0; i < idim; i++)
				T->a[i*odim+i] = 1;
		} else {  /*  idim > odim */
			for(i = 0; i < odim; i++)
				T->a[i*odim+i] = 1;
		}
	}
	return(T);
}


	/* Euclidean translations */
TransformN * 
TmNTranslateOrigin( TransformN *T, const HPointN *p )
{
/* old version, assumes no "w" coordinate at end of p  */
	short dim = p->dim;
	HPointN *pt = HPtNCreate( 1+dim, NULL );
	int i;

	for ( i=0; i< dim; i++)
		pt->v[i] = p->v[i];
	pt->v[dim] = 1;

	TmNTranslate(T,pt);
	return(T);
}

TransformN * 
TmNTranslate( TransformN *T, const HPointN *pt )
{
	register int i;
	HPointN *newpt;
	short dim = pt->dim;

	newpt = HPtNCreate(dim,pt->v);
	HPtNDehomogenize(newpt,newpt);

  if ( T == NULL ) {
	T = TmNCreate(dim,dim,NULL);
	TmNIdentity(T);
	for( i=0; i<dim - 1; i++)
		T->a[(dim-1)*dim+i] = newpt->v[i];
  } else {
	if (T->odim == dim) {
		TmNIdentity(T);
		for( i=0; i<dim - 1; i++)
			T->a[(dim-1)*dim+i] = newpt->v[i];
	} else if (T->odim > dim) {
		short idim = T->idim, odim = T->odim;
		TmNIdentity(T);
		for( i=0; i<dim - 1; i++)
			T->a[(idim-1)*odim+i] = newpt->v[i];
	} else if (T->odim < dim) {
		short idim = T->idim, odim = dim;
		TmNPad(T,idim,odim,T);
		for( i=0; i<dim - 1; i++)
			T->a[(idim-1)*odim+i] = newpt->v[i];
	}
  }

  HPtNDelete(newpt);
  return(T);
}

	/* Pad the transform, with ones down the main diagonal */
TransformN *
TmNPad(TransformN *Tin, short idim, short odim, TransformN *Tout)
{
	if(Tin == NULL) {
		Tout = TmNCreate(idim,odim,NULL);
		TmNIdentity(Tout);
	} else if ( odim <= 0 || idim <= 0 ) {
		;  /* do nothing */
	} else {


	short oldidim = Tin->idim, oldodim = Tin->odim;
	int i,j;
	
	if(Tin != Tout) {
		if (Tout == NULL) {
			Tout = TmNCreate(idim,odim,NULL);
		} else if(Tout->idim != idim || Tout->odim != odim) {
			Tout->a = OOGLRenewNE(HPtNCoord,Tout->a,idim*odim, "renew TransformN");
			Tout->idim = idim; Tout->odim = odim;
		}

	  if( oldidim < idim && oldodim < odim) {
		for ( i=0; i<oldidim; i++) {
			for( j=0; j<oldodim; j++) 
				Tout->a[i*odim+j] = Tin->a[i*oldodim+j];
			for( j=oldodim; j<odim; j++)
				if(i==j)
					Tout->a[i*odim+j] = 1;
				else
					Tout->a[i*odim+j] = 0;
		}
		for ( i=oldidim; i<idim; i++) {
			for( j=0; j<odim; j++)
				if(i==j)
					Tout->a[i*odim+j] = 1;
				else
					Tout->a[i*odim+j] = 0;
		}
	  } else if ( oldidim < idim && oldodim >= odim) {	
		for ( i=0; i<oldidim; i++) {
			for( j=0; j<odim; j++) 
				Tout->a[i*odim+j] = Tin->a[i*oldodim+j];
		}
		for ( i=oldidim; i<idim; i++) {
			for( j=0; j<odim; j++)
				if(i==j)
					Tout->a[i*odim+j] = 1;
				else
					Tout->a[i*odim+j] = 0;
		}
	  } else if ( oldidim >= idim && oldodim < odim) {	
		for ( i=0; i<idim; i++) {
			for( j=0; j<oldodim; j++) 
				Tout->a[i*odim+j] = Tin->a[i*oldodim+j];
			for( j=oldodim; j<odim; j++)
				if(i==j)
					Tout->a[i*odim+j] = 1;
				else
					Tout->a[i*odim+j] = 0;
		}
	  } else if ( oldidim >= idim && oldodim >= odim) {	
		for ( i=0; i<idim; i++)
			for( j=0; j<odim; j++) 
				Tout->a[i*odim+j] = Tin->a[i*oldodim+j];
	  }
	} else {
		TransformN *Tnew = TmNCreate(idim,odim,NULL);
	  if( oldidim < idim && oldodim < odim) {
		for ( i=0; i<oldidim; i++) {
			for( j=0; j<oldodim; j++) 
				Tnew->a[i*odim+j] = Tin->a[i*oldodim+j];
			for( j=oldodim; j<odim; j++)
				if(i==j)
					Tnew->a[i*odim+j] = 1;
		}
		for ( i=oldidim; i<idim; i++) {
			for( j=0; j<odim; j++)
				if(i==j)
					Tnew->a[i*odim+j] = 1;
		}
	  } else if ( oldidim < idim && oldodim >= odim) {	
		for ( i=0; i<oldidim; i++) {
			for( j=0; j<odim; j++) 
				Tnew->a[i*odim+j] = Tin->a[i*oldodim+j];
		}
		for ( i=oldidim; i<idim; i++) {
			for( j=0; j<odim; j++)
				if(i==j)
					Tnew->a[i*odim+j] = 1;
		}
	  } else if ( oldidim >= idim && oldodim < odim) {	
		for ( i=0; i<idim; i++) {
			for( j=0; j<oldodim; j++) 
				Tnew->a[i*odim+j] = Tin->a[i*oldodim+j];
			for( j=oldodim; j<odim; j++)
				if(i==j)
					Tnew->a[i*odim+j] = 1;
		}
	  } else if ( oldidim >= idim && oldodim >= odim) {	
		for ( i=0; i<idim; i++)
			for( j=0; j<odim; j++) 
				Tnew->a[i*odim+j] = Tin->a[i*oldodim+j];
	  }
		TmNCopy(Tnew,Tout);
	}
  }
	return(Tout);
}
	
	/* Pad with zeroes */
TransformN *
TmNPadZero(TransformN *Tin, short idim, short odim, TransformN *Tout)
{
	if(Tin == NULL) {
		Tout = TmNCreate(idim,odim,NULL);
		TmNIdentity(Tout);
	} else if ( odim <= 0 || idim <= 0 ) {
		;  /* do nothing */
	} else {


	short oldidim = Tin->idim, oldodim = Tin->odim;
	int i,j;
	
	if(Tin != Tout) {
		if (Tout == NULL) {
			Tout = TmNCreate(idim,odim,NULL);
		} else if(Tout->idim != idim || Tout->odim != odim) {
			Tout->a = OOGLRenewNE(HPtNCoord,Tout->a,idim*odim, "renew TransformN");
			Tout->idim = idim; Tout->odim = odim;
		}

	  bzero(Tout->a, idim*odim*sizeof(HPtNCoord) );

	  if( oldidim < idim && oldodim < odim) {
		for ( i=0; i<oldidim; i++) {
			for( j=0; j<oldodim; j++) 
				Tout->a[i*odim+j] = Tin->a[i*oldodim+j];
		}
	  } else if ( oldidim < idim && oldodim >= odim) {	
		for ( i=0; i<oldidim; i++) {
			for( j=0; j<odim; j++) 
				Tout->a[i*odim+j] = Tin->a[i*oldodim+j];
		}
	  } else if ( oldidim >= idim && oldodim < odim) {	
		for ( i=0; i<idim; i++) {
			for( j=0; j<oldodim; j++) 
				Tout->a[i*odim+j] = Tin->a[i*oldodim+j];
		}
	  } else if ( oldidim >= idim && oldodim >= odim) {	
		for ( i=0; i<idim; i++)
			for( j=0; j<odim; j++) 
				Tout->a[i*odim+j] = Tin->a[i*oldodim+j];
	  }
	} else {
		TransformN *Tnew = TmNCreate(idim,odim,NULL);
	  if( oldidim < idim && oldodim < odim) {
		for ( i=0; i<oldidim; i++) {
			for( j=0; j<oldodim; j++) 
				Tnew->a[i*odim+j] = Tin->a[i*oldodim+j];
		}
	  } else if ( oldidim < idim && oldodim >= odim) {	
		for ( i=0; i<oldidim; i++) {
			for( j=0; j<odim; j++) 
				Tnew->a[i*odim+j] = Tin->a[i*oldodim+j];
		}
	  } else if ( oldidim >= idim && oldodim < odim) {	
		for ( i=0; i<idim; i++) {
			for( j=0; j<oldodim; j++) 
				Tnew->a[i*odim+j] = Tin->a[i*oldodim+j];
		}
	  } else if ( oldidim >= idim && oldodim >= odim) {	
		for ( i=0; i<idim; i++)
			for( j=0; j<odim; j++) 
				Tnew->a[i*odim+j] = Tin->a[i*oldodim+j];
	  }
		TmNCopy(Tnew,Tout);
	}
  }
	return(Tout);
}
	

	/* Translations by the space of 'pt' */
TransformN * 
TmNSpaceTranslate( TransformN *T, HPointN *pt ) 
{
	/* this has not yet been implemented
	 *  for now, just call the euclidean one */
	return(TmNTranslate( T, pt ));
}

TransformN * 
TmNSpaceTranslateOrigin( TransformN *T, HPointN *pt )
{
	/* this has not yet been implemented
	 *  for now, just call the euclidean one */
	return(TmNTranslate( T, pt ));
}

	/* Scale by the components of 'amount' */
TransformN *
TmNScale( TransformN *T, const HPointN *amount )
{
	int i,j;
	short dim = amount->dim;

	if(T == NULL) {
		T = TmNCreate(dim,dim,NULL);
		TmNIdentity(T);
		for( i=0; i<dim-1; i++)
			T->a[i*dim+i] = amount->v[i];
	} else if( dim != T->idim || dim != T->odim ) {
		short idim = T->idim, odim = T->odim;
		switch( (dim > idim) + 2*(dim>odim) + 4*(idim>odim) ) {
			case 0:  /* odim > idim > dim */
				idim = odim;
				break;
			case 1:  /* odim > dim > idim */
				idim = odim;
				break;
			case 2:  /* impossible */
				OOGLError(1,"incorrect switch");
				break;
			case 3:  /* dim > odim > idim */
				odim = dim; idim = dim;
				break;
			case 4:  /* idim > odim > dim */
				odim = idim;
				break;
			case 5:  /* impossible */
				OOGLError(1,"incorrect switch");
				break;
			case 6:  /* idim > dim > odim */
				odim = idim;
				break;
			case 7:  /* dim > idim > odim */
				odim = dim; idim = dim;
				break;
		}
		TmNPadZero(T,idim,odim,T);
		TmNIdentity(T);
		for( i=0; i<dim-1; i++)
			T->a[i*odim+i] = amount->v[i];
	} else {
		TmNIdentity(T);
		for( i=0; i<dim-1; i++)
			T->a[i*dim+i] = amount->v[i];
	}

	return(T);
}

	/* Construct a geodesic rotation taking vector 'from' to 'toward' */
TransformN *
TmNRotate( TransformN *T, const HPointN *from, const HPointN *toward )
{
	register int i,j,k;
	short dim = from->dim;
	register HPtNCoord len, proj, cosA, sinA;
	HPtNCoord *planebasis1, *planebasis2; /* not homogeneous coords!*/
	HPtNCoord planecoord1, planecoord2;

	planebasis1 = (HPtNCoord *)malloc((dim-1)*sizeof(HPtNCoord));
	planebasis2 = (HPtNCoord *)malloc((dim-1)*sizeof(HPtNCoord));

#define NORMALIZE(v,vout)								\
	len = 0;											\
	for( k=0; k<dim-1; k++)								\
		len += v[k] * v[k];								\
	if ( len == 0 ) {									\
		len = 1;										\
	} else												\
		len = sqrt(len);								\
	for( k=0; k<dim-1; k++)								\
		vout[k] = v[k] / len;

#define DOTPRODUCT(v1,v2,result)						\
	result = 0;											\
	for( i = 0; i<dim -1; i++)                          \
		result += v1[i] * v2[i];
		

  if(dim != toward->dim) {
	if(toward->dim < dim) {
		HPointN *towardtmp = HPtNCreate(dim,NULL);
		HPtNPad((HPointN *)toward,dim,towardtmp);
		if (T == NULL)
			T = TmNCreate(dim,dim,NULL);
		else if (T->idim != dim || T->odim != dim) {
			T->a = OOGLRenewNE(HPtNCoord,T->a,dim*dim, "renew TransformN");
			T->idim = dim;  T->odim = dim;
		}
		bzero(T->a, dim*dim*sizeof(HPtNCoord));
		T->a[dim*dim-1] = 1;
		NORMALIZE(from->v,planebasis1);
		DOTPRODUCT(towardtmp->v,planebasis1,proj);
		for( i=0; i<dim-1; i++)
			planebasis2[i] = towardtmp->v[i] - proj * planebasis1[i];
		NORMALIZE(planebasis2, planebasis2);
		DOTPRODUCT(towardtmp->v,towardtmp->v,len);
		if( (len = sqrt(len)) == 0 ) {
			OOGLError(1,"zero vector unexpected");
			return(NULL);
		}
		cosA = proj / len;
		sinA = sqrt(1-cosA*cosA);
		for( i=0; i<dim-1; i++) {  /* each basis vector */
			planecoord1 = (cosA - 1) * planebasis1[i] + sinA * planebasis2[i];
			planecoord2 = -sinA * planebasis1[i] + (cosA - 1) * planebasis2[i];
			for( j=0; j<dim-1; j++)
				if ( i == j ) T->a[j*dim+i] = 1 + 
						planecoord1*planebasis1[j] + planecoord2*planebasis2[j];
				else T->a[j*dim+i] = planecoord1*planebasis1[j] + 
						planecoord2*planebasis2[j];
		}
	} else {
		HPointN *fromtmp = HPtNCreate(dim,NULL);
		dim = toward->dim;
		HPtNPad((HPointN *)from,dim,fromtmp);
		if (T == NULL)
			T = TmNCreate(dim,dim,NULL);
		else if (T->idim != dim || T->odim != dim) {
			T->a = OOGLRenewNE(HPtNCoord,T->a,dim*dim, "renew TransformN");
			T->idim = dim;  T->odim = dim;
		}
		bzero(T->a, dim*dim*sizeof(HPtNCoord));
		T->a[dim*dim-1] = 1;
		NORMALIZE(fromtmp->v,planebasis1);
		DOTPRODUCT(toward->v,planebasis1,proj);
		for( i=0; i<dim-1; i++)
			planebasis2[i] = toward->v[i] - proj * planebasis1[i];
		NORMALIZE(planebasis2, planebasis2);
		DOTPRODUCT(toward->v,toward->v,len);
		if( (len = sqrt(len)) == 0 ) {
			OOGLError(1,"zero vector unexpected");
			return(NULL);
		}
		cosA = proj / len;
		sinA = sqrt(1-cosA*cosA);
		for( i=0; i<dim-1; i++) {  /* each basis vector */
			planecoord1 = (cosA - 1) * planebasis1[i] + sinA * planebasis2[i];
			planecoord2 = -sinA * planebasis1[i] + (cosA - 1) * planebasis2[i];
			for( j=0; j<dim-1; j++)
				if ( i == j ) T->a[j*dim+i] = 1 + 
						planecoord1*planebasis1[j] + planecoord2*planebasis2[j];
				else T->a[j*dim+i] = planecoord1*planebasis1[j] + 
						planecoord2*planebasis2[j];
		}
	}
  } else {
	if (T == NULL)
		T = TmNCreate(dim,dim,NULL);
	else if (T->idim != dim || T->odim != dim) {
		T->a = OOGLRenewNE(HPtNCoord,T->a,dim*dim, "renew TransformN");
		T->idim = dim;  T->odim = dim;
	}
	bzero(T->a, dim*dim*sizeof(HPtNCoord));
	T->a[dim*dim-1] = 1;

	/* form orthonormal basis for plane */

	NORMALIZE(from->v,planebasis1);
	DOTPRODUCT(toward->v,planebasis1,proj);
	for( i=0; i<dim-1; i++)
		planebasis2[i] = toward->v[i] - proj * planebasis1[i];
	NORMALIZE(planebasis2, planebasis2);
	

	/* now that we have the basis vectors for the plane,
		calculate the rotation matrix */

	DOTPRODUCT(toward->v,toward->v,len);
	if( (len = sqrt(len)) == 0 ) {
		OOGLError(1,"zero vector unexpected");
		return(NULL);
	}
	cosA = proj / len;
	sinA = sqrt(1-cosA*cosA);
	
	for( i=0; i<dim-1; i++) {  /* each basis vector */
		planecoord1 = (cosA - 1) * planebasis1[i] + sinA * planebasis2[i];
		planecoord2 = -sinA * planebasis1[i] + (cosA - 1) * planebasis2[i];
		for( j=0; j<dim-1; j++)
			if ( i == j ) T->a[j*dim+i] = 1 + 
					planecoord1*planebasis1[j] + planecoord2*planebasis2[j];
			else T->a[j*dim+i] = planecoord1*planebasis1[j] + 
					planecoord2*planebasis2[j];
	}
  }
#undef DOTPRODUCT
#undef NORMALIZE

	free(planebasis1);
	free(planebasis2);

	return(T);
}


	/* special operation for computing n-d counterpart for movement
		in a subspace */
/* permute tells which dimensions are projected */
/* delta is the usual 4x4 matrix used in Geomview */
TransformN *
TmNApplyDN( TransformN *mat, int *perm, Transform3 delta)
{
	register HPtNCoord sum;
	short n = mat->idim, d = 4;  /* this is the dimension of the delta matrix */
	short nprime = mat->odim;
	HPtNCoord *tmp;
	int i,j,k;
	int permute[4];

	if( mat->odim != n) {
		if(nprime > n) {
			n = nprime;
			TmNPad(mat,n,n,mat);    /* Note: the input matrix is changed here */
		} else {
			TmNPad(mat,n,n,mat);
		}
	}

	/* Map "-1" in perm[] array to dimension N-1 (homogeneous divisor) */
	for(i = 0; i < 4; i++)
	   permute[i] = (perm[i] >= 0 && perm[i] < n) ? perm[i] : n-1;

	tmp = (HPtNCoord *) malloc(n*d*sizeof(HPtNCoord));

	for(i = 0; i<n; i++) {
		for(j= 0; j<d; j++) {
			tmp[i*d+j] = mat->a[i*n+permute[j] ];
		}
	}

	for(i = 0; i<n; i++) {
		for (j= 0; j<d; j++) {
			sum = 0;
			for(k = 0; k<d; k++) {
				sum += tmp[i*d+k] * delta[k][j];
			}
			mat->a[i*n + permute[j] ] = sum;
		}
	}

	free(tmp);
	return(mat);
}

    /* Return dimensions of a TransformN.  Value is first dimension. */
    /* idim and/or odim may be NULL, in which case they're not returned */
int 
TmNGetSize(const TransformN *T, int *idim, int *odim)
{
	if ( T == NULL)
		return(0);

	if(idim) *idim = T->idim;
	if(odim) *odim = T->odim;
	return(T->idim);
}

void
TmNPrint(FILE *f, const TransformN *T) 
{
    int i, j;
    short idim = T->idim, odim = T->odim;

	if( f == NULL)
		return;
    fprintf(f, "ntransform { %d %d\n", idim, odim);
    for(i = 0; i < idim; i++) {
        for(j = 0; j < odim; j++)
            fprintf(f, "%10.7f ", T->a[i*odim+j]);
        fprintf(f, "\n");
    }
    fprintf(f, "}\n");
}

TransformN *
TmNRead(FILE *f)
{
    HPtNCoord *a;
    int got, n, brack = 0;
    int idim, odim;

    fexpecttoken(f, "ntransform");

    if(fnextc(f,0) == '{')
	brack = fgetc(f);

    if(fgetni(f,1,&idim, 0) <= 0 || fgetni(f, 1, &odim, 0) <= 0 || idim <= 0 || odim <= 0) {
	OOGLSyntax(f, "Expected dimensions of N-D transform");
	return NULL;
    }
    n = idim*odim;
    a = OOGLNewNE(HPtNCoord, n, "new TransformN data");

    got = fgetnf(f, n, a, 0);
    if(got < n) {
	OOGLSyntax(f, "N-D transform: error reading %d'th of %d values.", got, n);
	return NULL;
    }

    if(brack)
	fexpecttoken(f, "}");

    return TmNCreate(idim,odim,a);
}

TransformN *
CtmNScale( HPtNCoord s, TransformN *in, TransformN *out)
{
	short idim = in->idim;
	short odim = in->odim;
	int i,j;

	if(out == in) {
		for( i=0; i<idim-1; i++)
			for( j = 0; j<odim-1; j++)
				out->a[i*odim+j] *= s;
	} else {
		out = TmNCopy(in, out);
		for( i=0; i<idim-1; i++)
			for( j = 0; j<odim-1; j++)
				out->a[i*odim+j] = s*in->a[i*odim+j];
	}
	return(out);
}

These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.