#pragma ident "@(#)varyQuant.c	1.2 98/08/25 SMI"

/*
 * varyQuant.c
 *
 *	@(#)varyQuant.c 1.8 96/11/04 
 *
 * Copyright (c) 1997 by Mike M. Chow - All Rights Reserved 
 *
 *
 *      Author: Mike M. Chow 
 *
 *      This module contains quantization assignments for compression.
 * 
 */

#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <sys/time.h>
#include "types.h"
#include "linalg.h"
#include "varyQuant.h"
#include "zdebug.h"
#include "utils.h"


#define zoff(statements) \
if (1)\
{int oldzzdebug; \
 oldzzdebug = zzdebug; \
 zzdebug = 0; \
 statements \
 zzdebug = oldzzdebug; \
			 } 


#define MIN_MESH_REGIONS     10
#if 0
/* To show smaller regions */
#define MIN_MESH_REGION_AREA 20
#define MAX_MESH_REGION_AREA 50
#endif

#define MIN_MESH_REGION_AREA 45
#define MAX_MESH_REGION_AREA 200
#define VerySmallRegionCutoff 10

static int MinMeshRegionCutoff;

static float VaryQuantError = 0.2;
static int VaryQuantLow = 7;
static int VaryQuantHigh = 15;
static int VaryQuant = 0;

typedef struct Norm2Quant Norm2Quant;

struct Norm2Quant {
  float low;
  float high;
  int   quant;
};

static Norm2Quant Normal2QuantTable[10] = {{-1.0, 10.0, 6}, 
					   {10.0, 35.0, 7},
					   {35.0, 50.0, 8},
					   {50.0, 70.0, 9},
					   {70.0, 100.0, 10},
					   {100.0, 700000, 11}};


int MaxQuant;


MeshRegion *
meshRegionCreate()
{
  MeshRegion *region;

  region = (MeshRegion *) malloc(sizeof(MeshRegion));
  region->flist = linkedListNewNode();
  region->refFacet = NULL;
  region->facetCnt = 0;
  region->adjRegions = linkedListNewNode();
  region->maxNormalChange = 0.0;
  region->quantLevel = 16.0;
  region->bbox = bboxCreate();
  region->maxSizeRatio = 1.0;
  region->boundaryRatio = 0.0;
  region->meshInfo = NULL;
  region->cmpSize = 0;
  region->cmpBuff = NULL;
  return region;
}


void
meshRegionDestroy(MeshRegion *region)
{
  linkedListDestroy(region->flist);
  meshInfoDestroy(region->meshInfo);
  free(region);
}


void
meshRegionAddFacet(MeshRegion *region, Facet *fptr)
{
  region->facetCnt++;
  linkedListAddNode(&region->flist, (void *) fptr);
  fptr->region = region;
}

void
meshRegionRemoveFacet(MeshRegion *region, Facet *fptr)
{
  linkedListRemoveNode(&region->flist, (void *) fptr);
  region->facetCnt--;
  fptr->region = NULL;
}


void
meshRegionFindBbox(MeshRegion *region)
{
  LinkedList *lptr;
  Facet *fptr;
  int i;

  for (lptr = region->flist; lptr != NULL; lptr = lptr->next){
    fptr = (Facet *) lptr->current;

    if (fptr != NULL){
      /* Extends bbox */
      for (i = 0; i < fptr->numVerts; i++) {
	Vertex *vert;
	/* Get vertex position and normal */
	vert = &region->owner->obj->varray[fptr->vindex[i]];
	bboxExtendPoint(region->bbox, vert->p);
      } /* End of for each vertex in the facet */
    }
  }
  for(i = 0; i < 3 ;i++){
    region->bboxSize[i] = region->bbox->max[i] - region->bbox->min[i];
  }
}




int
meshRegionCheck(void *node, void *data)
{
  MeshRegion *region = (MeshRegion *) node;
  LinkedList *lptr;


  zdo6 printf("Sanity checking region %x...\n", region);

  for (lptr = region->flist; lptr != NULL; lptr = lptr->next){
    Facet *fptr = (Facet *) lptr->current;
    zdo6 printf("facet %x region %x\n",
	       fptr, fptr->region);
    if (fptr->region != region) {
      printf("err!!! facet %x 's region %x is NOT owner region %x\n",
	     fptr, fptr->region, region);  
      exit(1);
    }
    
  }
  return 1;
}


void
mergeMeshRegion(MeshRegion *src, MeshRegion *dest)
{
  LinkedList *lptr;
  
  zdo6 printf("Merging region %x with dest %x \n", 
	     src, dest);

#if 0
  /* For debugging */

  for (lptr = dest->flist; lptr != NULL; lptr = lptr->next){
    Facet *facet;

    facet = (Facet *) lptr->current;
    if (facet != NULL){
      zdo6 printf("dest facet %x region owner %x\n",
		 facet, facet->region);
    }
  }
#endif

  for (lptr = dest->flist; lptr != NULL; lptr = lptr->next){
    Facet *facet;

    facet = (Facet *) lptr->current;
    if (facet != NULL){
      zdo6 printf("adding facet %x facet's old region %x, new region %x\n", facet,
		 facet->region, src);
      meshRegionAddFacet(src, facet);
      if (facet->boundaryFacet) {
	int i;
	Vertex *vert;

	facet->boundaryFacet = 0;
	for (i = 0; i < facet->numVerts; i++){
	  LinkedList *elptr;

	  vert = &src->owner->obj->varray[facet->vindex[i]];
	  vert->boundaryVertex = 0;

	  /* Check for boundary vertices */
	  
	  zdo6 printf("checking vert %x for hard edge. \n", vert);
	  for (elptr = vert->elist; elptr != NULL; elptr = elptr->next){
	    Edge *eptr;

	    eptr = (void *) elptr->current;
	    zdo6 printf("looking at vert edge %x hard %d\n", eptr, eptr->hard);
	    if (eptr->hard){
	      vert->boundaryVertex = 1;
	      break;
	    }
	  }
	  
	}
      }
      free(facet->pickedColor);
      facet->pickedColor = src->refFacet->pickedColor;
    }
  }
}



int
meshRegionPrint(void *node, void *data)
{
  MeshRegion *region = (MeshRegion *) node;

  printf("Region %x facets %d maxNCh %f maxSRatio %f boundRatio %f quant %d\n",
	 region, region->facetCnt, region->maxNormalChange, region->maxSizeRatio,
	 region->boundaryRatio, region->quantLevel);
  return 1;
}



int
meshRegionSetRegion(void *node, void *data)
{
  MeshRegion *region = (MeshRegion *) node;
  LinkedList *lptr;

  for (lptr = region->flist; lptr != NULL; lptr = lptr->next){
    Facet *fptr = (Facet *) lptr->current;
    fptr->region = region;    
  }
  return 1;
}



int
nukeOldAdjRegions(void *node, void *data)
{
  MeshRegion *region = (MeshRegion *) node;
  LinkedList *regions = (LinkedList *) data;
  LinkedList *lptr, *lptr2, *next;


  zdo6 printf("Nuke checking region %x...\n", region);

  for (lptr = region->adjRegions; lptr != NULL; lptr = next){
    
    MeshRegion *nbr = (MeshRegion *) lptr->current;
    int found;

    next = lptr->next;
    if (nbr != NULL) {
      found = 0;
      zdo6 printf("checking adj region %x..\n", nbr);
      for (lptr2 = regions; lptr2 != NULL; lptr2 = lptr2->next){
	MeshRegion *tmp = (MeshRegion *) lptr2->current;
	
	if (tmp == nbr){
	  found = 1;
	  zdo6 printf("found in part region %x.\n", tmp);
	}
      }
      if (!found){
	zdo6 printf("Bad region %x found not in part's regions! Nuking from adj region \n",
	       nbr);
	linkedListRemoveNode(&(region->adjRegions), nbr);
      }
    }
  }
  return 1;
}


int
testAddAdjRegion(Facet *adjFacet, Edge *edge, void *data)
{
  Facet *facet = (Facet *) data;
  
  if (adjFacet->region != facet->region){
    zdo6 printf("Adding adj region %x\n", adjFacet->region);
    facet->boundaryFacet = 1;
    adjFacet->boundaryFacet = 1;
    edge->v1->boundaryVertex = 1;
    edge->v2->boundaryVertex = 1;
    linkedListAddNode(&facet->region->adjRegions, (void *) adjFacet->region);
  }
  return 1;
}




 
int
meshRegionFindAdjRegions(void *node, void *data)
{
  LinkedList *lptr;
  MeshRegion *region = (MeshRegion *) node;

  zdo6 printf("Find adj regions of region %x\n", region);
  for (lptr = region->flist; lptr != NULL; lptr = lptr->next){
    Facet *facet;

    facet = (Facet *) lptr->current;
    if (facet != NULL){
      facetForEachAdjFacet(facet, testAddAdjRegion, (void *) facet);
    }
  }
  return 1;
}



int
printAdjRegions(MeshRegion *region)
{
  MeshRegion *adjRegion;
  LinkedList *rlptr2;

  zdo printf("Adj regions for %x are:\n", region); 
  /* See if any adjacent regions are mergeable */
  for (rlptr2 = region->adjRegions; rlptr2 != NULL; 
       rlptr2 = rlptr2->next)
    {
      adjRegion = (MeshRegion *) rlptr2->current;
      if (adjRegion != NULL){
	meshRegionPrint((void *) adjRegion, NULL);
      }
    }
  return 1;
}
      

int
mergeAdjFacets(Facet *adjFptr, Edge *adjEdge, void *data)
{
  MeshRegion *region = (MeshRegion *) data;
  Facet *fptr = region->refFacet;
  float diff;
  float dot, theta;


  /* Must not have a region already! */
  if (adjFptr->region != NULL)
    return 1;

  diff = fabs(adjFptr->radius - fptr->radius);
  
  /* If adjacent facet's radius is within 15% of fptr's radius
     then continue. */ 
  /*
    (diff < (0.25 * (adjFptr->radius + fptr->radius)/2.0))
    */
  dot = pointDot(adjFptr->normal, fptr->normal);
  theta = fabs(acos(dot)) * 180.0/M_PI;

  zdo6 printf("adj fptr rad %f diff %f \n", adjFptr->radius, diff);
  zdo6 printf("adj fptr theta %f, time stamp %ld \n", theta, adjFptr->timeStamp);
  if ((adjFptr->timeStamp != fptr->timeStamp) && 
      (theta < 10.0)) {
    zdo6 printf("adj fptr %x close, cont merge..\n", adjFptr);
    adjFptr->timeStamp = fptr->timeStamp;
    adjFptr->selected = 1;
    adjFptr->pickedColor = rgbaCreate(fptr->pickedColor);
    meshRegionAddFacet(region, adjFptr);
    facetForEachAdjFacet(adjFptr, mergeAdjFacets, (void *) region);
  }
  else 
    zdo6 printf("adj not close, cont\n");

  return 1;
}



int
normalChange2quant(float maxNChange, Point bboxSize, float *maxSizeRatio)
{
  int i;
  int quantLevel = 12;

  for (i = 0; i < 10; i++){
    if ((maxNChange <= Normal2QuantTable[i].high) &&
	(maxNChange >= Normal2QuantTable[i].low)){
      quantLevel = Normal2QuantTable[i].quant;
      break;
    }
  }
  
  if (quantLevel < 8){
    float xy, yz, xz, maxRatio;
    
    maxRatio = 0.0;
    /* Danger */

    xy = bboxSize[0]/bboxSize[1];
    if (xy < 1.0) 
      xy = 1.0/xy;
    xz = bboxSize[0]/bboxSize[2];
    if (xz < 1.0) 
      xz = 1.0/xz;
    yz = bboxSize[1]/bboxSize[2];
    if (yz < 1.0) 
      yz = 1.0/yz;
    if (xy > xz)
      maxRatio = xy;
    if (yz > maxRatio)
      maxRatio = yz;

    zdo6 printf("Danger: region < 8 bit, maxRatio %f\n",
	   maxRatio);
    if (maxRatio < 6.0){
      zdo6 printf("ratio too low, clamp to 8\n"); 
      quantLevel = 8;
    }
    *maxSizeRatio = maxRatio;
  }

  return quantLevel;
}
    
    


float
findMaxNChange(LinkedList *nodes, Facet *refFacet, ModelObject *obj)
{
  Facet *fptr, *maxDeltaFacet;
  float maxDeltaTheta;
  LinkedList *lptr, *lptr2;
  float deltaSum = 0.0;
  int totalDeltas = 0;
  int i;
  
  maxDeltaTheta = 0.0;
  maxDeltaFacet = NULL;

  for (lptr2 = nodes; lptr2 != NULL; lptr2 = lptr2->next){
    Facet *fptr2;

    fptr2 = (Facet *) lptr2->current;

    /* Check only 2/3 of the facets */
    if (randNum(0.0, 1.0) < 0.3)
      continue;
    for (lptr = nodes; lptr != NULL; lptr = lptr->next){
      fptr = (Facet *) lptr->current;
      
      if (fptr != NULL){
	float theta, dot;
	dot = pointDot(fptr2->normal, fptr->normal);
	theta = fabs(acos(dot)) * 180.0/M_PI;
	zdo6 printf("theta %f \n", theta);
	deltaSum += theta;
	totalDeltas++;
	if (theta > maxDeltaTheta) {
	  maxDeltaFacet = fptr;
	  maxDeltaTheta = theta;
	}
      }
    }
  }

  return maxDeltaTheta;
}
  

int
assignRegionQuant(void *node, void *data)
{
  MeshRegion *region = (MeshRegion *) node;
  Facet *fptr, *maxDeltaFacet;
  float maxDeltaTheta;
  LinkedList *lptr;
  float deltaSum = 0.0;
  int totalDeltas = 0;
  int i;
  
  maxDeltaTheta = 0.0;
  maxDeltaFacet = NULL;

  for (lptr = region->flist; lptr != NULL; lptr = lptr->next){
    fptr = (Facet *) lptr->current;

    if (fptr != NULL){
      float theta, dot;

      /* Extends bbox */
      for (i = 0; i < fptr->numVerts; i++) {
	Vertex *vert;
	/* Get vertex position and normal */
	vert = &region->owner->obj->varray[fptr->vindex[i]];
	bboxExtendPoint(region->bbox, vert->p);
      } /* End of for each vertex in the facet */

      dot = pointDot(region->refFacet->normal, fptr->normal);
      theta = fabs(acos(dot)) * 180.0/M_PI;
      zdo6 printf("theta %f \n", theta);
      deltaSum += theta;
      totalDeltas++;
      if (theta > maxDeltaTheta) {
	maxDeltaFacet = fptr;
	maxDeltaTheta = theta;
      }
    }
  }
    
  for(i = 0; i < 3 ;i++){
    region->bboxSize[i] = region->bbox->max[i] - region->bbox->min[i];
  }

  zdo6 printf("Max delta theta: %f\n", maxDeltaTheta);
  region->maxNormalChange = deltaSum / (float) totalDeltas;
  region->quantLevel = normalChange2quant(region->maxNormalChange,
					  region->bboxSize,
					  &region->maxSizeRatio);
  /* Snap to Max Quant level assigned earlier */
  if (region->quantLevel > MaxQuant)
    region->quantLevel = MaxQuant;

  /* Snap to user specified levels */
  if (region->quantLevel < VaryQuantLow)
    region->quantLevel = VaryQuantLow;
  if (region->quantLevel > VaryQuantHigh)
    region->quantLevel = VaryQuantHigh;
  return region->quantLevel;
}
 


int
regionIsSmall(MeshRegion *region)
{
  int regionInteriorVertices = 0, regionBoundaryVertices = 0;
  ModelObject *obj;
  LinkedList *lptr;
  int i;
  float boundaryRatio;

  /* Fast small region test */
  if (region->facetCnt < MinMeshRegionCutoff)
    return 1;

  /* Now, for the more expensive, boundary:interior ratio test */
  obj = region->owner->obj;
  for (lptr = region->flist; lptr != NULL; lptr = lptr->next)
    {
      Facet *fptr;

      fptr = (Facet *) lptr->current;

      if (fptr != NULL) {
	  for (i = 0; i < fptr->numVerts; i++) {
	    Vertex *vert;

	    /* Get vertex position and normal */
	    vert = &obj->varray[fptr->vindex[i]];
	    if (!vert->boundaryVertex){
	      regionInteriorVertices++;
	    }
	    else {
	      regionBoundaryVertices++;
	    }
	  }
      }
    }

  boundaryRatio = (float)regionBoundaryVertices 
    / (float)(regionBoundaryVertices + regionInteriorVertices);

  region->boundaryRatio = boundaryRatio;
  if (boundaryRatio > 0.99){
    zdo6 printf("boundaryRatio %f, small!\n", 
	       boundaryRatio);
    return 1;
  }
  return 0;
}
  
  
void
unlinkAdjRegions(LinkedList *allRegions, MeshRegion *killedRegion)
{
  LinkedList *lptr, *lptr2, *nextLink;

  for (lptr2 = allRegions; lptr2 != NULL; lptr2 = lptr2->next){
    MeshRegion *region = lptr2->current;
    
    for (lptr = region->adjRegions; lptr != NULL; lptr = nextLink){
      MeshRegion *adjRegion = (MeshRegion *) lptr->current;
      nextLink = lptr->next;

      if (adjRegion != NULL) 
	linkedListRemoveNode(&region->adjRegions, (void *) killedRegion);
    }
  }
}

  


int
regionsMergeable(MeshRegion *r1, MeshRegion *r2, int quantThreshold)
{
  Facet *fptr, *maxDeltaFacet;
  float maxDeltaTheta;
  LinkedList *lptr;
  float deltaSum = 0.0;
  int totalDeltas = 0;
  int i;
  Point bboxSize;
  BBox *combinedBbox;
  Facet *refFacet = r1->refFacet;
  int maxNormalChange;
  int quantLevel;
  float maxSizeRatio;
  int diff1, diff2;

  maxDeltaTheta = 0.0;
  maxDeltaFacet = NULL;

  combinedBbox = bboxCreate();
  bboxCopy(combinedBbox, r1->bbox);
  for (lptr = r2->flist; lptr != NULL; lptr = lptr->next){
    fptr = (Facet *) lptr->current;

    if (fptr != NULL){
      float theta, dot;

      /* Extends bbox */
      for (i = 0; i < fptr->numVerts; i++) {
	Vertex *vert;
	/* Get vertex position and normal */
	vert = &r2->owner->obj->varray[fptr->vindex[i]];
	bboxExtendPoint(combinedBbox, vert->p);
      } /* End of for each vertex in the facet */

      dot = pointDot(refFacet->normal, fptr->normal);
      theta = fabs(acos(dot)) * 180.0/M_PI;
      zdo6 printf("theta %f \n", theta);
      deltaSum += theta;
      totalDeltas++;
      if (theta > maxDeltaTheta) {
	maxDeltaFacet = fptr;
	maxDeltaTheta = theta;
      }
    }
  }
  
  for(i = 0; i < 3 ;i++){
    bboxSize[i] = combinedBbox->max[i] - combinedBbox->min[i];
  }

  zdo6 printf("Max delta theta: %f\n", maxDeltaTheta);
  maxNormalChange = deltaSum / (float) totalDeltas;
  quantLevel = normalChange2quant(maxNormalChange,
				  bboxSize,
				  &maxSizeRatio);

  zdo6 printf("r1 %d r2 %d Resulting quantLevel %d \n",
	     r1->quantLevel, r2->quantLevel, quantLevel);
  diff1 = (int) fabs((float)(quantLevel - r1->quantLevel));
  diff2 = (int) fabs((float )(quantLevel - r2->quantLevel));

  if ((diff1 > quantThreshold) || (diff2 > quantThreshold)){
    zdo6 printf("Large Regions %x and %x too different (%d %d), not merging\n",
	   r1, r2, diff1, diff2);
    return 0;
  }
  else {
    zdo6 printf("Merge Large Regions %x and %x (%d %d)\n",
	   r1, r2, diff1, diff2);
    return 1;
  }
}


void
mergeSmallRegions(LinkedList **regions)
{
  LinkedList *rlptr;
  MeshRegion *region;


  zdo printf("Merging small regions... \n");

  for(rlptr = *regions; rlptr != NULL; rlptr = rlptr->next){
    region = (MeshRegion *) rlptr->current;
    if (region != NULL){
      if (regionIsSmall(region)) {
	int change;
	int iterations = 0;

	zdo printf("Region %x has area %d, more than cutoff, start merging...\n\n",
		   region, region->facetCnt);

	change = 1;


	while ((change) &&
	       (regionIsSmall(region))) {
	  LinkedList *rlptr2, *nextLink;
	  MeshRegion *adjRegion;
	  
	  change = 0;
	  zdo printf("\nIteration %d maxNormalChange %f region facetCnt %d\n",
		     iterations, region->maxNormalChange, region->facetCnt);
	  iterations++;
	  zdo printAdjRegions(region);	  
	  /* See if any adjacent regions are mergeable */
	  for (rlptr2 = region->adjRegions; rlptr2 != NULL; 
	       rlptr2 = nextLink){
	    nextLink = rlptr2->next;
	    adjRegion = (MeshRegion *) rlptr2->current;
	    if ((adjRegion != NULL) &&
		(1)){
	      float diff;
	      float dot, theta;

	      diff = fabs(adjRegion->maxNormalChange -
			  region->maxNormalChange);

	      dot = pointDot(adjRegion->refFacet->normal, 
			     region->refFacet->normal);
	      theta = fabs(acos(dot)) * 180.0/M_PI;

	      zdo printf("Looking at Adj region %x, maxNchange %f facets %d\n", adjRegion,
			   adjRegion->maxNormalChange,
			   adjRegion->facetCnt); 

	      /*(diff < 3.0) ||  (region->maxNormalChange <
			  adjRegion->maxNormalChange)
			  */
	      /* Merge? */
	      if (regionIsSmall(adjRegion) && 
		  (regionsMergeable(region, adjRegion, 3))){
		LinkedList *lptr;

		zdo printf("Adj region %x, maxNchange %f facets %d merge.\n", adjRegion,
			   adjRegion->maxNormalChange,
			   adjRegion->facetCnt); 

		mergeMeshRegion(region, adjRegion);
		linkedListRemoveNode(regions, (void *) adjRegion);
		linkedListRemoveNode(&region->adjRegions, (void *) adjRegion);

		/* Unlink those neighbor regions that border adjRegion */ 
		unlinkAdjRegions(*regions, adjRegion);

		/* Check sanity */
		linkedListForEachNode(*regions, meshRegionCheck, NULL);
		change = 1;
	      }
	    }
	  }

	  /* Recompute adjacent regions of R by getting rid of old adj regions 
	   * R. 
	   */
	  linkedListDestroy(region->adjRegions);
	  region->adjRegions = linkedListNewNode();
	  zdo printf("Recomputing adjacent regions...\n");
	  meshRegionFindAdjRegions(region, NULL);

	  assignRegionQuant((void *) region, NULL);
	}
      }
    } 
  }
}




 
void
mergeLargeWithLargeRegions(LinkedList **regions)
{
  LinkedList *rlptr;
  MeshRegion *region;
  
  zdo printf("Merging Large wth Large regions... \n");

  for(rlptr = *regions; rlptr != NULL; rlptr = rlptr->next){
    region = (MeshRegion *) rlptr->current;
    if (region != NULL){
      if (region->facetCnt > MinMeshRegionCutoff) {
	int change;
	int iterations = 0;

	zdo printf("Large Region %x has area %d, more than cutoff, start merging...\n\n",
		   region, region->facetCnt);

	change = 1;


	while ((change)
	       ) {
	  LinkedList *rlptr2, *nextLink;
	  MeshRegion *adjRegion;
	  
	  zdo6 printAdjRegions(region);
	  change = 0;
	  zdo printf("\nIteration %d maxNormalChange %f region facetCnt %d\n",
		     iterations, region->maxNormalChange, region->facetCnt);
	  iterations++;
	  
	  /* See if any adjacent regions are mergeable */
	  for (rlptr2 = region->adjRegions; rlptr2 != NULL; 
	       rlptr2 = nextLink){
	    nextLink = rlptr2->next;
	    adjRegion = (MeshRegion *) rlptr2->current;
	    if ((adjRegion != NULL) &&
		(adjRegion->facetCnt > MinMeshRegionCutoff)){
	      float diff;

	      diff = fabs((float) (adjRegion->quantLevel -
				   region->quantLevel));

	      /* Merge if diff is small or large region is worse than smaller */
	      if (regionsMergeable(region, adjRegion, 2)){

		zdo printf("Adj region %x, quant diff %f facets %d merge.\n",
			   adjRegion,
			   diff,
			   adjRegion->facetCnt); 
		/* Check sanity */
		zdo6 linkedListForEachNode(*regions, meshRegionCheck, NULL);
		zdo printf("sanity check before merging.\n");
		mergeMeshRegion(region, adjRegion);
		linkedListRemoveNode(regions, (void *) adjRegion);
		linkedListRemoveNode(&region->adjRegions, (void *) adjRegion);
		/* Unlink those neighbor regions that border adjRegion */ 
		unlinkAdjRegions(*regions, adjRegion);

		/* Check sanity */
		linkedListForEachNode(*regions, meshRegionCheck, NULL);

		change = 1;
	      }
	    }
	  }

	  /* Recompute adjacent regions of R by getting rid of old adj regions 
	   * R. 
	   */
	  linkedListDestroy(region->adjRegions);
	  region->adjRegions = linkedListNewNode();
	  zdo printf("Recomputing adjacent regions...\n");
	  meshRegionFindAdjRegions(region, NULL);

	  assignRegionQuant((void *) region, NULL);
	}
      }
    } 
  }
}



void
mergeLargeWithSmallRegions(LinkedList **regions)
{
  LinkedList *rlptr;
  MeshRegion *region;
  
  zdo printf("Merging Large & Small regions... \n");

  for(rlptr = *regions; rlptr != NULL; rlptr = rlptr->next){
    region = (MeshRegion *) rlptr->current;
    if (region != NULL){
      if (region->facetCnt > MinMeshRegionCutoff) {
	int change;
	int iterations = 0;

	zdo printf("Region %x has area %d, more than cutoff, start merging...\n\n",
		   region, region->facetCnt);

	change = 1;

	while ((change)
	       ) {
	  LinkedList *rlptr2, *nextLink;
	  MeshRegion *adjRegion;

	  zdo6 printAdjRegions(region);
	  change = 0;
	  zdo printf("\nIteration %d maxNormalChange %f region facetCnt %d\n",
		     iterations, region->maxNormalChange, region->facetCnt);
	  iterations++;
	  
	  /* See if any adjacent regions are mergeable */
	  for (rlptr2 = region->adjRegions; rlptr2 != NULL; 
	       rlptr2 = nextLink){
	    nextLink = rlptr2->next;
	    adjRegion = (MeshRegion *) rlptr2->current;
	    if ((adjRegion != NULL) &&
		(regionIsSmall(adjRegion))){
	      float diff;

	      diff = fabs((float) (adjRegion->quantLevel -
				   region->quantLevel));


	      zdo printf("Looking Adj region %x, quant diff %f facets %d.\n",
			 adjRegion,
			 diff,
			 adjRegion->facetCnt); 

	      /* Merge if diff is small or large region is worse than smaller */
	      if ((region->quantLevel > adjRegion->quantLevel)
		  || regionsMergeable(region, adjRegion, 2))
		  {
		    zdo printf("Adj region is mergeable. \n");
		/* Check sanity */
		zdo linkedListForEachNode(*regions, meshRegionCheck, NULL);
		zdo printf("Sanity checked before merging.\n");
		mergeMeshRegion(region, adjRegion);
		linkedListRemoveNode(regions, (void *) adjRegion);

		/* Check sanity */
		linkedListForEachNode(*regions, meshRegionCheck, NULL);

		linkedListRemoveNode(&region->adjRegions, (void *) adjRegion);
		/* Unlink those neighbor regions that border adjRegion */ 
		unlinkAdjRegions(*regions, adjRegion);


		change = 1;
	      }
	    }
	  }

	  /* Recompute adjacent regions of R by getting rid of old adj regions 
	   * R. 
	   */
	  linkedListDestroy(region->adjRegions);
	  region->adjRegions = linkedListNewNode();
	  zdo printf("Recomputing adjacent regions...\n");
	  meshRegionFindAdjRegions(region, NULL);

	  assignRegionQuant((void *) region, NULL);
	}
      }
    } 
  }
}


void
mergeDisjointedRegions(LinkedList **regions)
{
  LinkedList *rlptr;
  MeshRegion *region;


  zdo printf("Merging disjointed regions... \n");

  for(rlptr = *regions; rlptr != NULL; rlptr = rlptr->next){
    region = (MeshRegion *) rlptr->current;
    if (region != NULL){
      if (1) {
	int change;
	int iterations = 0;

	zdo printf("Region %x has area %d, more than cutoff, start merging...\n\n",
		   region, region->facetCnt);

	change = 1;


	while ((change) &&
	       (regionIsSmall(region))) {
	  LinkedList *rlptr2, *nextLink;
	  MeshRegion *region2;
	  
	  change = 0;
	  zdo printf("\nIteration %d maxNormalChange %f region facetCnt %d\n",
		     iterations, region->maxNormalChange, region->facetCnt);
	  iterations++;
	  /* See if any adjacent regions are mergeable */
	  for (rlptr2 = *regions; rlptr2 != NULL; 
	       rlptr2 = nextLink){
	    nextLink = rlptr2->next;
	    region2 = (MeshRegion *) rlptr2->current;
	    zdo printf("Looking at another region %x, maxNchange %f facets %d\n", region2,
			   region2->maxNormalChange,
			   region2->facetCnt); 
	    /* Merge? */

	    if ((region2 != NULL) &&
		(region2 != region) &&
		(regionIsSmall(region2))
		){
		zdo printf("Merge.\n");

		mergeMeshRegion(region, region2);
		linkedListRemoveNode(regions, (void *) region2);
		linkedListRemoveNode(&region->adjRegions, (void *) region2);

		/* Unlink those neighbor regions that border adjRegion */ 
		unlinkAdjRegions(*regions, region2);

		/* Check sanity */
		linkedListForEachNode(*regions, meshRegionCheck, NULL);
		change = 1;
	      }
	    }
	}
	  /* Recompute adjacent regions of R by getting rid of old adj regions 
	   * R. 
	   */
	  linkedListDestroy(region->adjRegions);
	  region->adjRegions = linkedListNewNode();
	  zdo printf("Recomputing adjacent regions...\n");
	  meshRegionFindAdjRegions(region, NULL);

	  assignRegionQuant((void *) region, NULL);
	}
      }
    } 
}

void
findRegionCurvature(Part *part)
{
  int numVerts;
  Facet *fptr, *refFacet, *maxDeltaFacet;
  int i;
  Vertex *vert;
  float maxDeltaTheta;


  maxDeltaTheta = 0.0;
  maxDeltaFacet = NULL;
  refFacet = NULL;
  for (fptr = part->facets; fptr != NULL; fptr = fptr->next){
    if (fptr->selected != 0){
      float theta, dot;

      if (refFacet == NULL)
	refFacet = fptr;

      dot = pointDot(refFacet->normal, fptr->normal);
      theta = fabs(acos(dot)) * 180.0/M_PI;
      zdo6 printf("theta %f \n", theta);
      if (theta > maxDeltaTheta) {
	maxDeltaFacet = fptr;
	maxDeltaTheta = theta;
      }
    }
  }

}



void
mergeSimilarQuantRegions(LinkedList **regions)
{
  LinkedList *rlptr;
  MeshRegion *region;

  zdo printf("Merging similar quant regions... \n");

  for(rlptr = *regions; rlptr != NULL; rlptr = rlptr->next){
    region = (MeshRegion *) rlptr->current;
    if ((region != NULL) &&
	(region->facetCnt < MinMeshRegionCutoff)){
	int change;
	int iterations = 0;

	zdo printf("Region %x has area %d, less than cutoff, start merging...\n\n",
		   region, region->facetCnt);

	change = 1;
	while ((change)
	      ) {
	  LinkedList *rlptr2;
	  MeshRegion *adjRegion;

	  change = 0;
	  zdo printf("\nIteration %d region facetCnt %d\n",
		     iterations, region->facetCnt);
	  iterations++;
	  
	  /* See if any adjacent regions are mergeable */
	  for (rlptr2 = region->adjRegions; rlptr2 != NULL; 
	       rlptr2 = rlptr2->next){
	    adjRegion = (MeshRegion *) rlptr2->current;
	    if (adjRegion != NULL){
	      float diff;
	      float dot, theta;

	      dot = pointDot(adjRegion->refFacet->normal, 
			     region->refFacet->normal);
	      theta = fabs(acos(dot)) * 180.0/M_PI;


	      diff = fabs(adjRegion->maxNormalChange -
			  region->maxNormalChange);

	      /* Merge? */
	      if (theta < 10.0){

		zdo printf("Adj region %x diff %f, theta %f, merge.\n", adjRegion,
			   diff, theta, adjRegion->facetCnt); 

		mergeMeshRegion(region, adjRegion);
		linkedListRemoveNode(regions, (void *) adjRegion);
		linkedListRemoveNode(&region->adjRegions, (void *) adjRegion);
		change = 1;

		/* Assign quantization to this new region */
		assignRegionQuant((void *) region, NULL);
	      }
	    }
	  }

	  /* Recompute adjacent regions of R by getting rid of old adj regions 
	   * R. 
	   */
	  linkedListDestroy(region->adjRegions);
	  region->adjRegions = linkedListNewNode();
	  zdo printf("Recomputing adjacent regions...\n");
	  meshRegionFindAdjRegions(region, NULL);
	}
      }
  }
}




static int ObjQuantSum;
static int ObjVerts;

int
countQuant(void *node, void *data)
{
  MeshRegion *region = (MeshRegion *) node;
  LinkedList *lptr;
  
  for (lptr = region->flist; lptr != NULL; lptr = lptr->next){
    Facet *facet;

    facet = (Facet *) lptr->current;
    if (facet != NULL){
      if (facet->boundaryFacet)
	ObjQuantSum += facet->numVerts/2 * region->quantLevel +
	  facet->numVerts/2 * 12;
      else
	ObjQuantSum += facet->numVerts * region->quantLevel;
      ObjVerts+= facet->numVerts;
    }
  }
  return ObjQuantSum;
}
 

void
partCountQuant(Part *part, int *partQuantSum, int *partNumVerts)
{

  LinkedList *lptr, *lptr2;
  int sum = 0, numVerts= 1;

  for (lptr2 = part->regions; lptr2 != NULL; lptr2 = lptr2->next){
    MeshRegion *region = (MeshRegion *) lptr2->current;
    GTMesh *gtmesh;

    for (lptr = region->flist; lptr != NULL; lptr = lptr->next){
      Facet *facet;

      facet = (Facet *) lptr->current;
      if (facet != NULL){
	if (facet->boundaryFacet)
	  sum += facet->numVerts/2 * region->quantLevel +
	    facet->numVerts/2 * 12;
	else
	  sum += facet->numVerts * region->quantLevel;
	numVerts += facet->numVerts;
      }
    }

#if 0
    if (gtmesh->triStrips != NULL) {
      sum += gtmesh->triStrips->triCnt;
      numVerts += gtmesh->triStrips->triCnt + 2;
    }
#endif
  }
  *partQuantSum = sum;
  *partNumVerts = numVerts;
  zdo printf("Part %s ave bits %2.2f \n\n", part->name,
	     (float)sum/(float)numVerts);
}
  

void
groupPartOneRegion(Part *part, LinkedList **regions)
{
  MeshRegion *region;
  Facet *fptr;
  float c[4];
  struct timeval timeStamp;
  int i;
  int change = 1;

  zdo6 printf("***VaryQuant is OFF. Setting part to one region ****\n");

  /* Do some initializing. */
  gettimeofday(&timeStamp, NULL);


  region = meshRegionCreate();
  region->refFacet = part->facets;
  region->owner = part;
  
  /*Color this fptr */
  c[0] = randNum(0.0, 1.0);			       
  c[1] = randNum(0.0, 1.0);
  c[2] = randNum(0.0, 1.0);
  c[3] = 1.0;

  while (change){
    change = 0;
    for (fptr = part->facets; fptr != NULL; fptr = fptr->next){
      if ((fptr->timeStamp != timeStamp.tv_sec)){ 
	change = 1;
	/* Inlining meshRegionAddFacet (use fast linked list version) */
	region->facetCnt++;
	linkedListAddNodeFast(&region->flist, (void *) fptr);
	fptr->region = region;

	rgbaCopy(fptr->pickedColor, c);
	fptr->selected = 1.0;
	fptr->timeStamp = timeStamp.tv_sec;
      }
    }
  }
  linkedListAddNode(regions, (void *) region);
}



void
checkFacetsExistence(Part *part)
{
  int found = 0;
  Facet *fptr;

  for (fptr = part->facets; fptr != NULL; fptr = fptr->next){
    LinkedList *rptr;

    found = 0;
    fptr->timeStamp = -1;
    for (rptr = part->regions; rptr != NULL; rptr = rptr->next){
      MeshRegion *region = (MeshRegion *) rptr->current;
      LinkedList *lptr;

      for (lptr = region->flist; lptr != NULL; lptr = lptr->next){
	Facet *testFptr = (Facet *) lptr->current;

	if (fptr == testFptr) {
	  found++;
	}
      }
    }
    if (!found)
      printf("Error!! facet %x not found in regions!!!\n", fptr);
    if (found > 1) 
      printf("Error!! facet %x found mutiple times %d\n", fptr, found);
    zdo6 printf("facet %x found times %d\n", fptr, found);
  }
}



float
triVolume (Vertex *varray, int vindex1, int vindex2, 
	   int vindex3, Point center, float maxBboxLen)
{
  BBox *bbox;
  float volume;
  Point p1, p2, p3;
  float xsize, ysize, zsize;
  float minSize;

  bbox = bboxCreate();

  pointCopy(p1, varray[vindex1].p);
  pointSub(p1, p1, center);
  pointDivScalar(p1, p1, maxBboxLen);

  pointCopy(p2, varray[vindex2].p);
  pointSub(p2, p2, center);
  pointDivScalar(p2, p2, maxBboxLen);

  pointCopy(p3, varray[vindex3].p);
  pointSub(p3, p3, center);
  pointDivScalar(p3, p3, maxBboxLen);

  bboxExtendPoint(bbox, p1);
  bboxExtendPoint(bbox, p2);
  bboxExtendPoint(bbox, p3);
  
  xsize = (bbox->max[0] - bbox->min[0]);
  ysize = (bbox->max[1] - bbox->min[1]);
  zsize = (bbox->max[2] - bbox->min[2]);

  minSize = xsize;
  if (ysize < minSize)
    minSize = ysize;
  if (zsize < minSize)
    minSize = zsize;

  volume = minSize;
  zdo6 printf("tri volume %f \n", volume);
  return volume;
}



void
findMaxQuantLevel(Part *part)
{
    int i;
    TriStrip *s;
    int numStrips = 0, totalLen = 0;
    float aveStripLen = 0;
    Point partBboxCenter;
    float partMaxBboxLen = -1.0;
    int totalTris = 0;
    float triVolumeSum = 0.0;
    float aveTriVolume;
    int q;

    zdo printf("Calculating part %s max quant...\n", part->name);

    for (i = 0; i < 3; i++){
      partBboxCenter[i] = (part->bbox->max[i] + part->bbox->min[i])/ 2.0;
      zdo6 printf("partBboxCenter[%d] = %f\n", i, partBboxCenter[i]);
    }

    for (i = 0; i < 3; i++){
      float len;
      len = (part->bbox->max[i] - part->bbox->min[i]);
      if (len > partMaxBboxLen)
	partMaxBboxLen = len;
    }

    zdo6 printf("partMaxBboxLen %f\n", partMaxBboxLen);

    for (s = part->triStrips->strips; s != NULL; s = s->next){
      for (i = 0; i < s->length - 2 ; i++) {
	if (randNum(0.0, 20.0) < 10.0){
	  int vindex1, vindex2, vindex3;

	  if ((i % 2) == 0){
	    vindex1 = s->tri[i+1];
	    vindex2 = s->tri[i];
	    vindex3 = s->tri[i+2];
	  }
	  else{
	    vindex1 = s->tri[i];
	    vindex2 = s->tri[i+1];
	    vindex3 = s->tri[i+2];
	  }
	  triVolumeSum += triVolume(part->obj->varray, 
				    vindex1, 
				    vindex2, vindex3,
				    partBboxCenter, 
				    partMaxBboxLen);
	  totalTris++;
	}
      }
    }

    aveTriVolume = triVolumeSum / (float) totalTris;
    zdo6 printf("aveTriVolume %f\n", aveTriVolume);

    MaxQuant = 6;

    for (q = 15; q >= 6; q--){
      float cubeLen, cubeVol;
      float pow2 = 1.0;
      float diff;

      for (i = 0; i < q; i++)
	pow2 *= 2.0;

      cubeLen = 1.0 / (float) pow2;
      cubeVol = cubeLen;

      diff = aveTriVolume/cubeVol;
      zdo6 printf("checking quant %d cubeLen %f vol %f diff %f..", q,
		 cubeLen, cubeVol, diff);
      

      if (diff < 5.0){
	MaxQuant = q+1;
	if (MaxQuant > 16)
	  MaxQuant = 16;
	break;
      }
    } 
    part->maxQuant = MaxQuant;
    zdo printf("MaxQuant set to %d\n\n", MaxQuant);
}
    
    

    


    
void
groupObjRegions(ModelObject *obj)
{
  int i;
  LinkedList *regions;

  ObjQuantSum = 0;
  ObjVerts = 0;

  zzdebug = 1;
  for( i = 0; i < obj->partCnt; i++){
    regions = linkedListNewNode();
#if 0
    zdo6 groupPartOneRegion(obj->parts[i], &regions);
    zdo6 linkedListForEachNode(regions, meshRegionPrint, NULL);
    zdo6 findMaxQuantLevel(obj->parts[i]);
    zdo6 groupPartCurvatureRegions(obj->parts[i]);
#endif

#if 0
    groupQuantRegions(obj->parts[i]);
#endif
    groupPartOneRegion(obj->parts[i], &obj->parts[i]->regions);    
    /* Next is for debugging! */
    zdo6 checkFacetsExistence(obj->parts[i]);
    linkedListForEachNode(obj->parts[i]->regions, countQuant, NULL);
  }

  zdo6 printf("Object %s ave quant %f bits\n", obj->name, 
	     (float)ObjQuantSum/(float)ObjVerts);
}



int
countRegions(LinkedList *node)
{
  int cnt =0;
  LinkedList *rlptr2;

  /* See if any adjacent regions are mergeable */
  for (rlptr2 = node; rlptr2 != NULL; 
       rlptr2 = rlptr2->next){
    MeshRegion *region;

    region = (MeshRegion *) rlptr2->current;
    if (region != NULL){
      cnt++;
    }
  }
  return cnt;
}

void
groupPartCurvatureRegions(Part *part)
{
  struct timeval timeStamp;
  Facet *fptr;
  int totalFacets = 0;
  int i;

  /* Do some initializing. */
  gettimeofday(&timeStamp, NULL);

  zdo6 printf("Prior to groupPart old time stamp %ld \n", timeStamp.tv_sec);  
  timeStamp.tv_sec += (int) randNum(1.0, 100.0);
  zdo6 printf("new time stamp %ld \n", timeStamp.tv_sec);  

  if (!VaryQuant){
    groupPartOneRegion(part, &part->regions);
    return;
  }
#if 0
  zdo printf("Calc verts radii...\n");
  partCalculateCurvature(part);
#endif
  
  for (fptr = part->facets; fptr != NULL; fptr = fptr->next)
    totalFacets++;

  MinMeshRegionCutoff = totalFacets / MIN_MESH_REGIONS;
  if (MinMeshRegionCutoff > MAX_MESH_REGION_AREA)
    MinMeshRegionCutoff = MAX_MESH_REGION_AREA;

  if (MinMeshRegionCutoff < MIN_MESH_REGION_AREA)
    MinMeshRegionCutoff = MIN_MESH_REGION_AREA;

  zdo printf("\nVaryQuant part %s...\n", part->name);
  zdo printf("Total facets %d MinMeshRegionCutoff set to %d\n",
	     totalFacets, MinMeshRegionCutoff);

  for (fptr = part->facets; fptr != NULL; fptr = fptr->next){
    if (fptr->timeStamp != timeStamp.tv_sec){
      float c[4];
      MeshRegion *region;

      /*Color this fptr */
      c[0] = randNum(0.0, 1.0);			       
      c[1] = randNum(0.0, 1.0);
      c[2] = randNum(0.0, 1.0);
      c[3] = 1.0;
      
      region = meshRegionCreate();
      region->refFacet = fptr;
      region->owner = part;
      meshRegionAddFacet(region, fptr);

      fptr->pickedColor = rgbaCreate(c);
      fptr->selected = 1.0;
      fptr->timeStamp = timeStamp.tv_sec;
      zdo6 printf("Start merge on fptr %x timeStamp %ld\n",
		 fptr, fptr->timeStamp);

      facetForEachAdjFacet(fptr, mergeAdjFacets, (void *) region);
      linkedListAddNode(&part->regions, (void *) region);
    }
  }

  zdo printf("Finding adjacent regions...\n");
  /* For each of the regions, find adjacent regions of each */
  linkedListForEachNode(part->regions, meshRegionFindAdjRegions, NULL);

  zdo printf("Assigning quantizations to each region.\n"); 

  /* Assign quantization to each regions */
  linkedListForEachNode(part->regions, assignRegionQuant, NULL);

  /* Merge the small regions */
  zdo printf("Regions remaining %d \n", countRegions(part->regions));
  zdo printf("Merging small regions \n");
  zoff(
  mergeSmallRegions(&part->regions);
  )

    zdo printf("Regions remaining %d \n", countRegions(part->regions));
#if 0
  zdo printf("After...\n");
  for(i = 0; i <= part->obj->vertCnt; i++){
      Vertex *vert;

      vert = part->obj->varray+i;
      if ((fabs(vert->p[0]) > 20.0) ||
	  (fabs(vert->p[1]) > 20.0) ||
	  (fabs(vert->p[2]) > 20.0)){
	printf("Bad vert! ");
	printPoint(vert->p);
      }
  }
#endif
 
  zdo6 linkedListForEachNode(part->regions, nukeOldAdjRegions, (void *) part->regions);

  
  linkedListForEachNode(part->regions, meshRegionCheck, NULL);
  zdo6 printf("Sanity checked before merge large small\n");


  zdo printf("Merge large with small...\n");
zoff(
  /* Merge the large with small regions */
      mergeLargeWithSmallRegions(&part->regions);
)


  zdo printf("Regions remaining %d \n", countRegions(part->regions));
  zdo printf("Merge large with large...\n");
 /* Merge the large and large regions */
zoff(
  mergeLargeWithLargeRegions(&part->regions);
)
  zdo printf("Regions remaining %d \n", countRegions(part->regions));
  zdo printf("Merge disjointed regions...\n");
zoff(  
  mergeDisjointedRegions(&part->regions);
)
  zdo printf("Regions remaining %d \n", countRegions(part->regions));
  linkedListForEachNode(part->regions, meshRegionPrint, NULL);
  linkedListForEachNode(part->regions, meshRegionSetRegion, NULL);
  linkedListForEachNode(part->regions, meshRegionCheck, NULL);

#if 0
  if (UserState.testMode != 1)
    partCompressRegions(part);
#endif

}




#define POINT_DIST(p1, p2) (sqrt(((p1[0] - p2[0])*(p1[0] - p2[0])) \
				  + ((p1[1] - p2[1])*(p1[1] - p2[1])) \
				  + ((p1[2] - p2[2])*(p1[2] - p2[2]))))

float 
aveFacetWidth(Facet *fptr, Vertex *varray)
{
  int i;
  float widthSum = 0.0;
  float dist;
  float aveWidth = 0.0;


  for(i = 0; i < fptr->numVerts - 1; i++){
    int vindex1, vindex2;

    vindex1 = fptr->vindex[i];
    vindex2 = fptr->vindex[i+1];
    
    dist = POINT_DIST(varray[vindex1].p, varray[vindex2].p);
    widthSum += dist;
  }

  aveWidth = widthSum / (float) (fptr->numVerts-1);
  return aveWidth;
}

int
mergeAdjFacetsByWidth(Facet *adjFptr, Edge *adjEdge, void *data)
{
  MeshRegion *region = (MeshRegion *) data;
  Facet *fptr = region->refFacet;
  float ratio1, ratio2;
  Vertex *varray = region->owner->obj->varray;

  /* Must not have a region already! */
  if (adjFptr->region != NULL)
    return 1;

  ratio1 = aveFacetWidth(adjFptr, varray)/aveFacetWidth(fptr, varray);
  ratio2 = 1.0/ratio1;
  
  zdo6 printf("adj fptr width %f ratio1 %2.2f ratio2 %2.2f\n",
	     aveFacetWidth(adjFptr, varray), ratio1, ratio2);
  if ((adjFptr->timeStamp != fptr->timeStamp) && 
      ((ratio1 < 1.5) && (ratio2 < 1.5)))  {
    zdo6 printf("adj fptr %x close, cont merge..\n", adjFptr);
    adjFptr->timeStamp = fptr->timeStamp;
    adjFptr->selected = 1;
    adjFptr->pickedColor = rgbaCreate(fptr->pickedColor);
    meshRegionAddFacet(region, adjFptr);
    facetForEachAdjFacet(adjFptr, mergeAdjFacetsByWidth, (void *) region);
  }
  else 
    zdo6 printf("adj not close, cont\n");

  return 1;
}




typedef struct AdjFacetInfo AdjFacetInfo;
struct AdjFacetInfo {
  int depth;
  LinkedList *adjFacets;
  MeshRegion *ownerRegion;
};

int
getAdjFacets(Facet *adjFptr, Edge *adjEdge, void *data)
{
  AdjFacetInfo *adjFacetInfo = (AdjFacetInfo *) data;
  int depth = adjFacetInfo->depth;

  if (adjFptr->region != adjFacetInfo->ownerRegion)
    return 1;

  if (depth > 2) 
    return 1;
  linkedListAddNode(&adjFacetInfo->adjFacets, (void *) adjFptr);
  adjFacetInfo->depth++;
  facetForEachAdjFacet(adjFptr, getAdjFacets, (void *) adjFacetInfo);
  /* Put back the depth */
  adjFacetInfo->depth = depth;
  return 1;
}



void
regionFindMaxNormalChange(MeshRegion *region)
{
  LinkedList *lptr, *prevFlist = NULL;
  int iterations = 0;
  AdjFacetInfo adjFacetInfo;
  float maxNChangeSum = 0.0;
  float aveMaxNChange, maxNChange = -1.0;

  adjFacetInfo.depth = 0;
  adjFacetInfo.adjFacets = NULL;

  zdo6 printf("Find region %x Max Normal change...", region);
 
  for(lptr = region->flist; lptr != NULL; lptr = lptr->next){
    Facet *fptr = (Facet *) lptr->current;
    
    if ((fptr != NULL) && (randNum(0.0,1.0) < 0.5)){
      float normalChange; 
      
      if (iterations > 500) 
	break;

      if ((prevFlist != NULL)){
	if (linkedListFindNode(prevFlist, fptr)){
	  zdo6 printf("facet in prevFlist, cont\n");
	  continue;
	}
      }

      zdo6 printf("iteration %d \n", iterations);

      adjFacetInfo.depth = 0;
      /* Rid of old list */
      if (adjFacetInfo.adjFacets != NULL)
	linkedListDestroy(adjFacetInfo.adjFacets);

      adjFacetInfo.adjFacets = linkedListNewNode();
      adjFacetInfo.ownerRegion = region;
      facetForEachAdjFacet(fptr, getAdjFacets, (void *) &adjFacetInfo);
      normalChange = findMaxNChange(adjFacetInfo.adjFacets, fptr, region->owner->obj);
      
      zdo6 printf("max N change so far %2.2f\n", normalChange);
      maxNChangeSum += normalChange;

      if (normalChange > maxNChange)
	maxNChange = normalChange;

      iterations++;
      prevFlist = adjFacetInfo.adjFacets;
    }
  }

  if (iterations == 0)
    iterations = 1;
  aveMaxNChange = maxNChangeSum  / (float) iterations;
  zdo6 printf("Region's aveMaxNChange %2.2f\n", aveMaxNChange);
  region->maxNormalChange = aveMaxNChange;
}


#define MIN_QUANT_LEVEL 3
int
findQuantFromTriLen(float aveTriLen)
{
  int q;
  int i;
  int finalQuant = MIN_QUANT_LEVEL;
  double triBboxArea;

  triBboxArea = (double) aveTriLen * (double) aveTriLen;

  zdo6 printf("TriBboxArea %lf\n", triBboxArea);
  for (q = 15; q >= MIN_QUANT_LEVEL; q--){
    double cubeLen, cubeVol;
    double pow2 = 1.0;
    double diff;
  
    for (i = 0; i < q; i++)
      pow2 *= 2.0;
    
    cubeLen = 1.0 / (float) pow2;
    cubeVol = cubeLen;

    diff = aveTriLen/cubeVol;
    zdo6 printf("checking quant %d cubeLen %lf cube area %2.12lf diff %8.12lf..", q,
		 cubeLen, cubeVol, diff);

    /* See if aveTriLen is withing K units of grid length (K = 1/error) */
    if (diff < (1.0/VaryQuantError + 0.001)){
      zdo6 printf("Got it q=%d\n", q+1);
      finalQuant = q;
      break;
    }
  }
  return finalQuant;
}



int
assignRegionQuantByTriSize(void *node, void *data)
{
  MeshRegion *region = (MeshRegion *) node;
  Facet *fptr;
  int i;
  Point bboxCenter;
  float maxBboxLen = -1.0;
  int totalTris = 0;
  float triBboxLenSum = 0.0;
  float aveTriBboxLen;
  LinkedList *lptr;
  Vertex *varray = region->owner->obj->varray;

  zdo6 printf("assigning region %x quant...\n");
  meshRegionFindBbox(region);
  for (i = 0; i < 3; i++){
      bboxCenter[i] = (region->bbox->max[i] + region->bbox->min[i])/ 2.0;
  }

  zdo6 printf("bboxCenter: ");
  zdo6 printPoint(bboxCenter);

  for (i = 0; i < 3; i++){
    float len;
    len = region->bboxSize[i];
    if (len > maxBboxLen)
      maxBboxLen = len;
  }
  zdo6 printf("maxBboxLen %f\n", maxBboxLen);

  for (lptr = region->flist; lptr != NULL; lptr = lptr->next){
    fptr = (Facet *) lptr->current;

    if (fptr != NULL){
      if (fptr->numVerts == 3) {
	triBboxLenSum += triVolume(varray,
				   fptr->vindex[0], 
				   fptr->vindex[1], 
				   fptr->vindex[2], 
				   bboxCenter,
				   maxBboxLen);
	totalTris++;
      }
      else 
	if (fptr->numVerts == 4){
	triBboxLenSum += triVolume(varray,
				   fptr->vindex[0], 
				   fptr->vindex[1], 
				   fptr->vindex[2], 
				   bboxCenter,
				   maxBboxLen);
	totalTris++;
	triBboxLenSum += triVolume(varray,
				   fptr->vindex[2], 
				   fptr->vindex[3], 
				   fptr->vindex[0], 
				   bboxCenter,
				   maxBboxLen);
	totalTris++;
	}
    }
  }

  aveTriBboxLen = triBboxLenSum / (float) totalTris;
  zdo6 printf("aveTriBboxLen %f\n", aveTriBboxLen);
  region->quantLevel = findQuantFromTriLen(aveTriBboxLen);
  region->aveTriBboxLen = aveTriBboxLen;
  regionFindMaxNormalChange(region);

  /* Snap to user specified levels */
  if (region->quantLevel < VaryQuantLow)
    region->quantLevel = VaryQuantLow;
  if (region->quantLevel > VaryQuantHigh)
    region->quantLevel = VaryQuantHigh;

  zdo6 printf("got quant %d \n", region->quantLevel);
  
  /* (region->maxNormalChange > 90.0) */
  if (region->maxNormalChange > 90.0) {
    region->quantLevel += 1;
    zdo6 printf("up-ed by normal change, quant %d \n", region->quantLevel);
  }
  return region->quantLevel;
}
 

int
isRegionMergeableByTriSize(MeshRegion *r1, MeshRegion *r2, int quantThreshold)
{
  float ratio1, ratio2;
  float diff1,diff2;
  LinkedList *lptr;
  MeshRegion *rcombined;
  Facet *fptr;

  ratio1 = r1->aveTriBboxLen / r2->aveTriBboxLen;
  ratio2 = 1.0/ratio1;

  zdo printf("mergeable? r1len %2.2f ratio1 %2.2f r2len %2.2f ratio2 %2.2f\n",
	     r1->aveTriBboxLen,  ratio1, r2->aveTriBboxLen, ratio2);
  if (!((ratio1 < 3.0) && (ratio2 < 3.0)))
    return 0;

  rcombined = meshRegionCreate();
  rcombined->refFacet = r1->refFacet;
  rcombined->owner = r1->owner;
  for (lptr = r1->flist; lptr != NULL; lptr = lptr->next){
    fptr = (Facet *) lptr->current;
    meshRegionAddFacet(rcombined, fptr);
    fptr->region = r1;
  }
  for (lptr = r2->flist; lptr != NULL; lptr = lptr->next){
     fptr = (Facet *) lptr->current;
     meshRegionAddFacet(rcombined, (Facet *) lptr->current);
     fptr->region = r2;
  }

  assignRegionQuantByTriSize(rcombined, NULL);
  diff1 = fabs(r1->quantLevel - rcombined->quantLevel);
  diff2 = fabs(r2->quantLevel - rcombined->quantLevel);
 
  zdo printf("quant diff %2.2f %2.2f \n", diff1, diff2);

  meshRegionDestroy(rcombined);
  if ((diff1 <= quantThreshold) && (diff2 <= quantThreshold))
    return 1;
  
  return 0;
}

void
mergeSmallRegionsByTriSize(LinkedList **regions)
{
  LinkedList *rlptr;
  MeshRegion *region;
  
  zdo printf("Merging Small regions... \n");

  for(rlptr = *regions; rlptr != NULL; rlptr = rlptr->next){
    region = (MeshRegion *) rlptr->current;
    if (region != NULL){
      if ((region->facetCnt < MinMeshRegionCutoff)) {
	int change;
	int iterations = 0;

	zdo printf("Region %x has area %d, more than cutoff, start merging...\n\n",
		   region, region->facetCnt);

	change = 1;

	while ((change)
	       ) {
	  LinkedList *rlptr2, *nextLink;
	  MeshRegion *adjRegion;

	  zdo6 printAdjRegions(region);
	  change = 0;
	  zdo printf("\nIteration %d maxNormalChange %f region facetCnt %d\n",
		     iterations, region->maxNormalChange, region->facetCnt);
	  iterations++;
	  
	  /* See if any adjacent regions are mergeable */
	  for (rlptr2 = region->adjRegions; rlptr2 != NULL; 
	       rlptr2 = nextLink){
	    nextLink = rlptr2->next;
	    adjRegion = (MeshRegion *) rlptr2->current;
	    zdo 
	      if (adjRegion != NULL)
		printf("Looking at Adj region %x facetCnt %d \n", adjRegion,
		       adjRegion->facetCnt);
	    if ((adjRegion != NULL) &&
		(adjRegion->facetCnt < MinMeshRegionCutoff)){
	      zdo printf("Adj region is mergeable. \n");
	      /* Check sanity */
	      zdo linkedListForEachNode(*regions, meshRegionCheck, NULL);
	      zdo printf("Sanity checked before merging.\n");
	      mergeMeshRegion(region, adjRegion);
	      linkedListRemoveNode(regions, (void *) adjRegion);

	      /* Check sanity */
	      linkedListForEachNode(*regions, meshRegionCheck, NULL);
	      
	      linkedListRemoveNode(&region->adjRegions, (void *) adjRegion);
	      /* Unlink those neighbor regions that border adjRegion */ 
	      unlinkAdjRegions(*regions, adjRegion);
	      change = 1;
	      }
	    }
	  /* Recompute adjacent regions of R by getting rid of old adj regions 
	   * R. 
	   */
	  linkedListDestroy(region->adjRegions);
	  region->adjRegions = linkedListNewNode();
	  zdo printf("Recomputing adjacent regions...\n");
	  meshRegionFindAdjRegions(region, NULL);
	  assignRegionQuantByTriSize((void *) region, NULL);
	}
      }
    }
  } 
}


void
mergeLargeWithSmallRegionsByTriSize(LinkedList **regions)
{
  LinkedList *rlptr;
  MeshRegion *region;
  
  zdo printf("Merging By Width Large & Small regions... \n");

  for(rlptr = *regions; rlptr != NULL; rlptr = rlptr->next){
    region = (MeshRegion *) rlptr->current;
    if (region != NULL){
      if ((region->facetCnt > MinMeshRegionCutoff)) {
	int change;
	int iterations = 0;

	zdo printf("Region %x has area %d, more than cutoff, start merging...\n\n",
		   region, region->facetCnt);

	change = 1;

	while ((change)
	       ) {
	  LinkedList *rlptr2, *nextLink;
	  MeshRegion *adjRegion;

	  zdo6 printAdjRegions(region);
	  change = 0;
	  zdo printf("\nIteration %d maxNormalChange %f region facetCnt %d\n",
		     iterations, region->maxNormalChange, region->facetCnt);
	  iterations++;
	  
	  /* See if any adjacent regions are mergeable */
	  for (rlptr2 = region->adjRegions; rlptr2 != NULL; 
	       rlptr2 = nextLink){
	    nextLink = rlptr2->next;
	    adjRegion = (MeshRegion *) rlptr2->current;
	    zdo 
	      if (adjRegion != NULL)
		printf("Looking at Adj region %x facetCnt %d \n", adjRegion,
		       adjRegion->facetCnt);

	    if ((adjRegion != NULL)
		&& (((adjRegion->facetCnt < MinMeshRegionCutoff)
		    && (isRegionMergeableByTriSize(region, adjRegion, 2)))
#if 0
		    || (adjRegion->quantLevel < region->quantLevel) 
#endif
		    || (adjRegion->facetCnt < VerySmallRegionCutoff))){
	      zdo printf("Adj region is mergeable. \n");
	      /* Check sanity */
	      zdo linkedListForEachNode(*regions, meshRegionCheck, NULL);
	      zdo printf("Sanity checked before merging.\n");
	      mergeMeshRegion(region, adjRegion);
	      linkedListRemoveNode(regions, (void *) adjRegion);

	      /* Check sanity */
	      linkedListForEachNode(*regions, meshRegionCheck, NULL);
	      
	      linkedListRemoveNode(&region->adjRegions, (void *) adjRegion);
	      /* Unlink those neighbor regions that border adjRegion */ 
	      unlinkAdjRegions(*regions, adjRegion);
	      change = 1;
	      }
	    }

	  /* Recompute adjacent regions of R by getting rid of old adj regions 
	   * R. 
	   */
	  linkedListDestroy(region->adjRegions);
	  region->adjRegions = linkedListNewNode();
	  zdo printf("Recomputing adjacent regions...\n");
	  meshRegionFindAdjRegions(region, NULL);
	  assignRegionQuantByTriSize((void *) region, NULL);
	}
      }
    }
  } 
}



void
mergeLargeWithLargeRegionsByTriSize(LinkedList **regions)
{
  LinkedList *rlptr;
  MeshRegion *region;
  
  zdo printf("Merging By Width Large & Small regions... \n");

  for(rlptr = *regions; rlptr != NULL; rlptr = rlptr->next){
    region = (MeshRegion *) rlptr->current;
    if (region != NULL){
      if ((region->facetCnt > MinMeshRegionCutoff)) {
	int change;
	int iterations = 0;

	zdo printf("Region %x has area %d, more than cutoff, start merging...\n\n",
		   region, region->facetCnt);

	change = 1;

	while ((change)
	       ) {
	  LinkedList *rlptr2, *nextLink;
	  MeshRegion *adjRegion;

	  zdo printAdjRegions(region);
	  change = 0;
	  zdo printf("\nIteration %d maxNormalChange %f region facetCnt %d\n",
		     iterations, region->maxNormalChange, region->facetCnt);
	  iterations++;
	  
	  /* See if any adjacent regions are mergeable */
	  for (rlptr2 = region->adjRegions; rlptr2 != NULL; 
	       rlptr2 = nextLink){
	    nextLink = rlptr2->next;
	    adjRegion = (MeshRegion *) rlptr2->current;
	    zdo 
	      if (adjRegion != NULL)
		printf("Looking at Adj region %x facetCnt %d \n", adjRegion,
		       adjRegion->facetCnt);

	    if ((adjRegion != NULL)
		&& (((adjRegion->facetCnt > MinMeshRegionCutoff)
		    && (isRegionMergeableByTriSize(region, adjRegion, 1)))
		    )){
	      zdo printf("Adj region is mergeable. \n");
	      /* Check sanity */
	      zdo linkedListForEachNode(*regions, meshRegionCheck, NULL);
	      zdo printf("Sanity checked before merging.\n");
	      mergeMeshRegion(region, adjRegion);
	      linkedListRemoveNode(regions, (void *) adjRegion);

	      /* Check sanity */
	      linkedListForEachNode(*regions, meshRegionCheck, NULL);
	      
	      linkedListRemoveNode(&region->adjRegions, (void *) adjRegion);
	      /* Unlink those neighbor regions that border adjRegion */ 
	      unlinkAdjRegions(*regions, adjRegion);
	      change = 1;
	    }
	  }


	  /* Recompute adjacent regions of R by getting rid of old adj regions 
	   * R. 
	   */
	  linkedListDestroy(region->adjRegions);
	  region->adjRegions = linkedListNewNode();
	  zdo printf("Recomputing adjacent regions...\n");
	  meshRegionFindAdjRegions(region, NULL);
	  assignRegionQuantByTriSize((void *) region, NULL);
	}
      }
    }
  } 
}

void
mergeDisjointedRegionsByTriSize(LinkedList **regions)
{
  LinkedList *rlptr;
  MeshRegion *region;


  zdo printf("Merging disjointed regions... \n");

  for(rlptr = *regions; rlptr != NULL; rlptr = rlptr->next){
    region = (MeshRegion *) rlptr->current;
    if (region != NULL){
      if (
	  (regionIsSmall(region))
	  ) {
	int change;
	int iterations = 0;

	zdo printf("Region %x has area %d, more than cutoff, start merging...\n\n",
		   region, region->facetCnt);

	change = 1;


	while ((change)
	       ) {
	  LinkedList *rlptr2, *nextLink;
	  MeshRegion *region2;
	  
	  change = 0;
	  zdo printf("\nIteration %d maxNormalChange %f region facetCnt %d\n",
		     iterations, region->maxNormalChange, region->facetCnt);
	  iterations++;
	  /* See if any adjacent regions are mergeable */
	  for (rlptr2 = *regions; rlptr2 != NULL; 
	       rlptr2 = nextLink){
	    nextLink = rlptr2->next;
	    region2 = (MeshRegion *) rlptr2->current;
	    zdo printf("Looking at another region %x, maxNchange %f facets %d\n", region2,
			   region2->maxNormalChange,
			   region2->facetCnt); 
	    /* Merge? */

	    if ((region2 != NULL) &&
		(region2 != region) &&
		(regionIsSmall(region2))
		){
		zdo printf("Merge.\n");

		mergeMeshRegion(region, region2);
		linkedListRemoveNode(regions, (void *) region2);
		linkedListRemoveNode(&region->adjRegions, (void *) region2);

		/* Unlink those neighbor regions that border adjRegion */ 
		unlinkAdjRegions(*regions, region2);

		/* Check sanity */
		linkedListForEachNode(*regions, meshRegionCheck, NULL);
		change = 1;
	      }
	    }
	  }

	  /* Recompute adjacent regions of R by getting rid of old adj regions 
	   * R. 
	   */
	  linkedListDestroy(region->adjRegions);
	  region->adjRegions = linkedListNewNode();
	  zdo printf("Recomputing adjacent regions...\n");
	  meshRegionFindAdjRegions(region, NULL);

	  assignRegionQuantByTriSize((void *) region, NULL);
	}
      }
    } 
}




