#pragma ident "@(#)normals.c 1.1 97/10/31 SMI"

/*
 * normals.c
 *
 *	@(#)normals.c 1.1 96/06/17 21:54:22
 *
 * Copyright (c) 1996, Sun Microsystems Computer Company, Mountain View, CA
 *
 *			All Rights Reserved
 *
 *	This is a program for manipulating 3D objects.  It is currently being
 *	developed by Mike M. Chow and Scott R. Nelson. 
 *
 *      This module contains methods for computing normals, derived from 
 *      leotool.
 */


#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include "types.h"
#include "zdebug.h"

/*
 * Structure and data for computing correct facet normals.
 */

#define MAX_VERTS_PER_FACET	256	/* Way more than we'll likely see */
#define INVALID_NORMAL		0x0001	/* This normal is invalid */
#define DIVERGENT_NORMS		0x0002	/* This normal diverges from others */
#define REJECTED_NORM		0x0004	/* This normal has been rejected */
/* How much a nomral may differ from the divergence of the others in a facet */
#define REJECTION_THRESHOLD	1.5	
/* How much a normal may diverge before being suspect (1.0 is no divergence) */
#define DIVERGENCE_THRESHOLD	0.900
/* How small the square root of the sum of a normal's values may be */
#define ROOT_THRESHOLD		0.000000000000001

typedef struct vert_normal_info vert_normal_info;

struct vert_normal_info {
    int flag;				/* Flag information */
    Point	normal;			/* The computed normal */
    float divergence;			/* How this differs from the rest */
};

static vert_normal_info vni[MAX_VERTS_PER_FACET];




/*
 * dot_prod
 *
 *	Compute the dot product of the two specified normals.
 *	The dot product is the cosine of the angle between the vectors.
 */

double
dot_prod(Point v0,
    Point v1)
{
    double dp;				/* The resulting dot product */

    dp = v0[0] * v1[0] + v0[1] * v1[1] + v0[2] * v1[2];

    return (dp);
} /* End of dot_prod */



/*
 * cross_prod
 *
 *	Compute the cross product of the three specified vectors.
 *	The cross product produces a vector normal to two input vectors.
 */

void
cross_prod(Point cp,
	   Point v0,
	   Point v1,
	   Point v2)
{
    
    register float x01, y01, z01;	/* v0[0] - v1[0], etc. */
    register float x21, y21, z21;	/* v2[0] - v1[0], etc. */

    x01 = v0[0] - v1[0];
    y01 = v0[1] - v1[1];
    z01 = v0[2] - v1[2];
    x21 = v2[0] - v1[0];
    y21 = v2[1] - v1[1];
    z21 = v2[2] - v1[2];

    cp[0] = (y01 * z21) - (z01 * y21);
    cp[1] = (z01 * x21) - (x01 * z21);
    cp[2] = (x01 * y21) - (y01 * x21);

} /* End of cross_prod */



/*
 * get_edge_angle
 *
 *	Compute the angle between two facets.  The dot product of
 *	the facet normals is the cosine of the angle between the
 *	normals.  Take the arccosine of this and get the angle.
 *	Convert this to degrees for easier handling.  Store the
 *	result in the edge_angle field.
 */

static void
get_edge_angle(Edge * eptr)
{
    double dp;				/* Dot product of normals */

/* Dot product is cosine of angle between vectors, get this angle */

    /* Get the dot product (cosine of angle) */
    dp = dot_prod(eptr->f1->normal, eptr->f2->normal);

    /* Convert to angle in degrees */
    eptr->edgeAngle = (acos(dp)) * (180.0 / M_PI);

} /* End of get_edge_angle */



void
pointFindNormal(Point p1, Point p2, Point p3, Point normal, int *flag)
{
  Point cp;			/* Cross product of three vertices */
  double dp;				/* Dot product of two normals */
  float root;				/* Used to normalize cross product */
  
  /* Compute the cross product of the two edges at this vertex */
  cross_prod(cp, p2, p1, p3);
  dp = dot_prod(cp, cp);	/* Prepare to normalize the value */
  root = sqrt(dp);
  if (root < ROOT_THRESHOLD) {
    /* The root is too small to use */
    *flag |= INVALID_NORMAL;
  }
  else {
    /* Compute normalization value, throw in direction reversal, too */
    root = -1.0 / root;
  }

  /* Store the (normalized) normal for this vertex */
  normal[0] = cp[0] * root;
  normal[1] = cp[1] * root;
  normal[2] = cp[2] * root;
}


/*
 * compute_facet_normal
 *
 *	Given the pointer to a facet description, compute normals
 *	at each vertex, then compute an average to use as the facet
 *	normal.  Mark any that differ significantly from the rest.
 *	Uses a global storage area shared by a later routine.
 *
 *	Note: Variable names may be cryptic, but at least they are
 *	well commented.
 */

static int
compute_facet_normal(ModelObject *obj, Facet * facet_ptr)
{
    int dc;				/* Divergent count */
    int v, p, n;			/* Indices for vertex, previous, next */
    int nv;				/* Number of vertices in facet */
    Point v0, v1, v2;		/* Three vertex pointers */
    Point cp;			/* Cross product of three vertices */
    double dp;				/* Dot product of two normals */
    float root;				/* Used to normalize cross product */
    Point fn;			/* The computed facet normal */
    int vnc;				/* Valid normal count */
    Vertex *varray;
    Point tmpNorm;

    dc = 0;
    varray = obj->varray;
    /* Compute a normal at each vertex of the facet */
    nv = facet_ptr->numVerts;
    for (v = 0; v < nv; v++) {
	vni[v].flag = 0;		/* Normal is good so far... */

	p = (v + nv - 1) % nv;		/* Index to previous vertex */
	n = (v + 1) % nv;		/* Index to next vertex */

	pointFindNormal(
			varray[facet_ptr->vindex[v]].p,
			varray[facet_ptr->vindex[p]].p,
			varray[facet_ptr->vindex[n]].p, 
			vni[v].normal, &vni[v].flag);
	
    } /* End of for each normal */

    /* See if any of these normals diverge greatly */
    for (p = 0; (vni[p].flag != 0) && (p < nv); p++)
	;				/* Use a good normal for comparisons */
    for (v = 0; v < nv; v++) {
	dp = dot_prod(vni[p].normal, vni[v].normal);
	if (dp < DIVERGENCE_THRESHOLD)
	    dc += 1;
	vni[v].divergence = dp;
    }

    /* Use a simple average to get the facet normal */
    /* Questionable normals will be fixed later */
    fn[0] = 0.0;
    fn[1] = 0.0;
    fn[2] = 0.0;
    vnc = 0;
    for (v = 0; v < nv; v++) {
	if (vni[v].flag != 0)
	    continue;			/* Skip all bad normals */
	fn[0] += vni[v].normal[0];
	fn[1] += vni[v].normal[1];
	fn[2] += vni[v].normal[2];
	vnc += 1;			/* Count another valid normal */
    }
    tmpNorm[0] = fn[0] / (float) vnc;
    tmpNorm[1] = fn[1] / (float) vnc;
    tmpNorm[2] = fn[2] / (float) vnc;

    /* Prepare to normalize the value */
    dp = dot_prod(tmpNorm, tmpNorm);
    root = sqrt(dp);
    if ((root < ROOT_THRESHOLD)) {
/*	fprintf(stderr, "Facet normal root is too small: %g\n", root);*/
    }
    else {
	root = 1.0 / root;
    }

    tmpNorm[0] = tmpNorm[0] * root;
    tmpNorm[1] = tmpNorm[1] * root;
    tmpNorm[2] = tmpNorm[2] * root;


    if (1) {
      int cntr;
      
      for (cntr = 0; cntr < 3; cntr++){
	if (isnan(facet_ptr->normal[cntr])){
	  zdo fprintf(stdout, "Facet normal[i] is NaN!\n");
	  facet_ptr->normal[0] = 0.0;
	  facet_ptr->normal[1] = 0.0;
	  facet_ptr->normal[2] = 1.0;
	}
	else
	      /* Store the (normalized) normal for this vertex */
	  facet_ptr->normal[cntr] = tmpNorm[cntr];	  
      }
    }



    return (dc);			/* Non-zero = problems */

} /* End of compute_facet_normal */



/*
 * get_facet_normal
 *
 *	Takes a pointer to a facet and sets the normal field in the facet.
 */

 void
facetFindNormal(ModelObject *obj, Facet * facet_ptr)
{
    int diverge_count;			/* How many normals were divergent */
    float div_av;			/* Average divergence */
    int i, j;				/* Iterative counters */
    int nv;				/* Number of vertices */
    double dp;				/* Dot product */
    Point fn;			/* The computed facet normal */
    int vnc;				/* Valid normal count */
    float root;				/* For normalization of the normal */
    Point tmpNorm;

    nv = facet_ptr->numVerts;
    diverge_count = compute_facet_normal(obj, facet_ptr);

    if (diverge_count > 0) {
	/* This one isn't flat enough, clean it up a bit */

	div_av = 0.0;			/* Find average divergence */
	for (i = 0; i < nv; i++) {
	    vni[i].divergence = 0.0;
	    for (j = 0; j < nv; j++) {
		dp = dot_prod(vni[i].normal, vni[j].normal);
		/* Dot product of identical vectors is 1.0 */
		vni[i].divergence += 1.0 - dp;
	    }
	    div_av += vni[i].divergence;
	}
	div_av /= (float) nv;
	div_av *=  REJECTION_THRESHOLD;	/* The threshold for rejection */
	for (i = 0; i < nv; i++) {
	    if (vni[i].divergence > div_av) {
		vni[i].flag |= DIVERGENT_NORMS;
	    }
	}

	/* Recompute the facet normal without problem vertices */
	fn[0] = 0.0;
	fn[1] = 0.0;
	fn[2] = 0.0;
	vnc = 0;
	for (i = 0; i < nv; i++) {
	    if (vni[i].flag != 0)
		continue;		/* Skip all bad normals */
	    fn[0] += vni[i].normal[0];
	    fn[1] += vni[i].normal[1];
	    fn[2] += vni[i].normal[2];
	    vnc += 1;			/* Count another valid normal */
	}
	tmpNorm[0] = fn[0] / (float) vnc;
	tmpNorm[1] = fn[1] / (float) vnc;
	tmpNorm[2] = fn[2] / (float) vnc;

	/* Prepare to normalize the value */
	dp = dot_prod(tmpNorm, tmpNorm);
	root = sqrt(dp);
	if ((root < ROOT_THRESHOLD)){
/*	    fprintf(stderr, "Facet normal root 2 is too small: %g\n", root);*/
	    tmpNorm[0] = 0.0;
	    tmpNorm[1] = 0.0;
	    tmpNorm[2] = 1.0;
	}
	else {
	    root = 1.0 / root;

	    /* Store the (normalized) normal for this vertex */
	    tmpNorm[0] = tmpNorm[0] * root;
	    tmpNorm[1] = tmpNorm[1] * root;
	    tmpNorm[2] = tmpNorm[2] * root;
	}
	if (1) {
	  int cntr;
	  
	  for (cntr = 0; cntr < 3; cntr++){
	    if (isnan(tmpNorm[cntr])){
	      zdo fprintf(stdout, "Facet normal[i] is NaN!\n");
	      facet_ptr->normal[0] = 0.0;
	      facet_ptr->normal[1] = 0.0;
	      facet_ptr->normal[2] = 1.0;
	    }
	    else 
	      facet_ptr->normal[cntr] = tmpNorm[cntr];
	  }
	}

	/* Try to fix concave facets to have better chances of looking right */
	for (i = 0; i < nv; i++) {
	    if (vni[i].flag != 0)
		break;			/* Find first divergent normal */
	}
	if (i < nv) {
	    VertIndex verts[MAX_VERTS_PER_FACET];

	    /*-
	     * If the vertex right after the problem vertex is first,
	     * triangulation will be reasonable.  This still doesn't always
	     * work if there is more than one concave vertex.
	     */
	    for (j = 0; j < nv; j++)
		verts[j] = facet_ptr->vindex[(j + i + 1) % nv];
	    for (j = 0; j < nv; j++)
		facet_ptr->vindex[j] = verts[j];
	}

    } /* End of cleaning up a facet with divergent normals */

#if 0
    /* Copy facet normal into vertex normals, for now */
    for (i = 0; i < facet_ptr->numVerts; i++)
	pointCopy(obj->varray[facet_ptr->vindex[i]].normal, facet_ptr->normal);
#endif

} /* End of facetFindNormal */



/*
 * vertexGetNormal
 *
 *	Find the average normal of all adjoining facet to this vertex.
 */

void
vertexGetNormal(Vertex *vert)
{
  FacetList *flistPtr;       		/* Pointer to the current facet */
  int numFacets;			/* # facets sharing the point */
  float x, y, z;			/* The normal value */
  float root;				/* For normalizing the normal */
  
  flistPtr = vert->facets;
  numFacets = 0;
  x = 0.0;
  y = 0.0;
  z = 0.0;

  /* Go through each facet of the vertex to sum up their normals. */
  while (flistPtr != NULL) {
    if (flistPtr->current != NULL){
      x += flistPtr->current->normal[0];
      y += flistPtr->current->normal[1];
      z += flistPtr->current->normal[2];
      numFacets++;
    }
    flistPtr = flistPtr->next;
  }
  
  if (numFacets == 0) {
    /* Set to default value */
    vert->normal[0] = 0.0;
    vert->normal[1] = 0.0;
    vert->normal[2] = 1.0;
  }
  else {
    int i;

    x /= (float) numFacets;
    y /= (float) numFacets;
    z /= (float) numFacets;
    
    root = sqrt(x * x + y * y + z * z);
 
    if (root > 0.0)
      root = 1.0 / root;

    vert->normal[0] = x * root;
    vert->normal[1] = y * root;
    vert->normal[2] = z * root;
    
    for(i = 0;i < 3; i++){
      if (fabs(vert->normal[i]) > 1.0) {
	/* Set to default value */
	vert->normal[0] = x;
	vert->normal[1] = y;
	vert->normal[2] = z;
	break;
      }
    }
    
  }
  vert->flags |= TYPES_VTX_NORM;

} /* End of vertexGetNormal */





void
objGetVertexNormals (ModelObject *obj)
{
  int i;

  for(i = 0; i < obj->vertCnt; i++){
    vertexGetNormal((obj->varray)+i);
  }
}
 /* End of normals.c*/
