path: root/header.h

#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 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;

  for (int 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;
}

void report(io, markers, reporttriangles, reportneighbors, reportsegments,
            reportedges, reportnorms)
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");
  }
}

/*
 * ----------------------------------------------------------------------
 * 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;
}