insurfeval.cc

Mesa is an open-source implementation of the OpenGL specification - a system for rendering interacti

 * n: return normal, of dimension 3
 */
void OpenGLSurfaceEvaluator::inComputeNormal2(REAL *pu, REAL *pv, REAL *n)
{
  REAL mag; 

  n[0] = pu[1]*pv[2] - pu[2]*pv[1];
  n[1] = pu[2]*pv[0] - pu[0]*pv[2];
  n[2] = pu[0]*pv[1] - pu[1]*pv[0];  

  mag = sqrt(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);

  if (mag > 0.0) {
     n[0] /= mag; 
     n[1] /= mag;
     n[2] /= mag;
  }
}
 


/*Compute point and normal
 *see the head of inDoDomain2WithDerivs
 *for the meaning of the arguments
 */
void OpenGLSurfaceEvaluator::inDoEvalCoord2(REAL u, REAL v,
			   REAL *retPoint, REAL *retNormal)
{

  REAL du[4];
  REAL dv[4];

 
  assert(global_ev_k>=3 && global_ev_k <= 4);
  /*compute homegeneous point and partial derivatives*/
  inDoDomain2WithDerivs(global_ev_k, u, v, global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, retPoint, du, dv);

#ifdef AVOID_ZERO_NORMAL

  if(myabs(dv[0]) <= MYZERO && myabs(dv[1]) <= MYZERO && myabs(dv[2]) <= MYZERO)
    {

      REAL tempdu[4];
      REAL tempdata[4];
      REAL u1 = global_ev_u1;
      REAL u2 = global_ev_u2;
      if(u-MYDELTA*(u2-u1) < u1)
	u = u+ MYDELTA*(u2-u1);
      else
	u = u-MYDELTA*(u2-u1);
      inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, tempdu, dv);
    }
  if(myabs(du[0]) <= MYZERO && myabs(du[1]) <= MYZERO && myabs(du[2]) <= MYZERO)
    {
      REAL tempdv[4];
      REAL tempdata[4];
      REAL v1 = global_ev_v1;
      REAL v2 = global_ev_v2;
      if(v-MYDELTA*(v2-v1) < v1)
	v = v+ MYDELTA*(v2-v1);
      else
	v = v-MYDELTA*(v2-v1);
      inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, du, tempdv);
    }
#endif


  /*compute normal*/
  switch(global_ev_k){
  case 3:
    inComputeNormal2(du, dv, retNormal);

    break;
  case 4:
    inComputeFirstPartials(retPoint, du, dv);
    inComputeNormal2(du, dv, retNormal);
    /*transform the homegeneous coordinate of retPoint into inhomogenous one*/
    retPoint[0] /= retPoint[3];
    retPoint[1] /= retPoint[3];
    retPoint[2] /= retPoint[3];
    break;
  }
  /*output this vertex*/
/*  inMeshStreamInsert(global_ms, retPoint, retNormal);*/



  glNormal3fv(retNormal);
  glVertex3fv(retPoint);




  #ifdef DEBUG
  printf("vertex(%f,%f,%f)\n", retPoint[0],retPoint[1],retPoint[2]);
  #endif
  


}

/*Compute point and normal
 *see the head of inDoDomain2WithDerivs
 *for the meaning of the arguments
 */
void OpenGLSurfaceEvaluator::inDoEvalCoord2NOGE_BU(REAL u, REAL v,
			   REAL *retPoint, REAL *retNormal)
{

  REAL du[4];
  REAL dv[4];

 
  assert(global_ev_k>=3 && global_ev_k <= 4);
  /*compute homegeneous point and partial derivatives*/
//   inPreEvaluateBU(global_ev_k, global_ev_uorder, global_ev_vorder, (u-global_ev_u1)/(global_ev_u2-global_ev_u1), global_ev_ctlPoints);
  inDoDomain2WithDerivsBU(global_ev_k, u, v, global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, retPoint, du, dv);


#ifdef AVOID_ZERO_NORMAL

  if(myabs(dv[0]) <= MYZERO && myabs(dv[1]) <= MYZERO && myabs(dv[2]) <= MYZERO)
    {

      REAL tempdu[4];
      REAL tempdata[4];
      REAL u1 = global_ev_u1;
      REAL u2 = global_ev_u2;
      if(u-MYDELTA*(u2-u1) < u1)
	u = u+ MYDELTA*(u2-u1);
      else
	u = u-MYDELTA*(u2-u1);
      inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, tempdu, dv);
    }
  if(myabs(du[0]) <= MYZERO && myabs(du[1]) <= MYZERO && myabs(du[2]) <= MYZERO)
    {
      REAL tempdv[4];
      REAL tempdata[4];
      REAL v1 = global_ev_v1;
      REAL v2 = global_ev_v2;
      if(v-MYDELTA*(v2-v1) < v1)
	v = v+ MYDELTA*(v2-v1);
      else
	v = v-MYDELTA*(v2-v1);
      inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, du, tempdv);
    }
#endif

  /*compute normal*/
  switch(global_ev_k){
  case 3:
    inComputeNormal2(du, dv, retNormal);
    break;
  case 4:
    inComputeFirstPartials(retPoint, du, dv);
    inComputeNormal2(du, dv, retNormal);
    /*transform the homegeneous coordinate of retPoint into inhomogenous one*/
    retPoint[0] /= retPoint[3];
    retPoint[1] /= retPoint[3];
    retPoint[2] /= retPoint[3];
    break;
  }
}

/*Compute point and normal
 *see the head of inDoDomain2WithDerivs
 *for the meaning of the arguments
 */
void OpenGLSurfaceEvaluator::inDoEvalCoord2NOGE_BV(REAL u, REAL v,
			   REAL *retPoint, REAL *retNormal)
{

  REAL du[4];
  REAL dv[4];

 
  assert(global_ev_k>=3 && global_ev_k <= 4);
  /*compute homegeneous point and partial derivatives*/
//   inPreEvaluateBV(global_ev_k, global_ev_uorder, global_ev_vorder, (v-global_ev_v1)/(global_ev_v2-global_ev_v1), global_ev_ctlPoints);

  inDoDomain2WithDerivsBV(global_ev_k, u, v, global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, retPoint, du, dv);


#ifdef AVOID_ZERO_NORMAL

  if(myabs(dv[0]) <= MYZERO && myabs(dv[1]) <= MYZERO && myabs(dv[2]) <= MYZERO)
    {

      REAL tempdu[4];
      REAL tempdata[4];
      REAL u1 = global_ev_u1;
      REAL u2 = global_ev_u2;
      if(u-MYDELTA*(u2-u1) < u1)
	u = u+ MYDELTA*(u2-u1);
      else
	u = u-MYDELTA*(u2-u1);
      inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, tempdu, dv);
    }
  if(myabs(du[0]) <= MYZERO && myabs(du[1]) <= MYZERO && myabs(du[2]) <= MYZERO)
    {
      REAL tempdv[4];
      REAL tempdata[4];
      REAL v1 = global_ev_v1;
      REAL v2 = global_ev_v2;
      if(v-MYDELTA*(v2-v1) < v1)
	v = v+ MYDELTA*(v2-v1);
      else
	v = v-MYDELTA*(v2-v1);
      inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, du, tempdv);
    }
#endif

  /*compute normal*/
  switch(global_ev_k){
  case 3:
    inComputeNormal2(du, dv, retNormal);
    break;
  case 4:
    inComputeFirstPartials(retPoint, du, dv);
    inComputeNormal2(du, dv, retNormal);
    /*transform the homegeneous coordinate of retPoint into inhomogenous one*/
    retPoint[0] /= retPoint[3];
    retPoint[1] /= retPoint[3];
    retPoint[2] /= retPoint[3];
    break;
  }
}
 

/*Compute point and normal
 *see the head of inDoDomain2WithDerivs
 *for the meaning of the arguments
 */
void OpenGLSurfaceEvaluator::inDoEvalCoord2NOGE(REAL u, REAL v,
			   REAL *retPoint, REAL *retNormal)
{

  REAL du[4];
  REAL dv[4];

 
  assert(global_ev_k>=3 && global_ev_k <= 4);
  /*compute homegeneous point and partial derivatives*/
  inDoDomain2WithDerivs(global_ev_k, u, v, global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, retPoint, du, dv);


#ifdef AVOID_ZERO_NORMAL

  if(myabs(dv[0]) <= MYZERO && myabs(dv[1]) <= MYZERO && myabs(dv[2]) <= MYZERO)
    {

      REAL tempdu[4];
      REAL tempdata[4];
      REAL u1 = global_ev_u1;
      REAL u2 = global_ev_u2;
      if(u-MYDELTA*(u2-u1) < u1)
	u = u+ MYDELTA*(u2-u1);
      else
	u = u-MYDELTA*(u2-u1);
      inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, tempdu, dv);
    }
  if(myabs(du[0]) <= MYZERO && myabs(du[1]) <= MYZERO && myabs(du[2]) <= MYZERO)
    {
      REAL tempdv[4];
      REAL tempdata[4];
      REAL v1 = global_ev_v1;
      REAL v2 = global_ev_v2;
      if(v-MYDELTA*(v2-v1) < v1)
	v = v+ MYDELTA*(v2-v1);
      else
	v = v-MYDELTA*(v2-v1);
      inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, du, tempdv);
    }
#endif

  /*compute normal*/
  switch(global_ev_k){
  case 3:
    inComputeNormal2(du, dv, retNormal);
    break;
  case 4:
    inComputeFirstPartials(retPoint, du, dv);
    inComputeNormal2(du, dv, retNormal);
    /*transform the homegeneous coordinate of retPoint into inhomogenous one*/
    retPoint[0] /= retPoint[3];
    retPoint[1] /= retPoint[3];
    retPoint[2] /= retPoint[3];
    break;
  }
//  glNormal3fv(retNormal);
//  glVertex3fv(retPoint);
}
 
void OpenGLSurfaceEvaluator::inPreEvaluateBV(int k, int uorder, int vorder, REAL vprime, REAL *baseData)
{
  int j,row,col;
  REAL p, pdv;
  REAL *data;

  if(global_vprime != vprime || global_vorder != vorder) {      
    inPreEvaluateWithDeriv(vorder, vprime, global_vcoeff, global_vcoeffDeriv);
    global_vprime = vprime;
    global_vorder = vorder;
  }

  for(j=0; j<k; j++){
    data = baseData+j;
    for(row=0; row<uorder; row++){
      p = global_vcoeff[0] * (*data);
      pdv = global_vcoeffDeriv[0] * (*data);
      data += k;
      for(col = 1; col < vorder; col++){
	p += global_vcoeff[col] *  (*data);
	pdv += global_vcoeffDeriv[col] * (*data);
	data += k;
      }
      global_BV[row][j]  = p;
      global_PBV[row][j]  = pdv;
    }
  }
}

void OpenGLSurfaceEvaluator::inPreEvaluateBU(int k, int uorder, int vorder, REAL uprime, REAL *baseData)
{
  int j,row,col;
  REAL p, pdu;
  REAL *data;

  if(global_uprime != uprime || global_uorder != uorder) {      
    inPreEvaluateWithDeriv(uorder, uprime, global_ucoeff, global_ucoeffDeriv);
    global_uprime = uprime;
    global_uorder = uorder;
  }

  for(j=0; j<k; j++){
    data = baseData+j;
    for(col=0; col<vorder; col++){
      data = baseData+j + k*col;
      p = global_ucoeff[0] * (*data);
      pdu = global_ucoeffDeriv[0] * (*data);
      data += k*uorder;
      for(row = 1; row < uorder; row++){
	p += global_ucoeff[row] *  (*data);
	pdu += global_ucoeffDeriv[row] * (*data);
	data += k * uorder;
      }
      global_BU[col][j]  = p;
      global_PBU[col][j]  = pdu;
    }
  }
}
 
void OpenGLSurfaceEvaluator::inDoDomain2WithDerivsBU(int k, REAL u, REAL v,
						      REAL u1, REAL u2, int uorder,
						      REAL v1, REAL v2, int vorder,
						      REAL *baseData,
						      REAL *retPoint, REAL* retdu, REAL *retdv)
{
  int j, col;

  REAL vprime;


  if((u2 == u1) || (v2 == v1))
    return;

  vprime = (v - v1) / (v2 - v1);


  if(global_vprime != vprime || global_vorder != vorder) {
    inPreEvaluateWithDeriv(vorder, vprime, global_vcoeff, global_vcoeffDeriv);
    global_vprime = vprime;
    global_vorder = vorder;
  }


  for(j=0; j<k; j++)
    {
      retPoint[j] = retdu[j] = retdv[j] = 0.0;
      for (col = 0; col < vorder; col++)  {
	retPoint[j] += global_BU[col][j] * global_vcoeff[col];
	retdu[j] += global_PBU[col][j] * global_vcoeff[col];
	retdv[j] += global_BU[col][j] * global_vcoeffDeriv[col];
      }
    }
}    
   
void OpenGLSurfaceEvaluator::inDoDomain2WithDerivsBV(int k, REAL u, REAL v,
						      REAL u1, REAL u2, int uorder,
						      REAL v1, REAL v2, int vorder,
						      REAL *baseData,
						      REAL *retPoint, REAL* retdu, REAL *retdv)
{
  int j, row;
  REAL uprime;


  if((u2 == u1) || (v2 == v1))
    return;
  uprime = (u - u1) / (u2 - u1);


  if(global_uprime != uprime || global_uorder != uorder) {
    inPreEvaluateWithDeriv(uorder, uprime, global_ucoeff, global_ucoeffDeriv);
    global_uprime = uprime;
    global_uorder = uorder;
  }


  for(j=0; j<k; j++)
    {
      retPoint[j] = retdu[j] = retdv[j] = 0.0;
      for (row = 0; row < uorder; row++)  {
	retPoint[j] += global_BV[row][j] * global_ucoeff[row];
	retdu[j] += global_BV[row][j] * global_ucoeffDeriv[row];
	retdv[j] += global_PBV[row][j] * global_ucoeff[row];
      }
    }
}
  

/*
 *given a Bezier surface, and parameter (u,v), compute the point in the object space,
 *and the normal
 *k: the dimension of the object space: usually 2,3,or 4.
 *u,v: the paramter pair.
 *u1,u2,uorder: the Bezier polynomial of u coord is defined on [u1,u2] with order uorder.
 *v1,v2,vorder: the Bezier polynomial of v coord is defined on [v1,v2] with order vorder.
 *baseData: contrl points. arranged as: (u,v,k).
 *retPoint:  the computed point (one point) with dimension k.
 *retdu: the computed partial derivative with respect to u.
 *retdv: the computed partial derivative with respect to v.
 */
void OpenGLSurfaceEvaluator::inDoDomain2WithDerivs(int k, REAL u, REAL v, 
				REAL u1, REAL u2, int uorder, 
				REAL v1,  REAL v2, int vorder, 
				REAL *baseData,
				REAL *retPoint, REAL *retdu, REAL *retdv)
{
    int j, row, col;
    REAL uprime;
    REAL vprime;
    REAL p;
    REAL pdv;
    REAL *data;

    if((u2 == u1) || (v2 == v1))
	return;
    uprime = (u - u1) / (u2 - u1);
    vprime = (v - v1) / (v2 - v1);
    
    /* Compute coefficients for values and derivs */

    /* Use already cached values if possible */
    if(global_uprime != uprime || global_uorder != uorder) {
        inPreEvaluateWithDeriv(uorder, uprime, global_ucoeff, global_ucoeffDeriv);
	global_uorder = uorder;
	global_uprime = uprime;
    }
    if (global_vprime != vprime || 
	  global_vorder != vorder) {
	inPreEvaluateWithDeriv(vorder, vprime, global_vcoeff, global_vcoeffDeriv);
	global_vorder = vorder;
	global_vprime = vprime;
    }

    for (j = 0; j < k; j++) {
	data=baseData+j;
	retPoint[j] = retdu[j] = retdv[j] = 0.0;
	for (row = 0; row < uorder; row++)  {
	    /* 
	    ** Minor optimization.
	    ** The col == 0 part of the loop is extracted so we don't
	    ** have to initialize p and pdv to 0.
	    */
	    p = global_vcoeff[0] * (*data);
	    pdv = global_vcoeffDeriv[0] * (*data);
	    data += k;
	    for (col = 1; col < vorder; col++) {
		/* Incrementally build up p, pdv value */
		p += global_vcoeff[col] * (*data);
		pdv += global_vcoeffDeriv[col] * (*data);
		data += k;
	    }
	    /* Use p, pdv value to incrementally add up r, du, dv */
	    retPoint[j] += global_ucoeff[row] * p;
	    retdu[j] += global_ucoeffDeriv[row] * p;
	    retdv[j] += global_ucoeff[row] * pdv;
	}
    }  
}


/*
 *compute the Bezier polynomials C[n,j](v) for all j at v with 
 *return values stored in coeff[], where 
 *  C[n,j](v) = (n,j) * v^j * (1-v)^(n-j),
 *  j=0,1,2,...,n.
 *order : n+1
 *vprime: v
 *coeff : coeff[j]=C[n,j](v), this array store the returned values.
 *The algorithm is a recursive scheme:
 *   C[0,0]=1;
 *   C[n,j](v) = (1-v)*C[n-1,j](v) + v*C[n-1,j-1](v), n>=1
 *This code is copied from opengl/soft/so_eval.c:PreEvaluate
 */
void OpenGLSurfaceEvaluator::inPreEvaluate(int order, REAL vprime, REAL *coeff)
{
  int i, j;
  REAL oldval, temp;
  REAL oneMinusvprime;
  
  /*
   * Minor optimization
   * Compute orders 1 and 2 outright, and set coeff[0], coeff[1] to
     * their i==1 loop values to avoid the initialization and the i==1 loop.
     */