#include <cstdio>
#include <vector>

struct vector3
{
  vector3(float x_, float y_, float z_) :
    x(x_), y(y_), z(z_) {}

  float x, y, z;
};

struct vector5
{
  vector5(float x_, float y_, float z_, float w_, float q_) :
    x(x_), y(y_), z(z_), w(w_), q(q_) {}

  float x, y, z, w, q;
};

int dims[3] = {4, 4, 4};
int volume[4 * 4 * 4];

int f(int i, int j, int k)
{
  return volume[i + dims[0] * (j + dims[1] * k)];
}

void
greedy(const int *volume, const int *dims, std::vector<vector3> &vertices, std::vector<vector5> &faces)
{
  for (std::size_t axis = 0; axis < 3; ++axis)
  {
    // printf("axis: %ld\n", axis);

    const std::size_t u = (axis + 1) % 3;
    const std::size_t v = (axis + 2) % 3;

    // printf("u: %ld, v: %ld\n", u, v);

    int x[3] = {0}, q[3] = {0}, mask[dims[u] * dims[v]];
    // printf("dims[u]: %d, dims[v]: %d\n", dims[u], dims[v]);

    // printf("x: %d, %d, %d\n", x[0], x[1], x[2]);
    // printf("q: %d, %d, %d\n", q[0], q[1], q[2]);

    // Compute mask
    q[axis] = 1;
    for (x[axis] = -1; x[axis] < dims[axis];)
    {
      // printf("x: %d, %d, %d\n", x[0], x[1], x[2]);
      // printf("q: %d, %d, %d\n", q[0], q[1], q[2]);

      std::size_t counter = 0;
      for (x[v] = 0; x[v] < dims[v]; ++x[v])
        for (x[u] = 0; x[u] < dims[u]; ++x[u], ++counter)
        {
          const int a = 0 <= x[axis] ? f(x[0], x[1], x[2]) : 0;
          const int b = x[axis] < dims[axis] - 1 ? f(x[0] + q[0],
                                                     x[1] + q[1],
                                                     x[2] + q[2]) : 0;
          const bool ba = static_cast<bool>(a);
          if (ba == static_cast<bool>(b))
            mask[counter] = 0;
          else if (ba)
            mask[counter] = a;
          else
            mask[counter] = -b;
        }

      ++x[axis];

      // for (size_t a = 0; a < dims[u] * dims[v]; ++a)
      //   printf("%d, ", mask[a]);
      // printf("\n");

      // Generate mesh for mask using lexicographic ordering
      std::size_t width = 0, height = 0;

      counter = 0;
      for (std::size_t j = 0; j < dims[v]; ++j)
        for (std::size_t i = 0; i < dims[u];)
        {
          int c = mask[counter];
          if (c)
          {
            // Compute width
            for (width = 1; c == mask[counter + width] &&
                   i + width < dims[u]; ++width)
            {}

            // Compute height
            bool done = false;
            for (height = 1; j + height < dims[v]; ++height)
            {
              for (std::size_t k = 0; k < width; ++k)
                if (c != mask[counter + k + height * dims[u]])
                {
                  done = true;
                  break;
                }

              if (done)
                break;
            }

            // Add quad
            x[u] = i;
            x[v] = j;

            int du[3] = {0}, dv[3] = {0};

            if (c > 0)
            {
              dv[v] = height;
              du[u] = width;
            }
            else
            {
              c = -c;
              du[v] = height;
              dv[u] = width;
            }

            const std::size_t vertexSize = vertices.size();
            vertices.push_back(vector3(x[0], x[1], x[2]));
            vertices.push_back(vector3(x[0] + du[0], x[1] + du[1], x[2] + du[2]));
            vertices.push_back(vector3(x[0] + du[0] + dv[0], x[1] + du[1] + dv[1], x[2] + du[2] + dv[2]));
            vertices.push_back(vector3(x[0] + dv[0], x[1] + dv[1], x[2] + dv[2]));

            faces.push_back(vector5(vertexSize, vertexSize + 1, vertexSize + 2, vertexSize + 3, c));

            for (std::size_t b = 0; b < width; ++b)
              for (std::size_t a = 0; a < height; ++a)
                mask[counter + b + a * dims[u]] = 0;

            // Increment counters
            i += width; counter += width;
          }
          else
          {
            ++i;
            ++counter;
          }
        }
    }

    // printf("x: %d, %d, %d\n", x[0], x[1], x[2]);
    // printf("q: %d, %d, %d\n", q[0], q[1], q[2]);


    // printf("\n");

  }
}

int
main(void)
{
  std::vector<vector3> vertices;
  std::vector<vector5> faces;

  for (size_t a = 0; a < dims[0] * dims[1] * dims[2]; ++a)
    volume[a] = 255;

  greedy(volume, dims, vertices, faces);

  for (size_t a = 0; a < vertices.size(); ++a)
    printf("%ld: %f, %f, %f\n", a, vertices[a].x, vertices[a].y, vertices[a].z);
  printf("\n");

  for (size_t a = 0; a < faces.size(); ++a)
    printf("%ld: %f, %f, %f, %f, %f\n", a, faces[a].x, faces[a].y, faces[a].z, faces[a].w, faces[a].q);


  return 0;
}