/* See {salamic_closer_Stolfi.h} */
/* Last edited on 2015-11-16 01:26:11 by stolfilocal */

#define _GNU_SOURCE
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>

#include <bool.h>
#include <vec.h>
#include <jsfile.h>
#include <affirm.h>

#include <stmesh.h>

#include <salamic_closer.h>

#include <salamic_closer_Stolfi.h>

/* AUXILIARY TABLES */

typedef struct uxep_t { stmesh_edge_unx_t e[2]; } uxep_t;
  /* A pair of indices of unoriented edges. Used below 
    to represent two edges that are on the same face and are 
    crossed by the slicing plane. */

typedef struct jfp_t { uint32_t c[2]; } jfp_t;
  /* A pair of indices into the list of cut faces. */
  
#define UNX_NULL UINT32_MAX
  /* An index value used as {NULL}. */ 

/* 
  The procedures {salamic_closer_Stolfi_get_path} and
  {salamic_closer_Stolfi_extend_path} use the auxiliary tables
  {uxep[0..mf-1]}, {efirstf[0..ne-1]}, and {fnextf[0..mf-1]}, where 
  {ne=mesh.ne}.

  Expects that {uxep[0..mf-1]} be a list of pairs of indices of mesh
  edges, in arbitrary order. Each pair {p=uxep[jf]} in this
  list should correspond to a triangle of the mesh that intercepts the
  slicing plane, specifically through mesh edges with indices {p.e[0]} and
  {p.e[1]}.

  All the pairs {uxep[k]} that that include the same edge with index
  {uxe} are supposed to be linked in a list defined by the vectors
  {efirstf} and {fnextf}. Namely, {p=uxep[efirstf[uxe]]} is one of the
  pairs in {uxep} that use that edge, that is, that has {p.e[0]=uxe} or
  {p.e[1]=uxe}. If there is no such pair, then {efirstf[uxe]} is 
  {UNX_NULL}.
  
  Also, for all {jf} in {0..nf-1}, {uxep[fnextf[jf].c[r]]}
  is the next pair in {uxep} that uses the edge {uxep[jf].e[r]}.
  If there is no such next pair, {fnextf[jf].c[r]} is {UNX_NULL}.

  However, these lists should include only the pairs {uxep[jf]} that
  have not been included yet in any loop. Whenever these procedures use
  a pair {uxep[jf] = (e0,e1)}, they remove that entry from the lists of
  both edges {e0} and {e1}, and set {uxep[jf]} to an invalid index pair.
  However, {f[jf]} is left undisturbed. */


/* INTERNAL PROTOTYPES */

void salamic_closer_Stolfi_get_path
  ( stmesh_t mesh, 
    int32_t pZ, 
    uint32_t mf,
    stmesh_face_t f[],
    uxep_t uxep[],
    uint32_t efirstf[],
    jfp_t fnextf[],
    uint32_t kf, 
    uint32_t *mepP,
    stmesh_edge_t ep[],
    uint32_t mep_max, 
    salamic_stats_t *st
  );
  /* Extracts a closed loop or maximal open path from a cross-section of
    the mesh {mesh} with a slicing plane at quantized {Z}-coordinate
    {pZ}.
    
    The extracted path will include the segment defined by the pair
    {uxep[kf]}, that should be unused. The path is extended in both
    directions from that segment, with zero or more still-unused faces,
    until it closes or reaches an open edge of the mesh. The procedure
    removes from the lists defined by {efirstf,fnextf} any pair {uxep[jf]}
    that it included in the path from the linked lists..
    
    The path (or loop) is returned as a list {ep[0..mep-1]} of mesh edges
    in {mesh.e}. The intersection of each edge with the slicing plane is
    a vertex of the path. The path is a closed loop if and only if {e[0]
    == e[me-1]}. The value of {mep} is returned in {*mepP}, whose previous
    value is ignored.  The array {ep} is assumed to have been allocated with 
    {mep_max} elements. */

void salamic_closer_Stolfi_extend_path
  ( stmesh_t mesh, 
    int32_t pZ, 
    uint32_t mf,
    stmesh_face_t f[],
    uxep_t uxep[],
    uint32_t efirstf[],
    jfp_t fnextf[],
    uint32_t *mepP,
    stmesh_edge_t ep[],
    uint32_t mep_max
  );
  /* Extends a path from a cross-section of the mesh {mesh} with a
    slicing plane at quantized {Z}-coordinate {pZ}, until it closes or
    reaches an open edge.
    
    Assumes that the path's vertices obtained so far are the
    intersections of the slicing plane with the edges {ep[0..mep-1]} where
    {mep = (*mepP)}. The path must have at least two vertices. Expects
    that the linked lists defined by {efirstf} and {fnextf} contain all
    pairs in {uxep[0..mf-1]} that correspond to triangles of the mesh
    that intercept the slicing plane, and have not been used in previous
    loops or paths.
    
    The path is extended forwards, appending more edges to {ep} and
    incrementing {mep}, until it closes or reaches an open edge of the
    mesh. If it closes, it will return with {ep[0] == e[mep-1]}. The
    procedure removes from the lists defined by {efirstf,fnextf} any
    pair {uxep[jf]} that it added to the path. The array {ep} is assumed
    to have been allocated with {mep_max} elements. */

void salamic_closer_Stolfi_remove_uxep
  ( stmesh_t mesh,
    uint32_t mf, 
    uxep_t uxep[], 
    uint32_t efirstf[],
    jfp_t fnextf[], 
    uint32_t kf
  );
  /* Removes the edge pair {uxep[kf]} from the lists of its two edges.
    Also sets {uxep[kf]} to a pair of null indices. */
    
/* IMPLEMENTATIONS */

void salamic_closer_Stolfi_close
  ( stmesh_t mesh,
    int32_t pZ,
    uint32_t mf,
    stmesh_face_t f[], 
    uint32_t *ncP, 
    uint32_t estart[],
    stmesh_edge_t e[],
    salamic_stats_t *st
  )
  {
    uint32_t ne = stmesh_edge_count(mesh);
    
    /* Allocate the auxiliary tables: */
    uxep_t *uxep = notnull(malloc(mf*sizeof(uxep_t)), "no mem"); /* Edge pairs from each cut face. */
    uint32_t *efirstf = notnull(malloc(ne*sizeof(uint32_t)), "no mem"); /* First face with same edge. */
    jfp_t *fnextf = notnull(malloc(mf*sizeof(jfp_t)), "no mem"); /* Next face with same edge. */
    
    /* Initialize the lists to empty: */
    stmesh_edge_unx_t uxe;
    for (uxe = 0; uxe < ne; uxe++) { efirstf[uxe] = UNX_NULL; }
    
    /* Collect the edge pairs from {f} and build the lists: */
    uint32_t jf; /* Index into {f}, {uxep}, or {fnextf}. */
    for(jf = 0; jf < mf; jf++) 
      { /* Get the two edges {e[0..1]} of face {f[jf]} intercepted by the plane: */
        stmesh_edge_t e[2];
        stmesh_face_get_sliced_sides(mesh, f[jf], pZ, e);
        
        /* Prepare the pair {p = uxep[jf]}: */
        uxep_t *p = &(uxep[jf]);
        jfp_t *fnx = &(fnextf[jf]);
        int r;
        for (r = 0; r < 2; r++) 
          { /* Get the index of one of the two edges: */
            uxe = stmesh_edge_get_unx(mesh, e[r]);
            p->e[r] = uxe;
            /* Insert the pair {uxep[jf]} in the corresponding list: */
            fnx->c[r] = efirstf[uxe];
            efirstf[uxe] = jf;
          }
      }
    
    uint32_t me_max = 2*mf; /* Assumes that {e} was allocated with {me_max} elements. */
    uint32_t kf = 0;        /* Next segment that may be the start of a new path/loop is {f[kf]}. */
    uint32_t nc = 0;        /* Number of loops/paths found. */
    uint32_t ini = 0;       /* Start of next path in {e} array. */
    estart[0] = ini;
    while (TRUE) 
      { /* Look for the first unused segment: */
        while ((kf < mf) && (uxep[kf].e[0] == UNX_NULL)) { kf++; }
        if (kf >= mf) { break; }
        /* Collect another loop/path that includes this segment, mark used faces: */
        uint32_t mep;
        uint32_t mep_max = me_max - ini;
        stmesh_edge_t *ep = &(e[ini]);
        salamic_closer_Stolfi_get_path(mesh, pZ, mf, f, uxep, efirstf, fnextf, kf, &mep, ep, mep_max, st);
        assert((mep >= 2) && (mep <= mep_max));
        /* Got one more loop, prepare for the next one: */
        nc++;
        ini = ini + mep;
        estart[nc] = ini;
      }
    free(fnextf);
    free(efirstf);
    free(uxep);
    
    (*ncP) = nc;
  }

void salamic_closer_Stolfi_get_path
  ( stmesh_t mesh, 
    int32_t pZ, 
    uint32_t mf,
    stmesh_face_t f[],
    uxep_t uxep[],
    uint32_t efirstf[],
    jfp_t fnextf[],
    uint32_t kf, 
    uint32_t *mepP,
    stmesh_edge_t ep[], 
    uint32_t mep_max,
    salamic_stats_t *st
  )
  {
    bool_t verbose = FALSE;
    bool_t debug = FALSE;
    
    uint32_t mep = 0;
    
    /* Get the first face {fk}, its pair {pk}, two mesh edges {e[0],e[1]}: */
    assert(mep_max >= 2);
    demand(kf < mf, "invalid start face index");
    stmesh_face_t fk = f[kf];
    uxep_t pk = uxep[kf];
    int r;
    for (r = 0; r < 2; r++) { ep[r] = stmesh_get_edge(mesh, pk.e[r]); }
    mep = 2;

    if (debug) 
      { stmesh_face_unx_t uxfk = stmesh_face_get_unx(mesh, fk);
        fprintf(stderr, "starting path with face %u mesh edges %u and %u\n", uxfk, pk.e[0], pk.e[1]);
      }
    salamic_closer_Stolfi_remove_uxep(mesh, mf, uxep, efirstf, fnextf, kf);
    
    /* Collect a path whose first two edges are {e[0],e[1]}: */
    salamic_closer_Stolfi_extend_path(mesh, pZ, mf, f, uxep, efirstf, fnextf, &mep, ep, mep_max);
    
    if (ep[0] != ep[mep-1])
      { /* Path did not close: reverse it and extend from other end: */
        if (debug) { fprintf(stderr, "  got %d vertices, not closed: reversing...\n", mep); }
        uint32_t iv, jv;
        for (iv = 0, jv = mep-1; iv < jv; iv++, jv--)
          { stmesh_edge_t t = ep[iv]; ep[iv] = ep[jv]; ep[jv] = t; }
        if (debug) { fprintf(stderr, "  extending the other end...\n"); }
        salamic_closer_Stolfi_extend_path(mesh, pZ, mf, f, uxep, efirstf, fnextf, &mep, ep, mep_max);
        assert(ep[0] != ep[mep-1]);
      }
      
    if (verbose) 
      { assert(mep >= 2);
        stmesh_edge_unx_t uxe0 = stmesh_edge_get_unx(mesh, ep[0]);
        stmesh_edge_unx_t uxe1 = stmesh_edge_get_unx(mesh, ep[1]);
        stmesh_edge_unx_t uxef = stmesh_edge_get_unx(mesh, ep[mep-1]);
        fprintf(stderr, "got path with %6d corners - edges %u", mep, uxe0);
        if (mep >= 3) { fprintf(stderr, " %u", uxe1); }
        fprintf(stderr, " ... %u\n", uxef);
      }
      
    (*mepP) = mep;
  }
    
void salamic_closer_Stolfi_extend_path
  ( stmesh_t mesh, 
    int32_t pZ, 
    uint32_t mf,
    stmesh_face_t f[],
    uxep_t uxep[],
    uint32_t efirstf[],
    jfp_t fnextf[],
    uint32_t *mepP,
    stmesh_edge_t ep[],
    uint32_t mep_max
  )
  { 
    uint32_t mep = (*mepP);
    
    stmesh_edge_t e_ini = ep[0];    /* First vertex of path. */
    stmesh_edge_t e_fin = ep[mep-1]; /* Current last vertex of path. */
    while(TRUE)
      { 
        /* Look for another unused face {f[jf]} that uses edge {e_fin}: */
        stmesh_edge_unx_t uxe_fin = stmesh_edge_get_unx(mesh, e_fin);
        uint32_t jf = efirstf[uxe_fin];
        if (jf == UNX_NULL)
          { /* Reached border of mesh: */
            break;
          }
        else
          { /* Unused face {f[jf]} uses this edge, get the other edge {e_nxt}: */
            assert(jf < mf);
            uxep_t *pj = &(uxep[jf]);
            stmesh_edge_unx_t uxe_nxt; /* Index of edge that defines the next vertex of path: */
            if (pj->e[0] == uxe_fin)
              { uxe_nxt = pj->e[1]; }
            else if (pj->e[1] == uxe_fin)
              { uxe_nxt = pj->e[0]; }
            else
              { assert(FALSE); /* Edge {uxe_fin} must be in every pair of its list. */ }
            stmesh_edge_t e_nxt = stmesh_get_edge(mesh, uxe_nxt); /* Next vertex of path. */
            
            /* Remove the face {f[jf]} from the lists: */
            salamic_closer_Stolfi_remove_uxep(mesh, mf, uxep, efirstf, fnextf, jf);

            /* Append the new {e_nxt} to the path: */
            demand (mep < mep_max, "too many vertices in loop/path");
            ep[mep] = e_nxt;
            mep++;

            /* Prepare for next iteration: */
            e_fin = e_nxt;

            /* Check for closure: */
            if (e_fin == e_ini) { break; }
          }
     }
      
    (*mepP) = mep;
  }

void salamic_closer_Stolfi_remove_uxep
  ( stmesh_t mesh,
    uint32_t mf, 
    uxep_t uxep[], 
    uint32_t efirstf[],
    jfp_t fnextf[], 
    uint32_t kf
  )
  { 
    uint32_t ne = stmesh_edge_count(mesh);
    
    int r;
    for (r = 0; r < 2; r++)
      { stmesh_edge_unx_t uxer = uxep[kf].e[r];
        assert(uxer < ne); 
        assert(efirstf[uxer] != UNX_NULL);
        
        uint32_t jf = efirstf[uxer];
        if (jf == kf)
          { /* Pair {uxep[kf]} is the first in the list: */
            efirstf[uxer] = fnextf[kf].c[r];
          }
        else
          { /* Find the pair {uxep[jf_ant]} preceding {uxep[kf]} in the list of the edge {uxer}: */
            uint32_t jf_ant = UNX_NULL; /* Previous pair where this edge appeared. */
            int s_ant = -1; /* Which of the two edges was that. */
            while (jf != kf)
              { assert(jf != UNX_NULL); /* The pair {uxep[kf]} must be in the list. */
                assert(jf < mf);
                /* Move to the next pair: */
                jf_ant = jf;
                if (uxep[jf].e[0] == uxer)
                  { s_ant = 0; jf = fnextf[jf].c[0]; }
                else if (uxep[jf].e[1] == uxer)
                  { s_ant = 1; jf = fnextf[jf].c[1]; }
                else
                  { assert(FALSE); /* Edge {uxe} must be in every pair of its list. */ }
              }
            /* Now remove {uxep[kf].e[r]} from the list of edge {uxer}: */
            assert((jf_ant != UNX_NULL) && (s_ant != -1));
            fnextf[jf_ant].c[s_ant] = fnextf[kf].c[r];
          }
        fnextf[kf].c[r] = UNX_NULL;
        uxep[kf].e[r] = UNX_NULL;
      }
  }