#include <math.h>
#include "GR_assert.h"
#include "GR_events.h"
#include "GR_misc.h"
#include "GR_Cell.h"
#include "GR_Geometry.h"
#include "GR_Face.h"
#include "GR_Mesh.h"
#include "GR_Quality.h"
#include "GR_Vec.h"
#include "GR_Vertex.h"

#ifdef IRIX6
#include <values.h>
#endif

//@ Cell quality functions

//@@ Maximum dihedral angle
void vMaxDihed(const CellSkel* const pC, int * const piN, double adRes[])
     // Range: 0 - 180.  (Actually 70.5 - 180)
{
  assert(piN != NULL);
  assert(adRes != NULL);
  *piN = 1;
  adRes[0] = pC->dMaxDihed();
}

//@@ Minimum dihedral angle
void vMinDihed(const CellSkel* const pC, int * const piN, double adRes[])
     // Range: 0 - 180.  (Actually 0 - 70.5)
{
  assert(piN != NULL);
  assert(adRes != NULL);
  *piN = 1;
  adRes[0] = pC->dMinDihed();
}

//@@ All dihedral angles
void vAllDihed(const CellSkel* const pC, int * const piN, double adRes[])
     // Range: 0 - 180.  (Equilateral tet is 70.5)
{
  assert(piN != NULL);
  assert(adRes != NULL);
  pC->vAllDihed(adRes, piN);
}

//@@ Maximum solid angle
void vMaxSolid(const CellSkel* const pC, int * const piN, double adRes[])
     // Range: 0 - 360.  (Smallest max is actually around 31.6.)
{
  assert(piN != NULL);
  assert(adRes != NULL);
  *piN = 1;
  adRes[0] = pC->dMaxSolid();
}

//@@ Minimum solid angle
void vMinSolid(const CellSkel* const pC, int * const piN, double adRes[])
     // Range: 0 - 360.  (Largest min is actually around 31.6.)
{
  assert(piN != NULL);
  assert(adRes != NULL);
  *piN = 1;
  adRes[0] = pC->dMinSolid();
}

//@@ All solid angles
void vAllSolid(const CellSkel* const pC, int * const piN, double adRes[])
     // Range: 0 - 360.  (Equilateral tet is 31.6)
{
  assert(piN != NULL);
  assert(adRes != NULL);
  pC->vAllSolid(adRes, piN);
}

//@@ Ratio of incircle to circumcircle radius
void vInToCircumCircle(const CellSkel* const pC, int * const piN,
		       double adRes[])
     // Range:  Normalized to be 0 to 1
{
  assert(pC->eType() == Cell::eTriCell ||
	 pC->eType() == Cell::eTriCV);
  assert(piN != NULL);
  assert(adRes != NULL);
  assert(pC->qValid());
  const double *adLoc0 = pC->pVVert(0)->adCoords();
  const double *adLoc1 = pC->pVVert(1)->adCoords();
  const double *adLoc2 = pC->pVVert(2)->adCoords();
  double dArea;
  vNormal2D(adLoc0, adLoc1, adLoc2, &dArea);
  dArea /= 2;
  double dLenA = dDIST2D(adLoc0, adLoc1);
  double dLenB = dDIST2D(adLoc1, adLoc2);
  double dLenC = dDIST2D(adLoc2, adLoc0);
  *piN = 1;
  adRes[0] = ( 16*dArea*dArea / ((dLenA + dLenB + dLenC) *
				 (dLenA * dLenB * dLenC)) );
}

//@@ Ratio of insphere to circumsphere radius
void vInToCircumSphere(const CellSkel* const pCIn, int * const piN,
		       double adRes[])
     // Range:  Normalized to be 0 to 1
{
  assert(pCIn->eType() == Cell::eTet);
  assert(piN != NULL);
  assert(adRes != NULL);
  const TetCell * const pC = dynamic_cast<const TetCell*>(pCIn);
  assert(pC->qValid());
  double dVol = pC->dSize();
  double dAreaA = pC->pFFace(0)->dSize();
  double dAreaB = pC->pFFace(1)->dSize();
  double dAreaC = pC->pFFace(2)->dSize();
  double dAreaD = pC->pFFace(3)->dSize();
  const Vert *pVA = pC->pVVert(0);
  const Vert *pVB = pC->pVVert(1);
  const Vert *pVC = pC->pVVert(2);
  const Vert *pVD = pC->pVVert(3);
  double dLenAB = dDIST3D(pVA->adCoords(), pVB->adCoords());
  double dLenBC = dDIST3D(pVB->adCoords(), pVC->adCoords());
  double dLenCA = dDIST3D(pVC->adCoords(), pVA->adCoords());
  double dLenAD = dDIST3D(pVA->adCoords(), pVD->adCoords());
  double dLenBD = dDIST3D(pVB->adCoords(), pVD->adCoords());
  double dLenCD = dDIST3D(pVC->adCoords(), pVD->adCoords());
  double dProd1 = dLenAB*dLenCD;
  double dProd2 = dLenBC*dLenAD;
  double dProd3 = dLenCA*dLenBD;
  *piN = 1;
  adRes[0] = ( 216*dVol*dVol /
	   ( (dAreaA + dAreaB + dAreaC + dAreaD) *
	     sqrt((dProd1 + dProd2 + dProd3) * (dProd1 + dProd2 - dProd3) * 
		  (dProd2 + dProd3 - dProd1) * (dProd3 + dProd1 - dProd2))
	     + 1.e-100 )
	   );
}

//@@ Normalized ratio of shortest edge length to circumradius
void vShortestToRadius(const CellSkel* const pC, int * const piN,
		       double adRes[])
     // Range:  Normalized to be 0 to 1
{
  assert(pC->eType() == Cell::eTet || pC->eType() == Cell::eTetCV ||
	 pC->eType() == Cell::eTriCell || pC->eType() == Cell::eTriCV);
  assert(piN != NULL);
  assert(adRes != NULL);
  *piN = 1;
  if (pC->eType() == Cell::eTriCell ||
      pC->eType() == Cell::eTriCV) {
    double dLen01 = dDIST2D(pC->pVVert(0)->adCoords(),
			    pC->pVVert(1)->adCoords());
    double dLen12 = dDIST2D(pC->pVVert(1)->adCoords(),
			    pC->pVVert(2)->adCoords());
    double dLen20 = dDIST2D(pC->pVVert(2)->adCoords(),
			    pC->pVVert(0)->adCoords());
    double dShortest = min(min(dLen01, dLen12), dLen20);
    // For an equilateral triangle with edge length 1, the circumradius
    // is sqrt(1/3), so we need this factor to normalize.
    double dRadius = dynamic_cast<const TriCell*>(pC)->dCircumradius();
    adRes[0] = sqrt(1./3.) * dShortest / dRadius;
  }
  else {
//      double dLen01_Sq = dDIST_SQ_3D(pC->pVVert(0)->adCoords(),
//  				   pC->pVVert(1)->adCoords());
//      double dLen02_Sq = dDIST_SQ_3D(pC->pVVert(0)->adCoords(),
//  				   pC->pVVert(2)->adCoords());
//      double dLen03_Sq = dDIST_SQ_3D(pC->pVVert(0)->adCoords(),
//  				   pC->pVVert(3)->adCoords());
//      double dLen12_Sq = dDIST_SQ_3D(pC->pVVert(1)->adCoords(),
//  				   pC->pVVert(2)->adCoords());
//      double dLen13_Sq = dDIST_SQ_3D(pC->pVVert(1)->adCoords(),
//  				   pC->pVVert(3)->adCoords());
//      double dLen23_Sq = dDIST_SQ_3D(pC->pVVert(2)->adCoords(),
//  				   pC->pVVert(3)->adCoords());
    // For an equilateral tetrahedron with edge length 1, the circumradius
    // is sqrt(3/8), so we need this factor to normalize.
    assert(pC->qValid());
    const double * const adA = pC->pVVert(0)->adCoords();
    const double * const adB = pC->pVVert(1)->adCoords();
    const double * const adC = pC->pVVert(2)->adCoords();
    const double * const adD = pC->pVVert(3)->adCoords();
      
    double adCircCent[3];
    {
      double adRow1[] = {adA[0] - adB[0],
			 adA[1] - adB[1],
			 adA[2] - adB[2]};
      double adRow2[] = {adA[0] - adC[0],
			 adA[1] - adC[1],
			 adA[2] - adC[2]};
      double adRow3[] = {adA[0] - adD[0],
			 adA[1] - adD[1],
			 adA[2] - adD[2]};
      
      double dRHS1 = 0.5 * (dDOT3D(adA, adA) - dDOT3D(adB, adB));
      double dRHS2 = 0.5 * (dDOT3D(adA, adA) - dDOT3D(adC, adC));
      double dRHS3 = 0.5 * (dDOT3D(adA, adA) - dDOT3D(adD, adD));
      
      vSolve3By3(adRow1, adRow2, adRow3, dRHS1, dRHS2, dRHS3, adCircCent);
    }

    double dLenSq_AB = dDIST_SQ_3D(adA, adB);
    double dLenSq_AC = dDIST_SQ_3D(adA, adC);
    double dLenSq_AD = dDIST_SQ_3D(adA, adD);
    double dLenSq_BC = dDIST_SQ_3D(adB, adC);
    double dLenSq_BD = dDIST_SQ_3D(adB, adD);
    double dLenSq_CD = dDIST_SQ_3D(adC, adD);

    double dShortestSq = min(min(min(dLenSq_AB, dLenSq_AC),
				 min(dLenSq_AD, dLenSq_BC)),
			     min(dLenSq_BD, dLenSq_CD));
    double dShortestLen = sqrt(dShortestSq);
    double dRadius = dDIST3D(adCircCent, adA);
    
    adRes[0] = sqrt(0.375) * dShortestLen / dRadius;
  }
}

//@@ Marcum's sliver measure Q1 (AIAA 95-0212)
void vSliverQ1(const CellSkel* const pCIn, int * const piN, double adRes[])
     // Range:  Normalized to be 0 to 1
{
  assert(pCIn->eType() == CellSkel::eTet || pCIn->eType() == CellSkel::eTetCV);
  assert(piN != NULL);
  assert(adRes != NULL);
//   const TetCell * const pC = dynamic_cast<const TetCell*>(pCIn);
//   assert(pC->qValid());
  const CellSkel* pC = pCIn;
  double dVol = pC->dSize();
  const Vert *pVA = pC->pVVert(0);
  const Vert *pVB = pC->pVVert(1);
  const Vert *pVC = pC->pVVert(2);
  const Vert *pVD = pC->pVVert(3);
  double dLenAB = dDIST3D(pVA->adCoords(), pVB->adCoords());
  double dLenBC = dDIST3D(pVB->adCoords(), pVC->adCoords());
  double dLenCA = dDIST3D(pVC->adCoords(), pVA->adCoords());
  double dLenAD = dDIST3D(pVA->adCoords(), pVD->adCoords());
  double dLenBD = dDIST3D(pVB->adCoords(), pVD->adCoords());
  double dLenCD = dDIST3D(pVC->adCoords(), pVD->adCoords());
  // Shewchuk uses the RMS length
  double dAveLen = sqrt((dLenAB*dLenAB + dLenBC*dLenBC +
			 dLenCA*dLenCA + dLenAD*dLenAD +
			 dLenBD*dLenBD + dLenCD*dLenCD) / 6.);
  // Marcum used the arithmetic average length
//   double dAveLen = (dLenAB + dLenBC + dLenCA + dLenAD + dLenBD +
// 		  dLenCD) / 6.;
  *piN = 1;
  adRes[0] = 6 * M_SQRT2 * dVol / (pow(dAveLen, 3.));
}

//@@ Marcum's skewness measure Q2 (AIAA 95-0212)
void vSkewnessQ2(const CellSkel* const pCIn, int * const piN, double adRes[])
     // Range:  Normalized to be 0 to 1
{
  assert(pCIn->eType() == Cell::eTet);
  assert(piN != NULL);
  assert(adRes != NULL);
  const TetCell * const pC = dynamic_cast<const TetCell*>(pCIn);
  assert(pC->qValid());
  const double Rsqrt3 = sqrt(3.);
  double dVol = pC->dSize();
  double dCircRad = (dynamic_cast<const SimplexCell*>(pC))->dCircumradius();
  *piN = 1;
  adRes[0] = 1.125 * Rsqrt3 * dVol / (pow(dCircRad, 3.));
}

//@@ Ratio of area to perimeter squared
void vAreaPerim(const CellSkel* const pC, int * const piN, double adRes[])
     // Range:  Normalized to be 0 to 1
{
  assert(piN != NULL);
  assert(adRes != NULL);
  assert(pC->eType() == CellSkel::eTriCell ||
	 pC->eType() == CellSkel::eTriCV);
  assert(pC->qValid());
  const double RConst = 12. * sqrt(3.);  // about 20
  const double *adLoc0 = pC->pVVert(0)->adCoords();
  const double *adLoc1 = pC->pVVert(1)->adCoords();
  const double *adLoc2 = pC->pVVert(2)->adCoords();

  double dArea, dLenA, dLenB, dLenC;
  vNormal2D(adLoc0, adLoc1, adLoc2, &dArea);
  dArea /= 2;
  dLenA = dDIST2D(adLoc0, adLoc1);
  dLenB = dDIST2D(adLoc1, adLoc2);
  dLenC = dDIST2D(adLoc2, adLoc0);
  double dPerim = dLenA + dLenB + dLenC;
  *piN = 1;
  adRes[0] = RConst * dArea / (dPerim * dPerim);
}

//@@ Ratio of volume to total surface area^1.5
void vVolArea(const CellSkel* const pCIn, int * const piN, double adRes[])
     // Range:  Normalized to be 0 to 1
{
  assert(pCIn->eType() == Cell::eTet);
  assert(piN != NULL);
  assert(adRes != NULL);
  const TetCell * const pC = dynamic_cast<const TetCell*>(pCIn);
  assert(pC->qValid());
  const double RConst = 2 * M_SQRT2 * pow(3.,1.75);  // about 20
  double dVol = pC->dSize();
  double dAreaA = pC->pFFace(0)->dSize();
  double dAreaB = pC->pFFace(1)->dSize();
  double dAreaC = pC->pFFace(2)->dSize();
  double dAreaD = pC->pFFace(3)->dSize();
  double dTotArea = dAreaA + dAreaB + dAreaC + dAreaD;
  *piN = 1;
  adRes[0] = RConst * dVol / pow(dTotArea, 1.5);
}

static double dQKolSmirn(const double dLam)
  // Computes the Kolmogorov-Smirnov significance, Q_KS
  // 1 - Q_KS is the confidence with which we can reject the null
  // hypothesis that the distributions are the same.
{
  // If the argument lambda is small, then the two distributions are
  // statistically the same.
  if (dLam < 0.1) return 1;
  // If the argument is large, the two are unquestionably different.
  else if (dLam > 4) return 0;
  
  double dSum = 0, d2LamSq = dLam*dLam*2.;
  // Compute enough terms that the last will be smaller than 1.e-10
  int iMaxTerms = 1 + int(ceil(sqrt(25. / d2LamSq)));
  for (int i = 1; i <= iMaxTerms; i += 2) 
    dSum += exp(- i*i * d2LamSq) - exp(- (i+1)*(i+1) * d2LamSq);
  return (2*dSum);
}
  
Quality::Quality(Mesh* const pMIn, const int iMeas)
  : pM(pMIn), pvfQualityFunc(NULL), iQualEvals(0), iMaxEvals(20), iNBins(0), 
    dMaxQual(0), dMinQual(0), aiNumQualVals(NULL), adSignifImp(NULL),
    a2dBinFrac(NULL), adMinQual(NULL), adMaxQual(NULL), qMaxMin(true)
{
  static const int aiNBins[] = {30, 30, 30, 30, 16, 30, 25, 25, 25, 25, 25};
  static const double adMax[] = {180, 180, 180, 360, 32, 360, 1, 1, 1, 1, 1};
  static const double adMin[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
  static const bool aqMaxMin[] = {false, true, true, false, true, true,
				  true, true, true, true, true};
  assert (iMeas >= 0 && iMeas < 11);
  iNBins = aiNBins[iMeas];
  dMaxQual = adMax[iMeas];
  dMinQual = adMin[iMeas];
  switch (iMeas) {
  case 0: pvfQualityFunc = vMaxDihed; break;
  case 1: pvfQualityFunc = vMinDihed; break;
  case 2: pvfQualityFunc = vAllDihed; break;
  case 3: pvfQualityFunc = vMaxSolid; break;
  case 4: pvfQualityFunc = vMinSolid; break;
  case 5: pvfQualityFunc = vAllSolid; break;
  case 6:
    if (!pM->qSimplicial()) {
      vMessage(0, "%s\n%s\n%s\n",
	       "Can only use the requested quality criterion, ratio of inball",
	       "to circumball radius, for simplicial meshes.",
	       "Setting quality measure to all dihedral angles.");
      pvfQualityFunc = vAllDihed;
    }
    else if (pM->eType() == Mesh::eMesh2D)
      pvfQualityFunc = vInToCircumCircle;
    else
      pvfQualityFunc = vInToCircumSphere;
    break;
  case 7: 
    if ((!pM->qSimplicial()) || (pM->eType() != Mesh::eVolMesh)) {
      vMessage(0, "%s\n%s\n%s\n",
	       "Can only use the requested quality criterion, Marcum's sliver",
	       "measure, for tetrahedral meshes.",
	       "Setting quality measure to all dihedral angles.");
      pvfQualityFunc = vAllDihed;
    }
    else
      pvfQualityFunc = vSliverQ1;
    break;
  case 8: 
    if ((!pM->qSimplicial()) || (pM->eType() != Mesh::eVolMesh)) {
      vMessage(0, "%s\n%s\n%s\n",
	       "Can only use the requested quality criterion, Marcum's skewness",
	       "measure, for tetrahedral meshes.",
	       "Setting quality measure to all dihedral angles.");
      pvfQualityFunc = vAllDihed;
    }
    else
      pvfQualityFunc = vSkewnessQ2;
    break;
  case 9:
    if (pM->eType() == Mesh::eMesh2D)
      pvfQualityFunc = vAreaPerim;
    else
      pvfQualityFunc = vVolArea;
    break;
  case 10:
    pvfQualityFunc = vShortestToRadius;
    break;
  default:
    vFatalError("Invalid quality measure selected", "Quality::Quality()");
  }

  a2dBinFrac = new double*[iMaxEvals];
  adMinQual = new double[iMaxEvals];
  adMaxQual = new double[iMaxEvals];
  adSignifImp = new double[iMaxEvals];
  adSignifImp[0] = 1;
  aiNumQualVals = new int[iMaxEvals];
  qMaxMin = aqMaxMin[iMeas];
  for (int i = 0; i < iMaxEvals; i++) {
    adMinQual[i] = adMax[iMeas]; // This pre-sets the min/max measured 
    adMaxQual[i] = adMin[iMeas]; // quality appropriately for use.
    a2dBinFrac[i] = new double[iNBins];
    for (int ii = 0; ii < iNBins; ii++)
      a2dBinFrac[i][ii] = 0.;
  }
}

Quality::Quality(const Quality& Q) :
  pM(Q.pM), pvfQualityFunc(Q.pvfQualityFunc), iQualEvals(Q.iQualEvals),
  iMaxEvals(Q.iMaxEvals), iNBins(Q.iNBins), dMaxQual(Q.dMaxQual),
  dMinQual(Q.dMinQual), aiNumQualVals(), adSignifImp(), a2dBinFrac(),
  adMinQual(), adMaxQual(), qMaxMin(Q.qMaxMin)
{
  a2dBinFrac = new double*[iMaxEvals];
  adMinQual = new double[iMaxEvals];
  adMaxQual = new double[iMaxEvals];
  adSignifImp = new double[iMaxEvals];
  aiNumQualVals = new int[iMaxEvals];

  for (int i = 0; i < iMaxEvals; i++) {
    adMinQual[i] = Q.adMinQual[i]; // This pre-sets the min/max measured 
    adMaxQual[i] = Q.adMaxQual[i]; // quality appropriately for use.
    adSignifImp[i] = Q.adSignifImp[i];
    aiNumQualVals[i] = Q.aiNumQualVals[i];
    a2dBinFrac[i] = new double[iNBins];
    for (int ii = 0; ii < iNBins; ii++)
      a2dBinFrac[i][ii] = Q.a2dBinFrac[i][ii];
  }
}
  

// Quality::Quality(Mesh* const pMIn,
//                  void vQual(const Cell* const, int * const , double []),
//                  const int iNB, const double dMin,
//                  const double dMax, const int iMaxEv)
// {
//   pM = pMIn;
//   pvfQualityFunc = vQual;
//   iQualEvals = 0;
//   iMaxEvals = iMaxEv;
//   iNBins = iNB;
//   dMaxQual = dMax;
//   dMinQual = dMin;
//   a2dBinFrac = new double*[iMaxEvals];
//   adMinQual = new double[iMaxEvals];
//   adMaxQual = new double[iMaxEvals];
//   adSignifImp = new double[iMaxEvals];
//   adSignifImp[0] = 1;
//   aiNumQualVals = new int[iMaxEvals];
//   for (int i = 0; i < iMaxEvals; i++) {
//     adMinQual[i] = dMax; // This pre-sets the min/max measured
//     adMaxQual[i] = dMin; // quality appropriately for use.
//     a2dBinFrac[i] = new double[iNBins];
//     for (int ii = 0; ii < iNB; ii++)
//       a2dBinFrac[i][ii] = 0.;
//   }
// }

Quality::~Quality()
{
  for (int i = 0; i < iMaxEvals; i++) 
    delete [] a2dBinFrac[i];
  delete [] a2dBinFrac;
  delete [] adSignifImp;
  delete [] aiNumQualVals;
  delete [] adMinQual;
  delete [] adMaxQual;
}

int Quality::iComputeBin(const double dQual) const
{
  int iBin = int( iNBins * (dQual - dMinQual) / (dMaxQual - dMinQual) );
  assert(iBin <= iNBins && iBin >= 0);
  if (iBin == iNBins) iBin = iNBins - 1;
  return (iBin);
}

void Quality::vEvaluate()
{
  SUMAA_LOG_EVENT_BEGIN(QUAL_ASSESS);
  if (iQualEvals >= iMaxEvals) {
    vMessage(0, "%s\n %s %d\n",
	     "Too many evaluations of mesh mesh quality; out of data space",
	     "Try initializing the Quality object with iMaxEvals > ",
	     iMaxEvals);
    return;
  }
  vMessage(1, "Evaluating quality of volume mesh...");
  int iSum = 0;
  double dMin = adMinQual[iQualEvals];
  double dMax = adMaxQual[iQualEvals];
  for (int i = pM->iNumCells() - 1; i >= 0; i--) {
    Cell *pC = pM->pCCell(i);
    assert(pC->qValid());
    if (pC->qDeleted()) continue;
    assert(pC->iFullCheck() == 1);
    int iN;
    double adResult[24];
    pvfQualityFunc(pC, &iN, adResult);
    for (int ii = 0; ii < iN; ii++) {
      double dQual = adResult[ii];
      dMin = min(dMin, dQual);
      dMax = max(dMax, dQual);
      int iBin = iComputeBin(dQual);
      a2dBinFrac[iQualEvals][iBin] ++;
      iSum++;
    }
  }
  aiNumQualVals[iQualEvals] = iSum;

  double dFrac = 1 / double(iSum);
  int ii;
  for (ii = 0; ii < iNBins; ii ++)
    a2dBinFrac[iQualEvals][ii] *= dFrac;
  adMinQual[iQualEvals] = dMin;
  adMaxQual[iQualEvals] = dMax;

  // Assess the statistical significance of the change in the quality
  // distribution between this evaluation and the previous one using the
  // Kolmogorov-Smirnov test for binned data.
  if (iQualEvals >= 1) {
    double dDiff = 0, dMaxDiff = 0;
    // Compute the difference in the cumulative distribution bin-by-bin,
    // and keep track of the maximum value.
    for (ii = 0; ii < iNBins; ii++) {
      dDiff = dDiff + a2dBinFrac[iQualEvals-1][ii] // Old value
	- a2dBinFrac[iQualEvals][ii]; // New value
      if (fabs(dDiff) > dMaxDiff)
	dMaxDiff = fabs(dDiff);
    }
    int iN1 = aiNumQualVals[iQualEvals-1];
    int iN2 = aiNumQualVals[iQualEvals];
    adSignifImp[iQualEvals] =
      1 - dQKolSmirn(dMaxDiff * sqrt(double(iN1*iN2)/double(iN1+iN2)));
  }

  iQualEvals++;
  
  vMessage(1, "done.\n");
  SUMAA_LOG_EVENT_END(QUAL_ASSESS);
}

double Quality::dEvaluate(const CellSkel * const pC) const
{
  assert(pC->qValid());
  assert(!pC->qDeleted());
  assert(pC->iFullCheck() == 1);
  int iN;
  double adResult[24];
  pvfQualityFunc(pC, &iN, adResult);

  assert(iN == 1);
  return (adResult[0]);
}

void Quality::vWriteToFile(char strQualFileName[]) const
{
  SUMAA_LOG_EVENT_BEGIN(OUTPUT);
  FILE *pF = fopen(strQualFileName, "w");
  if (pF == NULL) 
    vFatalError("Couldn't open output file for quality data.",
		"writing quality data to disk");

  // Begin with a header telling which quality measure and what the
  // min/max were for each pass.
  if (pvfQualityFunc == vMaxDihed) 
    fprintf(pF, "# Quality measure:  maximum dihedral angle\n"); 
  else if (pvfQualityFunc == vMinDihed)
    fprintf(pF, "# Quality measure:  minimum dihedral angle\n"); 
  else if (pvfQualityFunc == vAllDihed)
    fprintf(pF, "# Quality measure:  all dihedral angles\n"); 
  else if (pvfQualityFunc == vMaxSolid)
    fprintf(pF, "# Quality measure:  maximum solid angle\n"); 
  else if (pvfQualityFunc == vMinSolid)
    fprintf(pF, "# Quality measure:  minimum solid angle\n"); 
  else if (pvfQualityFunc == vAllSolid)
    fprintf(pF, "# Quality measure:  all solid angles\n"); 
  else if (pvfQualityFunc == vInToCircumSphere)
    fprintf(pF, "# Quality measure:  ratio of radii of insphere and circumsphere\n");
  else if (pvfQualityFunc == vSliverQ1)
    fprintf(pF, "# Quality measure:  Marcum's sliver measure\n"); 
  else if (pvfQualityFunc == vSkewnessQ2)
    fprintf(pF, "# Quality measure:  Marcum's skewness measure\n"); 
  else if (pvfQualityFunc == vVolArea)
    fprintf(pF, "# Quality measure:  Ratio of volume to surface area^1.5\n"); 
//   else if (pvfQualityFunc == vAreaPerim)
//     fprintf(pF, "# Quality measure:  Ratio of area to perimeter^2\n"); 

  fprintf(pF, "#\n");
  fprintf(pF, "# Qual eval        Min qual      Max qual    Conf imp\n");
//   fprintf(pF, "# Qual eval        Min qual      Max qual\n");
  int ii;
  for (ii = 0; ii < iQualEvals; ii++) {
    fprintf(pF, "#    %3d       %11.6f   %11.6f   %11.6f\n",
 	    ii, adMinQual[ii], adMaxQual[ii], adSignifImp[ii]);
//     fprintf(pF, "#    %3d       %11.6f   %11.6f\n",
// 	    ii, adMinQual[ii], adMaxQual[ii]);
  }

  fprintf(pF, "#\n");
  for (int i = 0 ; i < iNBins; i++) {
    fprintf(pF, "%8.3f", dMinQual + (double(i) + 0.5) / iNBins *
	    (dMaxQual - dMinQual));  
    for (ii = 0; ii < iQualEvals; ii++) 
      fprintf(pF, "%13.6f", a2dBinFrac[ii][i]);
    fprintf(pF, "\n");
  }
  fclose(pF);
  SUMAA_LOG_EVENT_END(OUTPUT);
}

void Quality::vImportQualityData(const double adData[])
{
  SUMAA_LOG_EVENT_BEGIN(QUAL_ASSESS);
  if (iQualEvals >= iMaxEvals) {
    vMessage(0, "%s\n %s %d\n",
	     "Too many evaluations of mesh mesh quality; out of data space",
	     "Try initializing the Quality object with iMaxEvals > ",
	     iMaxEvals);
    return;
  }
  vMessage(1, "Importing quality data...");
  int iSum = 0;
  double dMin = adMinQual[iQualEvals];
  double dMax = adMaxQual[iQualEvals];
  for (int i = pM->iNumCells() - 1; i >= 0; i--) {
    double dQual = adData[i];
    dMin = min(dMin, dQual);
    dMax = max(dMax, dQual);
    int iBin = iComputeBin(dQual);
    a2dBinFrac[iQualEvals][iBin] ++;
    iSum++;
  }
  aiNumQualVals[iQualEvals] = iSum;

  double dFrac = 1 / double(iSum);
  int ii;
  for (ii = 0; ii < iNBins; ii ++)
    a2dBinFrac[iQualEvals][ii] *= dFrac;
  adMinQual[iQualEvals] = dMin;
  adMaxQual[iQualEvals] = dMax;

  // Assess the statistical significance of the change in the quality
  // distribution between this evaluation and the previous one using the
  // Kolmogorov-Smirnov test for binned data.
  if (iQualEvals >= 1) {
    double dDiff = 0, dMaxDiff = 0;
    // Compute the difference in the cumulative distribution bin-by-bin,
    // and keep track of the maximum value.
    for (ii = 0; ii < iNBins; ii++) {
      dDiff = dDiff + a2dBinFrac[iQualEvals-1][ii] // Old value
	- a2dBinFrac[iQualEvals][ii]; // New value
      if (fabs(dDiff) > dMaxDiff)
	dMaxDiff = fabs(dDiff);
    }
    int iN1 = aiNumQualVals[iQualEvals-1];
    int iN2 = aiNumQualVals[iQualEvals];
    adSignifImp[iQualEvals] =
      1 - dQKolSmirn(dMaxDiff * sqrt(double(iN1*iN2)/double(iN1+iN2)));
  }

  iQualEvals++;
  
  vMessage(1, "done.\n");
  SUMAA_LOG_EVENT_END(QUAL_ASSESS);
}

