ftp.nice.ch/pub/next/science/mathematics/hippoplotamus.2.0.s.tar.gz#/hippo2.0/h3Dmat.c

This is h3Dmat.c in view mode; [Download] [Up]

#include <string.h>
#include <math.h>

#include "h3D.h"

/*
 * Internal C functions 
 */
void h3D_translate(display disp, float dx, float dy, float dz);
void h3D_scale(display disp, float sx, float sy, float sz);
void h3D_xRotationCS(display disp);
void h3D_xRotationResetCS(display disp);
void h3D_yRotationCS(display disp);
void h3D_zRotationCS(display disp, float cosAngle, float sinAngle);
void h3D_perspective(display disp);
void h3D_setUnitMatrix(display disp);

void 
h3D_render3D(display disp, int n, float *x, float *y, float *z, float **results)
{
    threeD_t	       *threeD;    	 
    float               h, *u = results[0], *v = results[1], *w = results[2];
    register float     *mat;

    threeD = disp->threeDWorkArea;	/* set local pointer */

    mat = threeD->matrix;
    while (n--) {
	h = mat[12] * *x + mat[13] * *y + mat[14] * *z + mat[15];
	*u++ = (mat[0] * *x + mat[1] * *y + mat[2] * *z + mat[3]) / h;
	*v++ = (mat[4] * *x++ + mat[5] * *y++ + mat[6] * *z++ +
		mat[7]) / h;
	*w++ = h;
    }
    return;
}

/* same as above, but not interested in depth information */
void 
h3D_render2D(display disp, int n, float *x, float *y, float *z, float **results)
{
    threeD_t	       *threeD;    	 
    float		h, *u = results[0], *v = results[1];
    register float     *mat;

    threeD = disp->threeDWorkArea;	/* set local pointer */

    mat = threeD->matrix;
    while (n--) {
	h = mat[12] * *x + mat[13] * *y + mat[14] * *z + mat[15];
	*u++ = (mat[0] * *x + mat[1] * *y + mat[2] * *z + mat[3]) / h;
	*v++ = (mat[4] * *x++ + mat[5] * *y++ + mat[6] * *z++ +
		mat[7]) / h;
    }
    return;
}

void 
h3D_renderXYShow(display disp, float *coords)
{
    threeD_t	       *threeD;    	 
    float               h;
    float               prevX, prevY, newX, newY;
    register float     *mat, *res, *x, *y, *z;
    int 		nCoords;
        
    threeD = disp->threeDWorkArea;	/* set local pointer */

    nCoords = threeD->nScatterPoints;
    if (nCoords <= 0)
	return;
    x = &coords[0];
    y = &coords[1];
    z = &coords[2];

    mat = threeD->matrix;
    res = threeD->scatterResults;

 /* save first XY pair */
    h = mat[12] * *x + mat[13] * *y + mat[14] * *z + mat[15];
    prevX = (mat[0] * *x + mat[1] * *y + mat[2] * *z + mat[3]) / h;
    prevY = (mat[4] * *x + mat[5] * *y + mat[6] * *z + mat[7]) / h;
    *res++ = prevX;
    *res++ = prevY;
    x +=3;
    y +=3;
    z +=3;
    nCoords--;

 /* now save the increments from previous coord */
    while (nCoords--) {
	h = mat[12] * *x + mat[13] * *y + mat[14] * *z + mat[15];
	newX = (mat[0] * *x + mat[1] * *y + mat[2] * *z + mat[3]) / h;
	newY = (mat[4] * *x + mat[5] * *y + mat[6] * *z + mat[7]) / h;
	*res++ = newX - prevX;
	*res++ = newY - prevY;
	prevX = newX;
	prevY = newY;
    	x +=3;
    	y +=3;
    	z +=3;
    }
    *res++ = 0.0;
    *res++ = 0.0;
    return;
}

/* Same as above but just save coords, not differences. Used for PS output */
void 
h3D_renderXYPS(display disp, float *coords)
{
    threeD_t	       *threeD;    	 
    float               h;
    register float     *mat, *res, *x, *y, *z;
    int 		nCoords;
        
    threeD = disp->threeDWorkArea;	/* set local pointer */

    nCoords = threeD->nScatterPoints;
    if (nCoords <= 0)
	return;
    x = &coords[0];
    y = &coords[1];
    z = &coords[2];

    mat = threeD->matrix;
    res = threeD->scatterResults;

    while (nCoords--) {
	h = mat[12] * *x + mat[13] * *y + mat[14] * *z + mat[15];
	*res++ = (mat[0] * *x + mat[1] * *y + mat[2] * *z + mat[3]) / h;
	*res++ = (mat[4] * *x + mat[5] * *y + mat[6] * *z + mat[7]) / h;
    	x +=3;
    	y +=3;
    	z +=3;
    }
    return;
}

void 
h3D_translate(display disp, float dx, float dy, float dz)
{
    threeD_t	       *threeD;    	 
    register float     *mat, *thisrow, *lastrow;
    register short int  n = 4;

    threeD = disp->threeDWorkArea;	/* set local pointer */

    mat = threeD->matrix;
    thisrow = mat;
    lastrow = mat + 12;
    while (n--)
	*thisrow++ += dx * *lastrow++;
    n = 4;
    lastrow -= 4;
    while (n--)
	*thisrow++ += dy * *lastrow++;
    n = 4;
    lastrow -= 4;
    while (n--)
	*thisrow++ += dz * *lastrow++;
    return;
}

void 
h3D_scale(display disp, float sx, float sy, float sz)
{
    threeD_t	       *threeD;    	 
    register float     *mat, *thisrow;
    register short int  n = 4;

    threeD = disp->threeDWorkArea;	/* set local pointer */

    mat = threeD->matrix;
    thisrow = mat;
    while (n--)
	*thisrow++ *= sx;
    n = 4;
    while (n--)
	*thisrow++ *= sy;
    n = 4;
    while (n--)
	*thisrow++ *= sz;
    return;
}

void 
h3D_xRotationCS(display disp)
{
    threeD_t	       *threeD;    	 
    register float      temp, *mat, *thisrow, *otherrow;
    register short int  n = 4;

    threeD = disp->threeDWorkArea;	/* set local pointer */

    mat = threeD->matrix;
    thisrow = mat + 4;
    otherrow = mat + 8;
    while (n--) {
	temp = threeD->cosPhi * *thisrow - threeD->sinPhi * *otherrow;
	*otherrow++ = threeD->sinPhi * *thisrow + threeD->cosPhi * *otherrow;
	*thisrow++ = temp;
    }
    return;
}

void 
h3D_xRotationResetCS(display disp)
{
    threeD_t	       *threeD;    	 
    register float      temp, *mat, *thisrow, *otherrow;
    register short int  n = 4;

    threeD = disp->threeDWorkArea;	/* set local pointer */

    mat = threeD->matrix;
    thisrow = mat + 4;
    otherrow = mat + 8;
    while (n--) {
	temp = *otherrow;
	*otherrow++ = -*thisrow;
	*thisrow++ = temp;
    }
    return;
}

void 
h3D_yRotationCS(display disp)
{
    threeD_t	       *threeD;    	 
    register float      temp, *mat, *thisrow, *otherrow;
    register short int  n = 4;

    threeD = disp->threeDWorkArea;	/* set local pointer */

    mat = threeD->matrix;
    thisrow = mat;
    otherrow = mat + 8;
    while (n--) {
	temp = threeD->cosTheta * *thisrow + threeD->sinTheta * *otherrow;
	*otherrow++ = -threeD->sinTheta * *thisrow + threeD->cosTheta * *otherrow;
	*thisrow++ = temp;
    }
    return;
}

void 
h3D_zRotationCS(display disp, float cosAngle, float sinAngle)
{
    threeD_t	       *threeD;    	 
    register float      temp, *mat, *thisrow, *otherrow;
    register short int  n = 4;

    threeD = disp->threeDWorkArea;	/* set local pointer */

    mat = threeD->matrix;
    thisrow = mat;
    otherrow = mat + 4;
    while (n--) {
	temp = cosAngle * *thisrow - sinAngle * *otherrow;
	*otherrow++ = sinAngle * *thisrow + cosAngle * *otherrow;
	*thisrow++ = temp;
    }
    return;
}

void 
h3D_perspective(display disp)
{
    threeD_t	       *threeD;    	 
    register float 	inv_d, *mat, *thisrow, *otherrow;
    register short int  n = 4;

    threeD = disp->threeDWorkArea;	/* set local pointer */

    mat = threeD->matrix;
    thisrow = mat + 12;
    otherrow = mat + 8;
    inv_d = 1.0/disp->dist;
    while (n--) {
	*thisrow++ -= inv_d * *otherrow;
	*otherrow++ = 0.;
    }
    return;
}

/*
 * Reset matrix to unity
 */
void 
h3D_setUnitMatrix(display disp)
{
    threeD_t	       *threeD;    	 
    register float     *mat;
    register short int  n = 16;

    threeD = disp->threeDWorkArea;	/* set local pointer */

    mat = threeD->matrix;
    while (n--) *mat++ = 0.0;

    mat = threeD->matrix;
    mat[0] = mat[5] = mat[10] = mat[15] = 1.0;
    return;
}

/*
 * The starting point for rotations and dist changes
 */
void 
h3D_reScaleMatrix(display disp)
{
    threeD_t	       *threeD;    	 

    threeD = disp->threeDWorkArea;	/* set local pointer */

    h3D_setUnitMatrix(disp);
    h3D_translate(disp, -(threeD->limits[0] + threeD->limits[3]) / 2., -(threeD->limits[1] + threeD->limits[4]) / 2.,
-(threeD->limits[2] + threeD->limits[5]) / 2.);
    h3D_scale(disp, 1./ (threeD->limits[3] - threeD->limits[0]), 1./ (threeD->limits[4] - threeD->limits[1]),
1./ (threeD->limits[5] - threeD->limits[2]));
    h3D_xRotationResetCS(disp);
    memmove(threeD->thetaSave, threeD->matrix, 16 * sizeof(float));
    threeD->thetaIsSaved = 1;
    threeD->phiIsSaved = threeD->distIsSaved = 0;
    return;
}

/*
 * Starting point followed by theta, phi, and dist changes.
 * Can skip the starting point if theta is cached.
 */
void 
h3D_newTheta(display disp)
{
    threeD_t	       *threeD;    	 

    threeD = disp->threeDWorkArea;	/* set local pointer */

    if (threeD->thetaIsSaved)
	memmove(threeD->matrix, threeD->thetaSave, 16 * sizeof(float));
    else
	h3D_reScaleMatrix(disp);
    threeD->phiIsSaved = threeD->distIsSaved = 0;
    threeD->cosTheta = cos(disp->theta);
    threeD->sinTheta = sin(disp->theta);
    h3D_yRotationCS(disp); 	/* theta change */
    h3D_xRotationCS(disp);	/* phi change   */
    h3D_perspective(disp);	/* dist change  */

 /* Have to reset the 'near face' every time theta changes */
    h3D_setNearFace(disp);

 /* Show that display re-calculation will be needed */
    h3D_setDirty(disp);

    return;
}

/*
 * Starting point followed by theta, phi, and dist changes.
 * Can skip starting point + theta change, if phi is cached. 
 */
void 
h3D_newPhi(display disp)
{
    threeD_t	       *threeD;    	 

    threeD = disp->threeDWorkArea;	/* set local pointer */

    if (threeD->phiIsSaved)
	memmove(threeD->matrix, threeD->phiSave, 16 * sizeof(float));
    else {
	if (threeD->thetaIsSaved)
	    memmove(threeD->matrix, threeD->thetaSave, 16 * sizeof(float));
	else
	    h3D_reScaleMatrix(disp);
	h3D_yRotationCS(disp);	/* theta change */
	memmove(threeD->phiSave, threeD->matrix, 16 * sizeof(float));
	threeD->phiIsSaved = 1;
    }
    threeD->distIsSaved = 0;
    threeD->cosPhi = cos(disp->phi);
    threeD->sinPhi = sin(disp->phi);
    h3D_xRotationCS(disp);	/* phi change   */
    h3D_perspective(disp);	/* dist change  */

 /* Show that display re-calculation will be needed */
    h3D_setDirty(disp);

    return;
}

/*
 * Starting point followed by theta, phi, and dist changes.
 * Can skip the starting point + theta + phi changes, if dist is cached
 */
void 
h3D_newDist(display disp)
{
    threeD_t	       *threeD;    	 

    threeD = disp->threeDWorkArea;	/* set local pointer */

    if (threeD->distIsSaved)
	memmove(threeD->matrix, threeD->distSave, 16 * sizeof(float));
    else {
	if (threeD->thetaIsSaved)
	    memmove(threeD->matrix, threeD->thetaSave, 16 * sizeof(float));
	else
	    h3D_reScaleMatrix(disp);
	h3D_yRotationCS(disp);	/* theta change */
	h3D_xRotationCS(disp);	/* phi change   */
	memmove(threeD->distSave, threeD->matrix, 16 * sizeof(float));
	threeD->distIsSaved = 1;
    }
    h3D_perspective(disp);	/* dist change  */

 /* Show that display re-calculation will be needed */
    h3D_setDirty(disp);

    return;
}

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