/* See {salamic_planes.h}. */
/* Last edited on 2021-06-27 12:28:07 by jstolfi */
  
#define _GNU_SOURCE
#include <stdio.h>
#include <assert.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>

#include <bool.h>
#include <vec.h>
#include <jsmath.h>
#include <affirm.h>
#include <fget.h>

#include <salamic_planes.h>
/* #include <stmesh_utils.h> */

int_vec_t salamic_planes_get_uniform(float startZ, float deltaZ, float eps, int minZ, int maxZ)
  {
    demand(deltaZ >= 2*eps, "plane Z spacing is too small");
    demand((minZ & 1 ) == 0, "{minZ} must be even");
    demand((maxZ & 1 ) == 0, "{maxZ} must be even");
    int32_t sZ = (int32_t)iroundfrac(startZ, eps, 2, 1, INT32_MAX); /* Round to odd. */
    int32_t dZ = (int32_t)iroundfrac(deltaZ, eps, 2, 0, INT32_MAX); /* Round to even. */
    assert(dZ >= 2);
    
    /* Get the quantized {Z-coordinate {bZ} of the first plane in the range: */
    int32_t uZ = (sZ - minZ) % dZ;
    if (uZ < 0) { uZ += dZ; }
    int32_t bZ = minZ + uZ;
    assert((bZ - sZ) % dZ == 0);
    assert(bZ >= minZ);
    assert(bZ - minZ < dZ);
    
    /* Compute the number of planes: */
    int32_t np = (maxZ - 1 - bZ)/dZ + 1;
    int_vec_t planeZ = int_vec_new(np);

    int32_t i;
    for (i = 0; i < np; i++) 
      { planeZ.e[i] = bZ; 
        bZ += dZ;
      }
    assert(bZ > maxZ);

    fprintf(stderr, "generated %d planes with Z in %d .. %d", np, planeZ.e[0], planeZ.e[np-1]);
    fprintf(stderr, " [%.4f _ %.4f mm]\n", planeZ.e[0]*(double)eps, planeZ.e[np-1]*(double)eps);
    
    return planeZ;
  }

int_vec_t salamic_planes_get_adaptive(char *ZFile, float eps, int minZ, int maxZ)
  {
    FILE *rd = fopen(ZFile, "rt");
    
    /* Get the number of planes: */
    int32_t line = 1;
    int32_t np_file = (int32_t)fget_int(rd);  /* Number of planes in file. */
    fget_eol(rd); line++;
    int_vec_t planeZ = int_vec_new(np_file);
    

    /* Read the {Z}-coords: */
    int32_t minZ_file = INT32_MAX; /* Minimum quantized {Z} in file. */
    int32_t maxZ_file = INT32_MIN; /* Maximum quantized {Z} in file. */
    int32_t i;
    int32_t np = 0;  /* Number of planes in range. */
    for(i = 0; i < np_file; i++)
      { float fZ = (float)fget_double(rd);
        int32_t pZ = (int32_t)iroundfrac(fZ, eps, 2, 1, INT32_MAX);
        if (pZ < minZ_file) { minZ_file = pZ; }
        if (pZ > maxZ_file) { maxZ_file = pZ; }
        if ((pZ >= minZ) && (pZ <= maxZ))
          { planeZ.e[np] = pZ;
            if ((np > 0) && (planeZ.e[np-1] >= pZ))
              { fprintf(stderr, "%s:%d: ** error: slicing planes not increasing or too close\n", ZFile, line);
                exit(1);
              }
            np++;
          }
        fget_eol(rd); line++;
      }
    fclose(rd);
    int_vec_trim(&planeZ, np);
    
    fprintf(stderr, "read %d planes from file with Z in %d .. %d", np_file, minZ_file, maxZ_file);
    fprintf(stderr, " [%.4f _ %.4f mm]\n", minZ_file*(double)eps, maxZ_file*(double)eps);
    fprintf(stderr, "kept %d planes with Z in %d .. %d", np, planeZ.e[0], planeZ.e[np-1]);
    fprintf(stderr, " [%.4f _ %.4f mm]\n", planeZ.e[0]*(double)eps, planeZ.e[np-1]*(double)eps);
    
    return planeZ;
  }