#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#ifdef SINGLE
#define REAL float
#else
#define REAL double
#endif

#include "triangle.h"

#ifndef DEBUG_IO
#define DEBUG_IO 0
#endif

#ifndef VERBOSE
#define VERBOSE 0
#endif

/*
 * ----------------------------------------------------------------------
 * wavelet functions
 * ----------------------------------------------------------------------
 */

int get_metadata(char *datafilename, int *nrow, int *ncol, double *xllcorner, double *yllcorner, double *cellsize, double *nodata);

int get_data(char *datafilename, int nrow, int ncol, int ncol_padded, double *data);

int compress(double *data, double *m, int nrow, int ncol);

int decompress(double *m, int nrow, int ncol);

int hard_threshold(double *m, int *flag, int n, int nlevels, int nrow, int ncol, double threshold);

int multiresolution_analysis(int argc, char** argv);

int multiresolution_meshing(int argc, char** argv);

/*
 * ----------------------------------------------------------------------
 */

int sign(double x)
{
  int flag = (x > 0) ? 1 : -1;
  if (x == 0) {
    flag = 0;
  }
  return flag;
}

double max(double a, double b)
{
  double results = (a > b) ? a : b;
  return results;
}

double min(double a, double b)
{
  double results = (a < b) ? a : b;
  return results;
}

int ipow(int base, int exponent)
{
  int result = 1;
  for (int i = exponent; i > 0; i--) {
    result = result * base;
  }
  return result;
}

int get_wavelet_matrix(double *m, int dim)
{
  int len = (int) (dim / 2);
  
  int r1[len]; /* rows for low frequency */
  int r2[len]; /* rows for high frequency */
  int c1[len]; /* columns */
  int c2[len]; /* columns */

  for (int i = 0; i < len; i ++) {
    r1[i] = i;
    r2[i] = len + i;
    c1[i] = 2 * i;
    c2[i] = 2 * i + 1;
    
    m[c1[i] + dim * r1[i]] = 0.5 * sqrt(2.0);
    m[c2[i] + dim * r1[i]] = 0.5 * sqrt(2.0);
    
    m[c1[i] + dim * r2[i]] = 0.5 * sqrt(2.0);
    m[c2[i] + dim * r2[i]] = - 0.5 * sqrt(2.0);
  }
  
  return 0;
}

int matmul(double* matrix1, int nrow1, int ncol1, double* matrix2, int nrow2, int ncol2, double* result) {

  if (ncol1 != nrow2) {
    fscanf(stderr, "[ERR] Dimensions of matrices don't match: %d-by%d with %d-by%d\n", &nrow1, &ncol1, &nrow2, &ncol2);
    return -1;
  }
  
  double sum = 0.0;
  int i;

#pragma omp parallel
  {
#pragma omp for private(i) reduction(+:sum)
    for (i = 0; i < nrow1; i ++) {
      for (int j = 0; j < ncol2; j ++) {
	for (int k = 0; k < ncol1; k ++) {
	  sum += matrix1[k + i * ncol1] * matrix2[j + k * ncol2];
	}
	result[j + ncol2 * i] = sum;
	sum = 0.0;
      }
    }
  }

  return 0;
}

int transpose(double* matrix, int nrow, int ncol) {

  double *_tmp = (double*) calloc(ncol * nrow, sizeof(double));

  for (int j = 0; j < nrow; j ++) {
    for (int i = 0; i < ncol; i ++) {
      _tmp[j + nrow * i] = matrix[i + ncol * j];
    }
  }

  memcpy(matrix, _tmp, nrow * ncol * sizeof(double));

  free(_tmp);
  
  return 0;
}

/*
 * ----------------------------------------------------------------------
 * meshing functions
 * ----------------------------------------------------------------------
 */

int get_border(double *data, double nodata, int nrow, int ncol)
{
  /* Square tracing */

  FILE *fp = fopen("output/boundary.asc", "w");

  int i = 0; /* x direction */
  int j = 0; /* y direction */

  double val;

  /* Bootstrap */

  val = data[i + ncol * j];

  while (val <= nodata) {

    i = i + 1;
    if (i == ncol) {
      i = 0;
      j = j + 1;
    }
    
    val = data[i + ncol * j];
    
  }

  int i0 = i;
  int j0 = j;

  /*
   *        1
   *        ^
   *        |
   * 2 < - - - - > 0
   *        |
   *        v
   *        3
   */
  
  int dir = 0;
  
  /* Detected first pixel of the boundary */
  
  while (1) {
    
    if (i < 0 || j < 0 || i >= ncol || j >= nrow) {
      val = nodata;
    } else {
      val = data[i + ncol * j];
    }
    
    if (val > nodata) {
      dir = (dir + 1) % 4; /* go left */
      fprintf(fp, "%d %d %f\n", i, j, val);
    } else {
      dir = (dir + 3) % 4; /* go right */
    }

    switch (dir) {
    case 0:
      i = i + 1;
      break;
    case 1:
      j = j + 1;
      break;
    case 2:
      i = i - 1;
      break;
    case 3:
      j = j - 1;
    }
    
    if (i == i0 && j == j0) {
      break;
    }
    
  }

  fclose(fp);

  return 0;
}

void write_avs(struct triangulateio *io, int counter, double *elevations, int ncol, int nrow, double xllcorner, double yllcorner, double cellsize, double nodata)
{

  char fname[1024];
  sprintf(fname, "output/mesh%d.inp", counter);
  
  FILE *fp;
  fp = fopen(fname, "w");

  fprintf(fp, "%d %d %d %d %d\n", io->numberofpoints, io->numberoftriangles, 0, 0, 0);

  for (int i = 0; i < io->numberofpoints; ++ i) {

    /* map elevation to the grid */
    double x = io->pointlist[2*i];
    double y = io->pointlist[2*i+1];

    int ix = (int) floor((x - xllcorner) / cellsize);
    int iy = nrow - (int) floor((y - yllcorner) / cellsize);

    double z = elevations[ix + ncol * iy];
    if (z == nodata) {
      int ix0 = (ix + 1) < ncol ? ix + 1 : ix;
      int ix1 = (ix - 1) > 0 ? ix - 1 : ix;
      int iy0 = (iy + 1) < nrow ? iy + 1 : iy;
      int iy1 = (iy - 1) > 0 ? iy - 1 : iy;

      /* indices of von neumann neighborhood */
      int neumann[4] = {ix0 + ncol * iy, ix1 + ncol * iy, ix + ncol * iy0, ix + ncol * iy1};

      for (int j = 0; j < 4; j ++) {
	double znew = elevations[neumann[j]];
	if (znew != nodata) {
	  z = znew;
	  break;
	}
      }

      int moore[4] = {ix0 + ncol * iy0, ix0 + ncol * iy1, ix1 + ncol * iy0, ix1 + ncol * iy1};

      for (int j = 0; j < 4; j ++) {
	double znew = elevations[moore[j]];
	if (znew != nodata) {
	  z = znew;
	  break;
	}
      }
      
    }
    
    fprintf(fp, "%08d %f %f %f\n", i+1, io->pointlist[2*i], io->pointlist[2*i+1], z);

  }

  for (int i = 0; i < io->numberoftriangles; ++ i) {
    fprintf(fp, "%08d %8d tri %8d %8d %8d\n", i+1, 1, io->trianglelist[3*i], io->trianglelist[3*i+1], io->trianglelist[3*i+2]);
  }

  fclose(fp);
  
}

void report(struct triangulateio *io, int markers, int reporttriangles,
	    int reportneighbors, int reportsegments, int reportedges,
	    int reportnorms)
{
  
  int i, j;

  for (i = 0; i < io->numberofpoints; i++) {
    printf("Point %4d:", i);
    for (j = 0; j < 2; j++) {
      printf("  %.6g", io->pointlist[i * 2 + j]);
    }
    if (io->numberofpointattributes > 0) {
      printf("   attributes");
    }
    for (j = 0; j < io->numberofpointattributes; j++) {
      printf("  %.6g",
             io->pointattributelist[i * io->numberofpointattributes + j]);
    }
    if (markers) {
      printf("   marker %d\n", io->pointmarkerlist[i]);
    } else {
      printf("\n");
    }
  }
  printf("\n");

  if (reporttriangles || reportneighbors) {
    for (i = 0; i < io->numberoftriangles; i++) {
      if (reporttriangles) {
        printf("Triangle %4d points:", i);
        for (j = 0; j < io->numberofcorners; j++) {
          printf("  %4d", io->trianglelist[i * io->numberofcorners + j]);
        }
        if (io->numberoftriangleattributes > 0) {
          printf("   attributes");
        }
        for (j = 0; j < io->numberoftriangleattributes; j++) {
          printf("  %.6g", io->triangleattributelist[i *
                                         io->numberoftriangleattributes + j]);
        }
        printf("\n");
      }
      if (reportneighbors) {
        printf("Triangle %4d neighbors:", i);
        for (j = 0; j < 3; j++) {
          printf("  %4d", io->neighborlist[i * 3 + j]);
        }
        printf("\n");
      }
    }
    printf("\n");
  }

  if (reportsegments) {
    for (i = 0; i < io->numberofsegments; i++) {
      printf("Segment %4d points:", i);
      for (j = 0; j < 2; j++) {
        printf("  %4d", io->segmentlist[i * 2 + j]);
      }
      if (markers) {
        printf("   marker %d\n", io->segmentmarkerlist[i]);
      } else {
        printf("\n");
      }
    }
    printf("\n");
  }

  if (reportedges) {
    for (i = 0; i < io->numberofedges; i++) {
      printf("Edge %4d points:", i);
      for (j = 0; j < 2; j++) {
        printf("  %4d", io->edgelist[i * 2 + j]);
      }
      if (reportnorms && (io->edgelist[i * 2 + 1] == -1)) {
        for (j = 0; j < 2; j++) {
          printf("  %.6g", io->normlist[i * 2 + j]);
        }
      }
      if (markers) {
        printf("   marker %d\n", io->edgemarkerlist[i]);
      } else {
        printf("\n");
      }
    }
    printf("\n");
  }
}