#include "GR_assert.h"
#ifndef WIN32
#include <unistd.h>
#else /* WIN32 defined */
int getopt(int argc, char *argv[], char *opstring);
#endif

#include <vector>

#include <stdio.h>
#include <string.h>
#include "GR_Mesh2D.h"
#include "GR_events.h"
#include "GR_misc.h"
#include "GR_ADT.h"

#ifdef IRIX
#include <getopt.h>
#endif

static void vUsage()
{
  vMessage(0, "Usage: scat2d -i basefilename [options] \n");
  vMessage(0, "   -m        Output 2D mesh [not done by default] \n\n");
  exit(1);
}

int main(int iNArg, char *apcArgs[]) {
  char strBaseFileName[FILE_NAME_LEN], strDataFileName[FILE_NAME_LEN],
      strPointFileName[FILE_NAME_LEN], strOutputFileName[FILE_NAME_LEN];
  int iPChar;

  vGRUMMPInit("scat2d");
  
  bool qError = true, qSaveMesh = false;
  while ((iPChar = getopt(iNArg, apcArgs, "i:m")) != EOF) {
    switch (iPChar) {
    case 'i':
      // Grab data file name in a way that prevents buffer overruns
      strncpy(strBaseFileName, optarg, FILE_NAME_LEN-1);
      {
	int iDum;
	for (iDum = 0; iDum < FILE_NAME_LEN && strBaseFileName[iDum];
	     iDum++) {}
	if (iDum == FILE_NAME_LEN) // No null; file name too long
	  vFatalError("File name too long", "command line");
      }
      qError = false;
      break;
    case 'm':
      qSaveMesh = true;
      break;
    default:
      vUsage();
    }
  }
  if (qError) vUsage();

  // Open log file
  vOpenMessageFile(strBaseFileName);
  
  // Open input file
  vMakeFileName(strDataFileName, "%s.data", strBaseFileName,
		"main() [scat2d.C]");  
  FILE *pFileIn = fopen(strDataFileName, "r");
  if (pFileIn == NULL)
    vFatalError("Scattered data file does not exist or could not be opened",
		"reading 2D scattered data");

  // Read the scattered data points and the data
  int iNVerts, iDataPerPoint;
  const int iBufSize = 2048;
  char *acBuffer = new char[iBufSize];

  vGetLineOrAbort(acBuffer, iBufSize, pFileIn);
  if (2 != sscanf(acBuffer, "%d %d", &iNVerts, &iDataPerPoint))
    vFatalError("Problem with mesh size info",
		"reading 2D scattered data");

  VAVertex VAInput;
  VAInput.vSetup(iNVerts);

  double **a2dData = new double*[iNVerts];
  int i;
  for (i = 0; i < iNVerts; i++) {
    a2dData[i] = new double[iDataPerPoint];
  }
  for (i = 0; i < iNVerts; i++) {
    double adCoord[2];
    char *pcData = acBuffer;
    char strErrorString[200];
    sprintf(strErrorString, "Insufficient data on line %d of file %s.",
	    i, strDataFileName);

    vGetLineOrAbort(acBuffer, iBufSize, pFileIn);
    vGetDoubleFromBuffer(&pcData, &(adCoord[0]), strErrorString, "main()");
    vGetDoubleFromBuffer(&pcData, &(adCoord[1]), strErrorString, "main()");
    VAInput.pTEntry(i)->vSetCoords(2, adCoord);

    for (int iData = 0; iData < iDataPerPoint; iData++)
      vGetDoubleFromBuffer(&pcData, &(a2dData[i][iData]),
			   strErrorString, "main()");
  }

  // Create a 2D mesh from the data points
  Mesh2D M2D(VAInput);
  assert(M2D.qSimplicial());

  // Create a tree of the cells from that mesh to make interpolation
  // efficient
  double (*a2dBBox)[4] = new double[M2D.iNumCells()][4];
  int iCell;
  for (iCell = 0; iCell < M2D.iNumCells(); iCell++) {
    a2dBBox[iCell][0] =  1.e30;
    a2dBBox[iCell][1] = -1.e30;
    a2dBBox[iCell][2] =  1.e30;
    a2dBBox[iCell][3] = -1.e30;
    Cell *pC = M2D.pCCell(iCell);
    for (int iV = 0; iV < 3; iV++) {
      const double *adCoord = pC->pVVert(iV)->adCoords();
      a2dBBox[iCell][0] = min(a2dBBox[iCell][0], adCoord[0]);
      a2dBBox[iCell][1] = max(a2dBBox[iCell][1], adCoord[0]);
      a2dBBox[iCell][2] = min(a2dBBox[iCell][2], adCoord[1]);
      a2dBBox[iCell][3] = max(a2dBBox[iCell][3], adCoord[1]);
    }
  }
  ADT BBoxTree(M2D.iNumCells(), 2, ADT::eBBoxes,
	       reinterpret_cast<double*>(a2dBBox));
  delete [] a2dBBox;

  // Read the target data points, interpolate on the fly, and spit out
  // the results.
  vMakeFileName(strPointFileName, "%s.points", strBaseFileName,
		"main() [scat2d.C]");  
  FILE *pFilePts = fopen(strPointFileName, "r");
  if (pFilePts == NULL)
    vFatalError("Target point file does not exist or could not be opened",
		"reading 2D target points for scattered data interp");

  int iNTargetVerts;
  vGetLineOrAbort(acBuffer, iBufSize, pFilePts);
  if (1 != sscanf(acBuffer, "%d", &iNTargetVerts))
    vFatalError("Problem with mesh size info",
		"reading 2D target points for scattered data interp");

  FILE *pF = NULL;
  vMakeFileName(strOutputFileName, "%s.out", strBaseFileName,
		"main() [scat2d.C]");  
  
  if (NULL == (pF = fopen(strOutputFileName, "w"))) 
    vFatalError("Couldn't open output file", "2D scattered data interp");

  std::vector<int> veciQueryResult;
  for (i = 0; i < iNTargetVerts; i++) {
    double adCoord[2];
    char *pcData = acBuffer;
    char strErrorString[200];
    sprintf(strErrorString, "Insufficient data on line %d of file %s.",
	    i, strPointFileName);

    // Read the target data points
    vGetLineOrAbort(acBuffer, iBufSize, pFilePts);
    vGetDoubleFromBuffer(&pcData, &(adCoord[0]), strErrorString, "main()");
    vGetDoubleFromBuffer(&pcData, &(adCoord[1]), strErrorString, "main()");

    // Find cell (nearly?) containing the target point
    double dEps = 1.e-8;
    do {
      double adBB[] = {adCoord[0] - dEps,
		       adCoord[0] + dEps,
		       adCoord[1] - dEps,
		       adCoord[1] + dEps};
      veciQueryResult = BBoxTree.veciRangeQuery(adBB);
      dEps *= 10;
    } while (veciQueryResult.empty());
    // Find out which cell contains the vert, or which one comes closest
    // if none contain it.  "Closest" is defined in the sense of having
    // barycentric coordinates that violate the range (0,1) by as little
    // as possible.

    double dBestValue = 1.e50;
    int iBestCell = -1;
    double adBestBary[3];
    Cell *pC;
    int iQ;
    for (iQ = veciQueryResult.size() - 1; iQ >= 0 && dBestValue > 0.5; iQ--) {
      // Get barycentric coords.  If best_value <(=) 0.5 on completion
      // then the vert lies within (on an edge of) the cell in list
      iCell = veciQueryResult[iQ];
      pC = M2D.pCCell(iCell);
      double adBary[3];
      dynamic_cast<TriCell*>(pC)->vBarycentrics(adCoord, adBary);

      double dThisValue = max(max(fabs(adBary[0]-0.5),
				  fabs(adBary[1]-0.5)),
			      fabs(adBary[2]-0.5));

      if (dThisValue < dBestValue) {
	iBestCell = iCell;
	dBestValue = dThisValue;
	adBestBary[0] = adBary[0];
	adBestBary[1] = adBary[1];
	adBestBary[2] = adBary[2];
      }
    } // Done checking all cells which have bboxes overlapping with the
      // vertex 
    assert(iBestCell >= 0);
    pC = M2D.pCCell(iBestCell);

    // Spit out the points...
    fprintf(pF, "%g %g", adCoord[0], adCoord[1]);

    // Interpolate data to the target points
    for (int iData = 0; iData < iDataPerPoint; iData++) {
      double dTargetData = 0;
      for (int ii = 0; ii < 3; ii++) {
	int iVert = M2D.iVertIndex(pC->pVVert(ii));
	dTargetData += a2dData[iVert][iData] * adBestBary[ii];
      }
      fprintf(pF, " %g", dTargetData);
    }
    fprintf(pF, "\n");
  }
  
  // Now delete a bunch of arrays.
  for (i = 0; i < iNVerts; i++) {
    delete [] a2dData[i];
  }
  delete [] a2dData;

  fclose(pF);
  fclose(pFilePts);
  fclose(pFileIn);

  if (qSaveMesh) {
    vMessage(1,"Writing 2D mesh file...\n");
    vWriteFile_Mesh2D(M2D, strBaseFileName);
    vMessage(1,"Done\n");
  }
  
  return (0);
}

