proj3d.c

BCView - Bayes Classifier Visualization Download xbcview Linux executable (218 kb) wbcview.exe W

/*----------------------------------------------------------------------
  File    : proj3d.c
  Contents: projection of 3D points to a 2D plane
  Author  : Christian Borgelt
  History : 19.10.2000 file created from file drawing.c
            10.12.2000 function p3d_visib added
            15.06.2001 complete redesign
            09.07.2001 bug in function p3d_view removed
----------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "proj3d.h"

/*----------------------------------------------------------------------
  Preprocessor Definitions
----------------------------------------------------------------------*/
#define	M_PI   3.14159265358979323846

/*----------------------------------------------------------------------
  Auxiliary Functions
----------------------------------------------------------------------*/

void _rotmat (double mat[3][3], double vec[3], double phi)
{                               /* --- compute rotation matrix */
  double sin_p, cos_p;          /* sin(phi) and cos(phi) */
  double _1_cos_p;              /* 1 - cos(phi) */
  double t;                     /* temporary buffer */

  sin_p = sin(phi); cos_p = cos(phi); _1_cos_p = 1 -cos_p;
  mat[0][0] = _1_cos_p *vec[0]*vec[0] +cos_p;
  mat[1][1] = _1_cos_p *vec[1]*vec[1] +cos_p;
  mat[2][2] = _1_cos_p *vec[2]*vec[2] +cos_p;
  mat[1][0] = mat[0][1] = _1_cos_p *vec[0]*vec[1];
  mat[2][0] = mat[0][2] = _1_cos_p *vec[0]*vec[2];
  mat[2][1] = mat[1][2] = _1_cos_p *vec[1]*vec[2];
  mat[1][0] += t = sin_p *vec[2]; mat[0][1] -= t;
  mat[2][0] -= t = sin_p *vec[1]; mat[0][2] += t;
  mat[2][1] += t = sin_p *vec[0]; mat[1][2] -= t;
}  /* _rotmat() */

/*--------------------------------------------------------------------*/

void _rotate (double mat[3][3], double vec[3])
{                               /* --- rotate a vector */
  double x, y;                  /* buffers for coordinates */

  x      = mat[0][0]*vec[0] + mat[1][0]*vec[1] + mat[2][0]*vec[2];
  y      = mat[0][1]*vec[0] + mat[1][1]*vec[1] + mat[2][1]*vec[2];
  vec[2] = mat[0][2]*vec[0] + mat[1][2]*vec[1] + mat[2][2]*vec[2];
  vec[1] = y; vec[0] = x;       /* matrix multiplication from right */
}  /* _rotate() */

/*--------------------------------------------------------------------*/

void _rotcrd (double mat[3][3], double a[3], double b[3], double c[3])
{                               /* --- rotate a coordinate cross */
  double x, y, z;               /* buffers for coordinates */
  double t;                     /* recipocal length of vector */

  x = mat[0][0]*b[0] + mat[1][0]*b[1] + mat[2][0]*b[2];
  y = mat[0][1]*b[0] + mat[1][1]*b[1] + mat[2][1]*b[2];
  z = mat[0][2]*b[0] + mat[1][2]*b[1] + mat[2][2]*b[2];
  t = 1/sqrt(x*x +y*y +z*z);    /* matrix multiplication from right */
  b[0] = x*t; b[1] = y*t; b[2] = z*t;
  c[0] = a[1]*b[2] - a[2]*b[1]; /* renormalize the rotated vector */
  c[1] = a[2]*b[0] - a[0]*b[2]; /* and compute the third as the */
  c[2] = a[0]*b[1] - a[1]*b[0]; /* cross product of the other two */
}  /* _rotcrd() */            

/*----------------------------------------------------------------------
  Main Functions
----------------------------------------------------------------------*/

PROJ3D* p3d_create (void)
{                               /* --- create a 3D to 2D projection */
  PROJ3D *p3d;                  /* created projection */

  p3d = (PROJ3D*)malloc(sizeof(PROJ3D));
  if (!p3d) return NULL;        /* allocate memory for the projection */
  p3d->curr.eye[0] = 0; p3d->curr.eye[1] = -1; p3d->curr.eye[2] = 0;
  p3d->curr.cx[0]  = 1; p3d->curr.cx[1]  =  0; p3d->curr.cx[2]  = 0;
  p3d->curr.cy[0]  = 0; p3d->curr.cy[1]  =  1; p3d->curr.cy[2]  = 0;
  p3d->curr.cz[0]  = 0; p3d->curr.cz[1]  =  0; p3d->curr.cz[2]  = 1;
  p3d->curr.dnear  = 1e-6;      /* set default viewing parameters */
  p3d->curr.dfar   = 1e+6;      /* (look towards the origin) */
  p3d->curr.dpp    = 1/tan(0.125*M_PI);
  return p3d;                   /* return the created projection */
}  /* p3d_create() */
  
/*--------------------------------------------------------------------*/

void p3d_eye (PROJ3D *p3d, double x, double y, double z)
{                               /* --- set eye (point to look from) */
  p3d->curr.eye[0] = x; p3d->curr.eye[1] = y; p3d->curr.eye[2] = z;
}  /* p3d_eye() */

/*--------------------------------------------------------------------*/

void p3d_view (PROJ3D *p3d, double h, double p, double r)
{                               /* --- set direction of view */
  double sin_h, cos_h;          /* sine/cosine of heading angle */
  double sin_p, cos_p;          /* sine/cosine of pitch   angle */
  double sin_r, cos_r;          /* sine/cosine of roll    angle */

  sin_h = sin(h); cos_h = cos(h);
  sin_p = sin(p); cos_p = cos(p);  /* compute sines and cosines */
  sin_r = sin(r); cos_r = cos(r);  /* and rotate the unit vectors */
  p3d->curr.cx[0] =  cos_h*cos_r + sin_h*sin_p*sin_r;
  p3d->curr.cx[1] =  sin_h*cos_r - cos_h*sin_p*sin_r;
  p3d->curr.cx[2] =  cos_p*sin_r;  /* x-axis for viewing */
  p3d->curr.cy[0] = -sin_h*cos_p;
  p3d->curr.cy[1] =  cos_h*cos_p;
  p3d->curr.cy[2] =  sin_p;        /* y-axis for viewing */
  p3d->curr.cz[0] = -cos_h*sin_r + sin_h*sin_p*cos_r;
  p3d->curr.cz[1] = -sin_h*sin_r - cos_h*sin_p*cos_r;
  p3d->curr.cz[2] =  cos_p*cos_r;  /* z-axis for viewing */
}  /* p3d_view() */

/*--------------------------------------------------------------------*/

void p3d_frust (PROJ3D *p3d, double fov, double dnear, double dfar)
{                               /* --- set viewing frustrum */
  p3d->curr.dpp = 1/tan(0.5*fov);  /* compute distance to proj. plane */
  p3d->curr.dnear = dnear;      /* note the distances */
  p3d->curr.dfar  = dfar;       /* to the near and far planes */
}  /* p3d_frust() */            /* of the viewing frustrum */

/*--------------------------------------------------------------------*/

void p3d_move (PROJ3D *p3d, double x, double y, double z)
{                               /* --- change eye position */
  if (x != 0) { p3d->curr.eye[0] += x *p3d->curr.cx[0];
                p3d->curr.eye[1] += x *p3d->curr.cx[1];
                p3d->curr.eye[2] += x *p3d->curr.cx[2]; }
  if (y != 0) { p3d->curr.eye[0] += y *p3d->curr.cy[0];
                p3d->curr.eye[1] += y *p3d->curr.cy[1];
                p3d->curr.eye[2] += y *p3d->curr.cy[2]; }
  if (z != 0) { p3d->curr.eye[0] += z *p3d->curr.cz[0];
                p3d->curr.eye[1] += z *p3d->curr.cz[1];
                p3d->curr.eye[2] += z *p3d->curr.cz[2]; }
}  /* p3d_move() */

/*--------------------------------------------------------------------*/

void p3d_yaw (PROJ3D *p3d, double phi, int mode)
{                               /* --- change heading angle */
  double mat[3][3];             /*     (rotate around vertical axis) */

  _rotmat(mat, p3d->curr.cz, phi);
  if (mode & P3D_VIEW)          /* rotate coordinate cross */
    _rotcrd(mat, p3d->curr.cz, p3d->curr.cx, p3d->curr.cy);
  if (mode & P3D_EYE)           /* rotate vector to the eye */
    _rotate(mat, p3d->curr.eye);
}  /* p3d_yaw() */

/*--------------------------------------------------------------------*/

void p3d_pitch (PROJ3D *p3d, double phi, int mode)
{                               /* --- change pitch angle */
  double mat[3][3];             /*     (rotate around horiz. axis) */

  _rotmat(mat, p3d->curr.cx, phi);
  if (mode & P3D_VIEW)          /* rotate coordinate cross */
    _rotcrd(mat, p3d->curr.cx, p3d->curr.cy, p3d->curr.cz);
  if (mode & P3D_EYE)           /* rotate vector to the eye */
    _rotate(mat, p3d->curr.eye);
}  /* p3d_pitch() */

/*--------------------------------------------------------------------*/

void p3d_roll (PROJ3D *p3d, double phi, int mode)
{                               /* --- change roll angle */
  double mat[3][3];             /*     (rotate around depth axis) */

  _rotmat(mat, p3d->curr.cy, phi);
  if (mode & P3D_VIEW)          /* rotate coordinate cross */
    _rotcrd(mat, p3d->curr.cy, p3d->curr.cz, p3d->curr.cx);
  if (mode & P3D_EYE)           /* rotate vector to the eye */
    _rotate(mat, p3d->curr.eye);
}  /* p3d_roll() */

/*--------------------------------------------------------------------*/

int p3d_proj (PROJ3D *p3d, double x, double y, double z,
                           double *px, double *py)
{                               /* --- project a point */
  double dx, dy, dz;            /* vector to object */
  double t;                     /* length of projection, factor */

  dx = x -p3d->curr.eye[0];     /* compute the vector from the eye */
  dy = y -p3d->curr.eye[1];     /* to the object and its projection */
  dz = z -p3d->curr.eye[2];     /* to the normal of the proj. plane */
  t  = dx*p3d->curr.cy[0] + dy*p3d->curr.cy[1] + dz*p3d->curr.cy[2];
  if ((t < p3d->curr.dnear) || (t > p3d->curr.dfar))
    return -1;                  /* if outside viewing frustrum, abort */
  t  = 1 /t;                    /* compute factor to object */
  dx = dx *t -p3d->curr.cy[0];  /* determine projection */
  dy = dy *t -p3d->curr.cy[1];  /* and the coordinates */
  dz = dz *t -p3d->curr.cy[2];  /* in the projection plane */
  *px = (dx*p3d->curr.cx[0] + dy*p3d->curr.cx[1] + dz*p3d->curr.cx[2])
      * p3d->curr.dpp;          /* compute x-coordinate of projection */
  *py = (dx*p3d->curr.cz[0] + dy*p3d->curr.cz[1] + dz*p3d->curr.cz[2])
      * p3d->curr.dpp;          /* compute y-coordinate of projection */
  return ((fabs(*px) > 1) || (fabs(*py) > 1))
       ? 1 : 0;                 /* return whether the point is */
}  /* p3d_proj() */             /* inside the viewing frustrum */

/*--------------------------------------------------------------------*/

int p3d_visib (PROJ3D *p3d, double x,  double y,  double z,
                            double nx, double ny, double nz)
{                               /* --- determine visibility */
  return ((x -p3d->curr.eye[0]) *nx  /* check normal vector against */
         +(y -p3d->curr.eye[1]) *ny  /* the vector from the eye */
         +(z -p3d->curr.eye[2]) *nz < 0) ? 1 : 0;
}  /* p3d_visib() */