// Given a mesh that was created for 3D modelling, this code will simplify it into a mesh
// suitable for use in a wireframe display, i.e. with the spurious lines removed that slash
// across flat planes planes. (Those lines are caused by the triangulation of quads or
// larger polygons from the original object's surface which is done when the 3D model is
// converted to a mesh.  A mesh in 3D modelling is almost always triangle-based because a
// triangle *must* be planar but a quad or polygon does not automatically have to be...)

// Note that this code can be used in conjunction with https://gtoal.com/vectrex/tailgunner/experiments/edges2faces.c
// which determines a list of faces from a list of only edges.

// Also helpful may be ~/github/uparse-main/lang/obj3d/obj3d.g on my home machine which
// is a parser for .obj files that outputs a C header file containg the vertices and edges
// of the model. (Needs my 'uparse' parser)  This converter will eventually be uploaded
// to github.

// see images  maxresdefault.jpg and triangulated_polygon.png in this directory for
// examples of why we might want to do this (images from https://www.youtube.com/watch?v=fu5HYNu-lw4 )

// code for converting to and from spherical coordinates (and incidentally projecting 3D coords
// to 2D screen coords using lookup tables rather than arithmetic) is in ~/src/3d/spherical

/*
   To convert a triangular mesh to remove the lines that cut through coplanar adjacent surfaces:

   * pick any triangle.  For each edge:
     * find the 3rd point of an adjacent edge-sharing triangle, if present.
     * Test if that point is coplanar with the original triangle.
     * If it is, remove the shared edge and extend the original triangle to a quad, including the new point
   * Recurse on the two added edges of the added triangle.  Except the surface that is being extended is now a polygon.
     * If the added edge is already part of a polygon rather than a triangle, pick any other point in that polygon that
       is not colinear with the common edge, and if it is also coplanar with the original triangle, merge the two polygons
       by removing the common edge and joining the two polygons at the endpoints of that removed edge.
   * We do not yet check for a planar donut shape where the depth-first search starting at the first triangle
     ends up reconnecting with that triangle after going around a loop.  Let's assume for now there are none
     in the models we are interested in. 
 */


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

#define EPSILON 0.001

static int debug = 0;

// determine equation of a plane from 3 points in space
void get_plane_equation(const float x1, const float y1, const float z1,
                        const float x2, const float y2, const float z2,
                        const float x3, const float y3, const float z3,
                        float *ax, float *by, float *cz, float *d) {
  float a1 = x2 - x1, a2 = x3 - x1, b1 = y2 - y1, b2 = y3 - y1, c1 = z2 - z1, c2 = z3 - z1;
  *ax = b1*c2 - b2*c1;
  *by = a2*c1 - a1*c2;
  *cz = a1*b2 - b1*a2;
  *d = -(*ax*x1 + *by*y1 + *cz*z1);
}

// is tested point coplanar with plane described by given equation of a plane?
int on_plane(const float A, const float B, const float C, const float D,
             const float x, const float y, const float z,
             const float epsilon) {
  float result = A*x + B*y + C*z + D;
  return -epsilon < result && result < epsilon;
}

/*
  Our test cube, and axes.


                                                 Z
                                                 |
                                                 |
                                                 |
                                                 |
                                                 |             Y
                                                 |            /
                                                 |           /
                                      4          |          /      1
                                       *---------|---------------*
                                      /:         |        /     /|
                                     / :         |       /     / |
                                    /  :         |      /     /  |
                                   /   :         |     /     /   |
                                  /    :         |    /     /    |
                                 /     :         |   /     /     |
                             3  /      :         |  /     /2     |
                               *-------------------------*       |
                               |       :         |/      |       |
                               |       :         +-------|--------------------------- X
                               |       :        0        |       |
                               |     8 *.................|.......* 5
                               |      /                  |      /
                               |     /                   |     /
                               |    /                    |    /
                               |   /                     |   /
                               |  /                      |  /
                               | /                       | /
                               |/                        |/
                               *-------------------------*
                              7                            6

 */
const float vertex[8+1][3] = {
  { 0.0,  0.0,  0.0}, // unused vertex, because a) convention is vertices are numbered from 1;
                      // and b) I'm using vertex[0] to store coord of rotational center of model.
  { 1.0,  1.0,  1.0}, // Vertex 1
  { 1.0, -1.0,  1.0}, // Vertex 2
  {-1.0, -1.0,  1.0}, // Vertex 3
  {-1.0,  1.0,  1.0}, // Vertex 4
  { 1.0,  1.0, -1.0}, // Vertex 5
  { 1.0, -1.0, -1.0}, // Vertex 6
  {-1.0, -1.0, -1.0}, // Vertex 7
  {-1.0,  1.0, -1.0}  // Vertex 8
};

// test whether three vertices of original seed triangle plus vertex that we're considering adding are in the same plane:
int coplanar(int v1, int v2, int v3, int test) {
  float A, B, C, D,
        x1 = vertex[v1][0],
        y1 = vertex[v1][1],
        z1 = vertex[v1][2],
        x2 = vertex[v2][0],
        y2 = vertex[v2][1],
        z2 = vertex[v2][2],
        x3 = vertex[v3][0],
        y3 = vertex[v3][1],
        z3 = vertex[v3][2];
  get_plane_equation(x1,y1,z1, x2,y2,z2, x3,y3,z3, &A, &B, &C, &D);
  return on_plane(A, B, C, D,
                  vertex[test][0] /*x*/,
                  vertex[test][1] /*y*/,
                  vertex[test][2] /*z*/,
                  EPSILON);
}

const int triangles[6*2][3] = {
  {1, 2, 3}, {3, 4, 1},  // Top face
  {5, 6, 7}, {7, 8, 5},  // Bottom face
  {3, 4, 8}, {8, 7, 3},  // Left face
  {1, 5, 6}, {6, 2, 1},  // Right face
  {3, 2, 6}, {6, 7, 3},  // Front face
  {4, 1, 5}, {5, 8, 4}   // Back face
};
int merged[12]               = {0,0,0,0,0,0,0,0,0,0,0,0};

int polygon_vertices[12][12];
int polygon_vertex_count[12] = {0,0,0,0,0,0,0,0,0,0,0,0};
int next_poly;

int try(int poly, int vertex_tri, int vertex_a, int vertex_b); // some cross-recursion going on here!

// given that we now know an added triangle is part of the same plane as the existing triangle or polygon,
// add it to the existing structure and remove it from the pool (after recursing on any triangles connected
// to the other two sides of the new triangle).
int merge_tri(int poly, int vertex_tri, int new_tri, int TA, int a, int TB, int b) {
  int TC, d;
  if (debug) fprintf(stderr, "Match found.  Merging triangles %d and %d with common edge between vertices %d==%d and %d==%d\n", vertex_tri, new_tri, TA, a, TB, b);
  // edges a - b on both triangles are shared!
  // a, b, c may not be the original order of vertices but they've now been shuffled so that the nomenclature works!

  /*

                      triangles (or polys)       quads (or polys)
              
                              d                         d
                              *                         *
                             / \                       / \
                            /   \                     /   \
                           /     \                   /     \
                          /       \                 /       \
     common edge -->   a *---------* b           a *         * b   <-- edge is eliminated
                          \       /                 \       /
                           \     /                   \     /
                            \   /                     \   /
                             \ /                       \ /
                              *                         *
                              c                         c


       given inputs a - c - b ( - a) and a - d - b ( - a ) , the joined polygon will be  a - c - b - d ( - a )
       even if c or d are multiple nodes (ie the inputs are polygons, not triangles).

   */
  if (!merged[vertex_tri]) { // actually I think this has to record a combination of triangle and edge
                             // since the triangle can be merged with up to 3 others.  This simplified
                             // test only works when the original input was in the form of quads.
    
    // initialise polygon with first triangle
    if ((triangles[vertex_tri][0] == TA || triangles[vertex_tri][1] == TA) && (triangles[vertex_tri][0] == TB || triangles[vertex_tri][1] == TB)) TC = triangles[vertex_tri][2];
    if ((triangles[vertex_tri][0] == TA || triangles[vertex_tri][2] == TA) && (triangles[vertex_tri][0] == TB || triangles[vertex_tri][2] == TB)) TC = triangles[vertex_tri][1];
    if ((triangles[vertex_tri][1] == TA || triangles[vertex_tri][2] == TA) && (triangles[vertex_tri][1] == TB || triangles[vertex_tri][2] == TB)) TC = triangles[vertex_tri][0];
    polygon_vertices[poly][polygon_vertex_count[poly]++] = TA;
    polygon_vertices[poly][polygon_vertex_count[poly]++] = TC;
    polygon_vertices[poly][polygon_vertex_count[poly]++] = TB;
    merged[vertex_tri] = 1; // this might need to be completed after the recursion.  Maybe make it tri-state rather than boolean: (unexamined, in_process, merged)
  }

  if ((triangles[new_tri][0] == a || triangles[new_tri][1] == a) && (triangles[new_tri][0] == b || triangles[new_tri][1] == b)) d = triangles[new_tri][2];
  if ((triangles[new_tri][0] == a || triangles[new_tri][2] == a) && (triangles[new_tri][0] == b || triangles[new_tri][2] == b)) d = triangles[new_tri][1];
  if ((triangles[new_tri][1] == a || triangles[new_tri][2] == a) && (triangles[new_tri][1] == b || triangles[new_tri][2] == b)) d = triangles[new_tri][0];
  polygon_vertices[poly][polygon_vertex_count[poly]++] = d; // (if it was a poly being added there will be more lines like this one...)
    // additional triangle/polygon to be merged.
    // if it's a polygon, the order is important! Fortunately only adding triangles for now.
  
  // recurse on the two unexamined edges of the added triangle. a - b has been done.
  merged[new_tri] = 1;
  if (!try(poly, new_tri, a, d)) try(poly, new_tri, d, a);
  if (!try(poly, new_tri, b, d)) try(poly, new_tri, d, b);
  return 1;
}

int common_edge(int v1, int a, int v2, int b, int v3, int c) {
  if (debug) fprintf(stderr, "Compare v%d == v%d && v%d == v%d (%s)\n", v1, a, v2, b, (v1==a) && (v2==b) ? "true" : "false");
  if ((v1==a) && (v2==b)) {
    return coplanar(v1, v2, v3, c); // might as well include that test here and simplify the calling code.
  } else return 0;
}

int try(int poly, int vertex_tri, int vertex_a, int vertex_b) {
  int found = 0;
  int vertex_c;

  if (debug) fprintf(stderr, "trying triangle %d edge between vertices %d and %d\n", vertex_tri, vertex_a, vertex_b);
  if (merged[vertex_tri]) return 0; // This triangle has already been added to a polygon.
  if (debug) fprintf(stderr, "triangle %d not merged already\n", vertex_tri);

  vertex_c = 3-vertex_a-vertex_b; // 0 + 1 + 2 = 3 so subtracting any two vertex numbers from 3 gives the index of the third vertex...

  // let's look at all other triangles that haven't yet been assimilated...
  for (int new_tri = vertex_tri + 1; new_tri < 12; new_tri++) {
    if (debug) fprintf(stderr, "comparing against triangle %d\n", new_tri);
    if (merged[new_tri]) continue;
    if (debug) fprintf(stderr, "triangle %d not merged already\n", new_tri);
    // Is the given edge present on any other triangle?
    int a, b, c, TA, TB, TC;
    a = triangles[new_tri][0];
    b = triangles[new_tri][1];
    c = triangles[new_tri][2];
    TA = triangles[vertex_tri][vertex_a];
    TB = triangles[vertex_tri][vertex_b];
    TC = triangles[vertex_tri][vertex_c];
    // 6x6 == 36 comparisons
    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// a b
         if (common_edge(TA, a, TB, b, TC, c)) found |= merge_tri(poly, vertex_tri, new_tri, TA, a, TB, b);
    else if (common_edge(TA, a, TB, c, TC, b)) found |= merge_tri(poly, vertex_tri, new_tri, TA, a, TB, c);
    else if (common_edge(TA, b, TB, a, TC, c)) found |= merge_tri(poly, vertex_tri, new_tri, TA, b, TB, a);
    else if (common_edge(TA, b, TB, c, TC, a)) found |= merge_tri(poly, vertex_tri, new_tri, TA, b, TB, c);
    else if (common_edge(TA, c, TB, a, TC, b)) found |= merge_tri(poly, vertex_tri, new_tri, TA, c, TB, a);
    else if (common_edge(TA, c, TB, b, TC, a)) found |= merge_tri(poly, vertex_tri, new_tri, TA, c, TB, b);
    else if (common_edge(TB, a, TA, b, TC, c)) found |= merge_tri(poly, vertex_tri, new_tri, TB, a, TA, b);
    else if (common_edge(TB, a, TA, c, TC, b)) found |= merge_tri(poly, vertex_tri, new_tri, TB, a, TA, c);
    else if (common_edge(TB, b, TA, a, TC, c)) found |= merge_tri(poly, vertex_tri, new_tri, TB, b, TA, a);
    else if (common_edge(TB, b, TA, c, TC, a)) found |= merge_tri(poly, vertex_tri, new_tri, TB, b, TA, c);
    else if (common_edge(TB, c, TA, a, TC, b)) found |= merge_tri(poly, vertex_tri, new_tri, TB, c, TA, a);
    else if (common_edge(TB, c, TA, b, TC, a)) found |= merge_tri(poly, vertex_tri, new_tri, TB, c, TA, b);
    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// a c
         if (common_edge(TA, a, TC, b, TB, c)) found |= merge_tri(poly, vertex_tri, new_tri, TA, a, TC, b);
    else if (common_edge(TA, a, TC, c, TB, b)) found |= merge_tri(poly, vertex_tri, new_tri, TA, a, TC, c);
    else if (common_edge(TA, b, TC, a, TB, c)) found |= merge_tri(poly, vertex_tri, new_tri, TA, b, TC, a);
    else if (common_edge(TA, b, TC, c, TB, a)) found |= merge_tri(poly, vertex_tri, new_tri, TA, b, TC, c);
    else if (common_edge(TA, c, TC, a, TB, b)) found |= merge_tri(poly, vertex_tri, new_tri, TA, c, TC, a);
    else if (common_edge(TA, c, TC, b, TB, a)) found |= merge_tri(poly, vertex_tri, new_tri, TA, c, TC, b);
    else if (common_edge(TC, a, TA, b, TB, c)) found |= merge_tri(poly, vertex_tri, new_tri, TC, a, TA, b);
    else if (common_edge(TC, a, TA, c, TB, b)) found |= merge_tri(poly, vertex_tri, new_tri, TC, a, TA, c);
    else if (common_edge(TC, b, TA, a, TB, c)) found |= merge_tri(poly, vertex_tri, new_tri, TC, b, TA, a);
    else if (common_edge(TC, b, TA, c, TB, a)) found |= merge_tri(poly, vertex_tri, new_tri, TC, b, TA, c);
    else if (common_edge(TC, c, TA, a, TB, b)) found |= merge_tri(poly, vertex_tri, new_tri, TC, c, TA, a);
    else if (common_edge(TC, c, TA, b, TB, a)) found |= merge_tri(poly, vertex_tri, new_tri, TC, c, TA, b);
    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// b c
         if (common_edge(TB, a, TC, b, TA, c)) found |= merge_tri(poly, vertex_tri, new_tri, TB, a, TC, b);
    else if (common_edge(TB, a, TC, c, TA, b)) found |= merge_tri(poly, vertex_tri, new_tri, TB, a, TC, c);
    else if (common_edge(TB, b, TC, a, TA, c)) found |= merge_tri(poly, vertex_tri, new_tri, TB, b, TC, a);
    else if (common_edge(TB, b, TC, c, TA, a)) found |= merge_tri(poly, vertex_tri, new_tri, TB, b, TC, c);
    else if (common_edge(TB, c, TC, a, TA, b)) found |= merge_tri(poly, vertex_tri, new_tri, TB, c, TC, a);
    else if (common_edge(TB, c, TC, b, TA, a)) found |= merge_tri(poly, vertex_tri, new_tri, TB, c, TC, b);
    else if (common_edge(TC, a, TB, b, TA, c)) found |= merge_tri(poly, vertex_tri, new_tri, TC, a, TB, b);
    else if (common_edge(TC, a, TB, c, TA, b)) found |= merge_tri(poly, vertex_tri, new_tri, TC, a, TB, c);
    else if (common_edge(TC, b, TB, a, TA, c)) found |= merge_tri(poly, vertex_tri, new_tri, TC, b, TB, a);
    else if (common_edge(TC, b, TB, c, TA, a)) found |= merge_tri(poly, vertex_tri, new_tri, TC, b, TB, c);
    else if (common_edge(TC, c, TB, a, TA, b)) found |= merge_tri(poly, vertex_tri, new_tri, TC, c, TB, a);
    else if (common_edge(TC, c, TB, b, TA, a)) found |= merge_tri(poly, vertex_tri, new_tri, TC, c, TB, b);
  }
  return found; // if any were found.  may be up to 3 for a given triangle plus any found recursively.
}

int main(int argc, char **argv) {
  int found;
  int tri, vertex, poly;

  next_poly = 0;
  for (tri = 0; tri < 12; tri++) {
    found = 0;
    if (debug) fprintf(stderr, "Examine base triangle %d\n", tri);
    if (!merged[tri]) found |= try(next_poly, tri, 0, 1); // An edge might not be shared with any other triangles in the same
    if (!merged[tri]) found |= try(next_poly, tri, 1, 2); // plane, so we have to check all three edges until we find one that
    if (!merged[tri]) found |= try(next_poly, tri, 2, 0); // is coplanar with another triangle. Some thought will be needed if
    if (!merged[tri]) found |= try(next_poly, tri, 1, 0); // 2 or 3 adjacent triangles are coplanar in a tessellated surface!
    if (!merged[tri]) found |= try(next_poly, tri, 2, 1);
    if (!merged[tri]) found |= try(next_poly, tri, 0, 2);

    if (polygon_vertex_count[next_poly] > 0) { // i.e. found == 1
      if (debug) fprintf(stderr, "Base triangle %d - edges were added to polygon[%d]\n", tri, next_poly);
      next_poly += 1; // we've fully explored the data and added every triangle to this polygon that we can...
    }
  }

  // Output resulting polygons
  printf("Simplified mesh polygons:\n");
  for (poly = 0; poly < next_poly; poly++) {
    printf("Polygon %d: ", poly);
    for (vertex = 0; vertex < polygon_vertex_count[poly]; vertex++) printf("%d ", polygon_vertices[poly][vertex]);
    printf("%d\n", polygon_vertices[poly][0]);
  }

  exit(EXIT_SUCCESS);
  return EXIT_FAILURE;
}
