surface_mesh.cpp

// Copyright 2017-2023, Nicholas Sharp and the Polyscope contributors. https://polyscope.run

#include "polyscope/surface_mesh.h"

#include "glm/fwd.hpp"
#include "polyscope/combining_hash_functions.h"
#include "polyscope/pick.h"
#include "polyscope/polyscope.h"
#include "polyscope/render/engine.h"

#include "imgui.h"
#include "polyscope/types.h"
#include "polyscope/utilities.h"

#include <unordered_map>
#include <utility>

namespace polyscope {

// Initialize statics
const std::string SurfaceMesh::structureTypeName = "Surface Mesh";


SurfaceMesh::SurfaceMesh(std::string name_)
    : QuantityStructure<SurfaceMesh>(name_, typeName()),
      // clang-format off

// == managed quantities

// positions
vertexPositions(           this, uniquePrefix() + "vertexPositions",     vertexPositionsData),

// connectivity / indices
// (triangle and face inds are always computed initially when we triangulate the mesh)
triangleVertexInds(        this, uniquePrefix() + "triangleVertexInds",          triangleVertexIndsData),
triangleFaceInds(          this, uniquePrefix() + "triangleFaceInds",            triangleFaceIndsData),
triangleCornerInds(        this, uniquePrefix() + "triangleCornerInds",          triangleCornerIndsData,         std::bind(&SurfaceMesh::computeTriangleCornerInds, this)),
triangleAllVertexInds(     this, uniquePrefix() + "triangleAllVertexInds",       triangleAllVertexIndsData,      std::bind(&SurfaceMesh::computeTriangleAllVertexInds, this)),
triangleAllEdgeInds(       this, uniquePrefix() + "triangleAllEdgeInds",         triangleAllEdgeIndsData,        std::bind(&SurfaceMesh::computeTriangleAllEdgeInds, this)),
triangleAllHalfedgeInds(   this, uniquePrefix() + "triangleHalfedgeInds",     triangleAllHalfedgeIndsData,    std::bind(&SurfaceMesh::computeTriangleAllHalfedgeInds, this)),
triangleAllCornerInds(     this, uniquePrefix() + "triangleAllCornerInds",    triangleAllCornerIndsData,      std::bind(&SurfaceMesh::computeTriangleAllCornerInds, this)),

// internal triangle data for rendering
baryCoord(              this, uniquePrefix() + "baryCoord",           baryCoordData),
edgeIsReal(             this, uniquePrefix() + "edgeIsReal",          edgeIsRealData),

// other internally-computed geometry
faceNormals(            this, uniquePrefix() + "faceNormals",         faceNormalsData,        std::bind(&SurfaceMesh::computeFaceNormals, this)),
faceCenters(            this, uniquePrefix() + "faceCenters",         faceCentersData,        std::bind(&SurfaceMesh::computeFaceCenters, this)),         
faceAreas(              this, uniquePrefix() + "faceAreas",           faceAreasData,          std::bind(&SurfaceMesh::computeFaceAreas, this)),
vertexNormals(          this, uniquePrefix() + "vertexNormals",       vertexNormalsData,      std::bind(&SurfaceMesh::computeVertexNormals, this)),
vertexAreas(            this, uniquePrefix() + "vertexAreas",         vertexAreasData,        std::bind(&SurfaceMesh::computeVertexAreas, this)),

// tangent spaces
defaultFaceTangentBasisX(   this, uniquePrefix() + "defaultFaceTangentBasisX",  defaultFaceTangentBasisXData,  std::bind(&SurfaceMesh::computeDefaultFaceTangentBasisX, this)),
defaultFaceTangentBasisY(   this, uniquePrefix() + "defaultFaceTangentBasisY",  defaultFaceTangentBasisYData,  std::bind(&SurfaceMesh::computeDefaultFaceTangentBasisY, this)),

// == persistent options
surfaceColor(           uniquePrefix() + "surfaceColor",    getNextUniqueColor()),
edgeColor(              uniquePrefix() + "edgeColor",       glm::vec3{0., 0., 0.}), material(uniquePrefix() + "material", "clay"),
edgeWidth(              uniquePrefix() + "edgeWidth",       0.),
backFacePolicy(         uniquePrefix() + "backFacePolicy",  BackFacePolicy::Different),
backFaceColor(          uniquePrefix() + "backFaceColor",   glm::vec3(1.f - surfaceColor.get().r, 1.f - surfaceColor.get().g, 1.f - surfaceColor.get().b)),
shadeStyle(             uniquePrefix() + "shadeStyle",      MeshShadeStyle::Flat)

// clang-format on
{}

SurfaceMesh::SurfaceMesh(std::string name_, const std::vector<glm::vec3>& vertexPositions_,
                         const std::vector<uint32_t>& faceIndsEntries_, const std::vector<uint32_t>& faceIndsStart_)
    : SurfaceMesh(name_) {

  vertexPositionsData = vertexPositions_;
  faceIndsEntries = faceIndsEntries_;
  faceIndsStart = faceIndsStart_;

  vertexPositions.checkInvalidValues();
  computeConnectivityData();
  updateObjectSpaceBounds();
}

SurfaceMesh::SurfaceMesh(std::string name_, const std::vector<glm::vec3>& vertexPositions_,
                         const std::vector<std::vector<size_t>>& facesIn)
    : SurfaceMesh(name_) {

  vertexPositionsData = vertexPositions_;
  nestedFacesToFlat(facesIn);

  vertexPositions.checkInvalidValues();
  computeConnectivityData();
  updateObjectSpaceBounds();
}

void SurfaceMesh::nestedFacesToFlat(const std::vector<std::vector<size_t>>& nestedInds) {

  faceIndsStart.clear();
  faceIndsEntries.clear();
  faceIndsStart.push_back(0);

  for (const std::vector<size_t>& face : nestedInds) {
    for (size_t iV : face) {
      faceIndsEntries.push_back(iV);
    }
    faceIndsStart.push_back(faceIndsEntries.size());
  }
}

void SurfaceMesh::computeConnectivityData() {

  // some number-of-elements arithmetic
  size_t numFaces = faceIndsStart.size() - 1;
  nCornersCount = faceIndsEntries.size();
  nFacesTriangulationCount = nCornersCount - 2 * numFaces;

  // fill out these buffers as we construct the triangulation
  triangleVertexIndsData.clear();
  triangleVertexIndsData.resize(3 * nFacesTriangulationCount);
  triangleFaceIndsData.clear();
  triangleFaceIndsData.resize(3 * nFacesTriangulationCount);
  baryCoordData.clear();
  baryCoordData.resize(3 * nFacesTriangulationCount);
  edgeIsRealData.clear();
  edgeIsRealData.resize(3 * nFacesTriangulationCount);

  // validate the face-vertex indices
  for (size_t iV : faceIndsEntries) {
    if (iV >= vertexPositions.size())
      exception("SurfaceMesh " + name + " has face vertex index " + std::to_string(iV) +
                " out of bounds for number of vertices " + std::to_string(vertexPositions.size()));
  }

  // construct the triangualted draw list and all other related data
  size_t iTriFace = 0;
  for (size_t iF = 0; iF < numFaces; iF++) {
    size_t D = faceIndsStart[iF + 1] - faceIndsStart[iF];

    size_t iStart = faceIndsStart[iF];
    uint32_t vRoot = faceIndsEntries[iStart];

    // implicitly triangulate from root
    for (size_t j = 1; (j + 1) < D; j++) {
      uint32_t vB = faceIndsEntries[iStart + j];
      uint32_t vC = faceIndsEntries[iStart + ((j + 1) % D)];

      // triangle vertex indices
      triangleVertexIndsData[3 * iTriFace + 0] = vRoot;
      triangleVertexIndsData[3 * iTriFace + 1] = vB;
      triangleVertexIndsData[3 * iTriFace + 2] = vC;

      // triangle face indices
      for (size_t k = 0; k < 3; k++) triangleFaceIndsData[3 * iTriFace + k] = iF;

      // barycentric coordinates
      baryCoordData[3 * iTriFace + 0] = glm::vec3{1., 0., 0.};
      baryCoordData[3 * iTriFace + 1] = glm::vec3{0., 1., 0.};
      baryCoordData[3 * iTriFace + 2] = glm::vec3{0., 0., 1.};

      // internal edges for triangulated polygons
      glm::vec3 edgeRealV{0., 1., 0.};
      if (j == 1) {
        edgeRealV.x = 1.;
      }
      if (j + 2 == D) {
        edgeRealV.z = 1.;
      }
      for (size_t k = 0; k < 3; k++) edgeIsRealData[3 * iTriFace + k] = edgeRealV;

      iTriFace++;
    }
  }

  vertexDataSize = nVertices();
  faceDataSize = nFaces();
  // edgeDataSize = ... we don't know this yet, gets set below
  halfedgeDataSize = nHalfedges();
  cornerDataSize = nCorners();

  triangleVertexInds.markHostBufferUpdated();
  triangleFaceInds.markHostBufferUpdated();
  baryCoord.markHostBufferUpdated();
  edgeIsReal.markHostBufferUpdated();
}

// =================================================
// =====    Lazily-Populated Connectivity   ========
// =================================================

void SurfaceMesh::computeTriangleAllEdgeInds() {

  // WARNING: logic duplicated in countEdges()

  if (edgePerm.empty())
    exception("SurfaceMesh " + name +
              " performed an operation which requires edge indices to be specified, but none have been set. "
              "Call setEdgePermutation().");

  triangleVertexInds.ensureHostBufferPopulated();
  triangleAllEdgeInds.data.resize(3 * 3 * nFacesTriangulation());
  halfedgeEdgeCorrespondence.resize(nHalfedges());

  // used to loop over edges
  std::unordered_map<std::pair<size_t, size_t>, size_t, polyscope::hash_combine::hash<std::pair<size_t, size_t>>>
      seenEdgeInds;

  auto createEdgeKey = [&](size_t a, size_t b) -> std::pair<size_t, size_t> {
    return std::make_pair(std::min(a, b), std::max(a, b));
  };

  size_t psEdgeInd = 0; // polyscope's edge index, iterated according to Polyscope's canonical ordering
  for (size_t iF = 0; iF < nFaces(); iF++) {
    size_t start = faceIndsStart[iF];
    size_t D = faceIndsStart[iF + 1] - start;

    // TODO why can't we use edges on non triangular meshes? Implement it.

    if (D != 3) {
      exception("SurfaceMesh " + name +
                " attempted to access triangle-edge indices, but it has non-triangular faces. These indices are "
                "only well-defined on a pure-triangular mesh.");
    }

    glm::uvec3 thisTriInds{0, 0, 0};
    for (size_t j = 0; j < 3; j++) {
      size_t vA = triangleVertexInds.data[3 * iF + j];
      size_t vB = triangleVertexInds.data[3 * iF + ((j + 1) % 3)];

      std::pair<size_t, size_t> key = createEdgeKey(vA, vB);

      size_t thisEdgeInd;
      if (seenEdgeInds.find(key) == seenEdgeInds.end()) {
        // process a new edge in the canonical ordering
        if (psEdgeInd >= edgePerm.size()) {
          exception("SurfaceMesh " + name +
                    " edge indexing out of bounds. Did you pass an edge ordering that is too short?");
        }
        thisEdgeInd = edgePerm[psEdgeInd];
        seenEdgeInds[key] = thisEdgeInd;
        psEdgeInd++;
      } else {
        // we've processed this edge before, retrieve its assigned index
        thisEdgeInd = seenEdgeInds[key];
      }

      halfedgeEdgeCorrespondence[start + j] = thisEdgeInd;
      thisTriInds[j] = thisEdgeInd;
    }

    for (size_t j = 0; j < 3; j++) {
      for (size_t k = 0; k < 3; k++) {
        triangleAllEdgeInds.data[9 * iF + 3 * j + k] = thisTriInds[k];
      }
    }
  }

  nEdgesCount = psEdgeInd;
  triangleAllEdgeInds.markHostBufferUpdated();
}

void SurfaceMesh::countEdges() {

  // WARNING: logic duplicated in computeTriangleAllEdgeInds()

  // used to loop over edges
  std::unordered_map<std::pair<size_t, size_t>, size_t, polyscope::hash_combine::hash<std::pair<size_t, size_t>>>
      seenEdgeInds;

  auto createEdgeKey = [&](size_t a, size_t b) -> std::pair<size_t, size_t> {
    return std::make_pair(std::min(a, b), std::max(a, b));
  };

  size_t psEdgeInd = 0; // polyscope's edge index, iterated according to Polyscope's canonical ordering
  for (size_t iF = 0; iF < nFaces(); iF++) {
    size_t start = faceIndsStart[iF];
    size_t D = faceIndsStart[iF + 1] - start;

    if (D != 3) {
      exception("SurfaceMesh " + name +
                " attempted to count edges, but mesh has non-triangular faces. Edge functions are only implemented on "
                "a pure-triangular mesh.");
    }

    for (size_t j = 0; j < 3; j++) {
      size_t vA = triangleVertexInds.data[3 * iF + j];
      size_t vB = triangleVertexInds.data[3 * iF + ((j + 1) % 3)];

      std::pair<size_t, size_t> key = createEdgeKey(vA, vB);

      if (seenEdgeInds.find(key) == seenEdgeInds.end()) {
        size_t thisEdgeInd = psEdgeInd;
        seenEdgeInds[key] = thisEdgeInd;
        psEdgeInd++;
      }
    }
  }

  nEdgesCount = psEdgeInd;
}

size_t SurfaceMesh::nEdges() {
  if (nEdgesCount == INVALID_IND) countEdges();
  return nEdgesCount;
}

void SurfaceMesh::computeTriangleCornerInds() {

  triangleCornerInds.data.clear();
  triangleCornerInds.data.reserve(3 * nFacesTriangulation());

  for (size_t iF = 0; iF < nFaces(); iF++) {
    size_t iStart = faceIndsStart[iF];
    size_t D = faceIndsStart[iF + 1] - iStart;

    // emit the data for triangles triangulating this face
    for (size_t j = 1; (j + 1) < D; j++) {
      uint32_t c0 = iStart;
      uint32_t c1 = iStart + j;
      uint32_t c2 = iStart + j + 1;

      triangleCornerInds.data.push_back(c0);
      triangleCornerInds.data.push_back(c1);
      triangleCornerInds.data.push_back(c2);
    }
  }

  triangleCornerInds.markHostBufferUpdated();
}

void SurfaceMesh::computeTriangleAllVertexInds() {

  triangleAllVertexInds.data.clear();
  triangleAllVertexInds.data.reserve(3 * 3 * nFacesTriangulation());

  for (size_t iF = 0; iF < nFaces(); iF++) {
    size_t iStart = faceIndsStart[iF];
    size_t D = faceIndsStart[iF + 1] - iStart;
    uint32_t vRoot = faceIndsEntries[iStart];

    // implicitly triangulate from root
    for (size_t j = 1; (j + 1) < D; j++) {
      uint32_t vB = faceIndsEntries[iStart + j];
      uint32_t vC = faceIndsEntries[iStart + ((j + 1) % D)];

      // triangle vertex indices, all three values-each
      for (size_t k = 0; k < 3; k++) {
        triangleAllVertexInds.data.push_back(vRoot);
        triangleAllVertexInds.data.push_back(vB);
        triangleAllVertexInds.data.push_back(vC);
      }
    }
  }

  triangleAllVertexInds.markHostBufferUpdated();
}

void SurfaceMesh::computeTriangleAllHalfedgeInds() {

  triangleAllHalfedgeInds.data.clear();
  triangleAllHalfedgeInds.data.reserve(3 * 3 * nFacesTriangulation());

  bool haveCustomIndex = !halfedgePerm.empty();

  for (size_t iF = 0; iF < nFaces(); iF++) {
    size_t iStart = faceIndsStart[iF];
    size_t D = faceIndsStart[iF + 1] - iStart;

    // emit the data for triangles triangulating this face
    for (size_t j = 1; (j + 1) < D; j++) {

      // FORNOW: for polygonal faces, substitute the opposite-edge value for all internal edges of the triangulation

      uint32_t he0 = iStart + j; // this is a dummy value due to triangulation of polygons
      uint32_t he1 = iStart + j; // this is the actual right value for the opposite edge
      uint32_t he2 = iStart + j; // this is a dummy value due to triangulation of polygons

      // substitute non-dummy values for first and last edge if this is not an internal tri
      if (j == 1) he0 = iStart;
      if (j + 2 == D) he2 = iStart + D - 1;

      if (haveCustomIndex) {
        he0 = halfedgePerm[he0];
        he1 = halfedgePerm[he1];
        he2 = halfedgePerm[he2];
      }

      for (size_t k = 0; k < 3; k++) {
        triangleAllHalfedgeInds.data.push_back(he0);
        triangleAllHalfedgeInds.data.push_back(he1);
        triangleAllHalfedgeInds.data.push_back(he2);
      }
    }
  }

  triangleAllHalfedgeInds.markHostBufferUpdated();
}

void SurfaceMesh::computeTriangleAllCornerInds() {

  triangleAllCornerInds.data.clear();
  triangleAllCornerInds.data.reserve(3 * nFacesTriangulation());

  bool haveCustomIndex = !cornerPerm.empty();

  for (size_t iF = 0; iF < nFaces(); iF++) {
    size_t iStart = faceIndsStart[iF];
    size_t D = faceIndsStart[iF + 1] - iStart;

    // emit the data for triangles triangulating this face
    for (size_t j = 1; (j + 1) < D; j++) {
      uint32_t c0 = iStart;
      uint32_t c1 = iStart + j;
      uint32_t c2 = iStart + j + 1;

      if (haveCustomIndex) {
        c0 = cornerPerm[c0];
        c1 = cornerPerm[c1];
        c2 = cornerPerm[c2];
      }

      for (size_t k = 0; k < 3; k++) {
        triangleAllCornerInds.data.push_back(c0);
        triangleAllCornerInds.data.push_back(c1);
        triangleAllCornerInds.data.push_back(c2);
      }
    }
  }

  triangleAllCornerInds.markHostBufferUpdated();
}


// =================================================
// ========    Geometric Quantities      ==========
// =================================================


size_t SurfaceMesh::nVertices() { return vertexPositions.size(); }

void SurfaceMesh::computeFaceNormals() {

  vertexPositions.ensureHostBufferPopulated();

  faceNormals.data.resize(nFaces());

  for (size_t iF = 0; iF < nFaces(); iF++) {
    size_t iStart = faceIndsStart[iF];
    size_t D = faceIndsStart[iF + 1] - iStart;

    glm::vec3 fN{0., 0., 0.};
    if (D == 3) {
      glm::vec3 pA = vertexPositions.data[faceIndsEntries[iStart + 0]];
      glm::vec3 pB = vertexPositions.data[faceIndsEntries[iStart + 1]];
      glm::vec3 pC = vertexPositions.data[faceIndsEntries[iStart + 2]];
      fN = glm::cross(pB - pA, pC - pA);
    } else {
      for (size_t j = 0; j < D; j++) {
        glm::vec3 pA = vertexPositions.data[faceIndsEntries[iStart + j]];
        glm::vec3 pB = vertexPositions.data[faceIndsEntries[iStart + (j + 1) % D]];
        glm::vec3 pC = vertexPositions.data[faceIndsEntries[iStart + (j + 2) % D]];
        fN += glm::cross(pC - pB, pA - pB);
      }
    }
    fN = glm::normalize(fN);
    faceNormals.data[iF] = fN;
  }

  faceNormals.markHostBufferUpdated();
}

void SurfaceMesh::computeFaceCenters() {

  vertexPositions.ensureHostBufferPopulated();

  faceCenters.data.resize(nFaces());

  for (size_t iF = 0; iF < nFaces(); iF++) {
    size_t start = faceIndsStart[iF];
    size_t D = faceIndsStart[iF + 1] - start;
    glm::vec3 faceCenter{0., 0., 0.};
    for (size_t j = 0; j < D; j++) {
      glm::vec3 pA = vertexPositions.data[faceIndsEntries[start + j]];
      faceCenter += pA;
    }
    faceCenter /= D;
    faceCenters.data[iF] = faceCenter;
  }

  faceCenters.markHostBufferUpdated();
}

void SurfaceMesh::computeFaceAreas() {

  vertexPositions.ensureHostBufferPopulated();

  faceAreas.data.resize(nFaces());

  // Loop over faces to compute face-valued quantities
  for (size_t iF = 0; iF < nFaces(); iF++) {
    size_t start = faceIndsStart[iF];
    size_t D = faceIndsStart[iF + 1] - start;

    // Compute a face normal
    double fA;
    if (D == 3) {
      glm::vec3 pA = vertexPositions.data[faceIndsEntries[start + 0]];
      glm::vec3 pB = vertexPositions.data[faceIndsEntries[start + 1]];
      glm::vec3 pC = vertexPositions.data[faceIndsEntries[start + 2]];
      glm::vec3 fN = glm::cross(pB - pA, pC - pA);
      fA = 0.5 * glm::length(fN);
    } else {
      fA = 0;
      glm::vec3 pRoot = vertexPositions.data[faceIndsEntries[start]];
      for (size_t j = 1; j + 1 < D; j++) {
        glm::vec3 pA = vertexPositions.data[faceIndsEntries[start + j]];
        glm::vec3 pB = vertexPositions.data[faceIndsEntries[start + j + 1]];
        fA += 0.5 * glm::length(glm::cross(pA - pRoot, pB - pRoot));
      }
    }

    faceAreas.data[iF] = fA;
  }

  faceAreas.markHostBufferUpdated();
}


void SurfaceMesh::computeVertexNormals() {

  faceNormals.ensureHostBufferPopulated();
  faceAreas.ensureHostBufferPopulated();

  vertexNormals.data.resize(nVertices());

  const glm::vec3 zero{0., 0., 0.};

  std::fill(vertexNormals.data.begin(), vertexNormals.data.end(), zero);

  // Accumulate quantities from each face
  for (size_t iF = 0; iF < nFaces(); iF++) {
    size_t start = faceIndsStart[iF];
    size_t D = faceIndsStart[iF + 1] - start;
    for (size_t j = 0; j < D; j++) {
      size_t iV = faceIndsEntries[start + j];
      vertexNormals.data[iV] += faceNormals.data[iF] * static_cast<float>(faceAreas.data[iF]);
    }
  }

  // Normalize
  for (size_t iV = 0; iV < nVertices(); iV++) {
    vertexNormals.data[iV] = glm::normalize(vertexNormals.data[iV]);
  }

  vertexNormals.markHostBufferUpdated();
}

void SurfaceMesh::computeVertexAreas() {

  faceAreas.ensureHostBufferPopulated();

  vertexAreas.data.resize(nVertices());
  std::fill(vertexAreas.data.begin(), vertexAreas.data.end(), 0.);

  // Accumulate quantities from each face
  for (size_t iF = 0; iF < nFaces(); iF++) {
    size_t D = faceIndsStart[iF + 1] - faceIndsStart[iF];
    size_t start = faceIndsStart[iF];
    for (size_t j = 0; j < D; j++) {
      size_t iV = faceIndsEntries[start + j];
      vertexAreas.data[iV] += faceAreas.data[iF] / D;
    }
  }

  vertexAreas.markHostBufferUpdated();
}

void SurfaceMesh::computeDefaultFaceTangentBasisX() {

  // NOTE: this function is weirdly duplicated into an 'X' and 'Y' paradigm to fit the compute-function-per-buffer
  // paradigm

  vertexPositions.ensureHostBufferPopulated();
  faceNormals.ensureHostBufferPopulated();

  defaultFaceTangentBasisX.data.resize(nFaces());

  for (size_t iF = 0; iF < nFaces(); iF++) {
    size_t D = faceIndsStart[iF + 1] - faceIndsStart[iF];
    if (D != 3) exception("Default face tangent spaces only available for pure-triangular meshes");

    size_t start = faceIndsStart[iF];

    glm::vec3 pA = vertexPositions.data[faceIndsEntries[start + 0]];
    glm::vec3 pB = vertexPositions.data[faceIndsEntries[start + 1]];
    glm::vec3 N = faceNormals.data[iF];

    glm::vec3 basisX = pB - pA;
    basisX = glm::normalize(basisX - N * glm::dot(N, basisX));

    glm::vec3 basisY = glm::normalize(-glm::cross(basisX, N));

    defaultFaceTangentBasisX.data[iF] = basisX;
  }

  defaultFaceTangentBasisX.markHostBufferUpdated();
}

void SurfaceMesh::computeDefaultFaceTangentBasisY() {

  // NOTE: this function is weirdly duplicated into an 'X' and 'Y' paradigm to fit the compute-function-per-buffer
  // paradigm

  vertexPositions.ensureHostBufferPopulated();
  faceNormals.ensureHostBufferPopulated();

  defaultFaceTangentBasisY.data.resize(nFaces());

  for (size_t iF = 0; iF < nFaces(); iF++) {
    size_t D = faceIndsStart[iF + 1] - faceIndsStart[iF];
    if (D != 3) exception("Default face tangent spaces only available for pure-triangular meshes");

    size_t start = faceIndsStart[iF];

    glm::vec3 pA = vertexPositions.data[faceIndsEntries[start + 0]];
    glm::vec3 pB = vertexPositions.data[faceIndsEntries[start + 1]];
    glm::vec3 N = faceNormals.data[iF];

    glm::vec3 basisX = pB - pA;
    basisX = glm::normalize(basisX - N * glm::dot(N, basisX));

    glm::vec3 basisY = glm::normalize(-glm::cross(basisX, N));

    defaultFaceTangentBasisY.data[iF] = basisY;
  }

  defaultFaceTangentBasisY.markHostBufferUpdated();
}

// === Edge Lengths ===

// void SurfaceMesh::computeEdgeLengths() {
//
//   vertexPositions.ensureHostBufferPopulated();
//
//   edgeLengths.data.resize(nEdges());
//
//   // Compute edge lengths
//   for (size_t iF = 0; iF < nFaces(); iF++) {
//     size_t D = faceIndsStart[iF + 1] - faceIndsStart[iF];
//     size_t start = faceIndsStart[iF];
//     for (size_t j = 0; j < D; j++) {
//       size_t iA = faceIndsEntries[start + j];
//       size_t iB = faceIndsEntries[start + (j + 1) % D];
//       glm::vec3 pA = vertexPositions.data[iA];
//       glm::vec3 pB = vertexPositions.data[iB];
//       edgeLengths.data[edgeIndices[iF][j]] = glm::length(pA - pB);
//     }
//   }
//
//   edgeLengths.markHostBufferUpdated();
// }

void SurfaceMesh::checkTriangular() {
  if (nFacesTriangulation() != nFaces()) {
    exception("Cannot proceed, SurfaceMesh " + name + " is not a triangular mesh.");
  }
}

void SurfaceMesh::ensureHaveManifoldConnectivity() {
  if (!twinHalfedge.empty()) return; // already populated

  triangleVertexInds.ensureHostBufferPopulated();

  twinHalfedge.resize(nHalfedges());

  // Maps from edge (sorted) to all halfedges incident on that edge
  std::unordered_map<std::pair<size_t, size_t>, std::vector<size_t>,
                     polyscope::hash_combine::hash<std::pair<size_t, size_t>>>
      edgeInds;

  // Fill out faceForHalfedge and populate edge lookup map
  for (size_t iF = 0; iF < nFacesTriangulation(); iF++) {
    for (size_t j = 0; j < 3; j++) {
      size_t iV = triangleVertexInds.data[3 * iF + j];
      size_t iVNext = triangleVertexInds.data[3 * iF + ((j + 1) % 3)];
      size_t iHe = 3 * iF + j;


      std::pair<size_t, size_t> edgeKey(std::min(iV, iVNext), std::max(iV, iVNext));

      // Make sure the key is populated
      auto it = edgeInds.find(edgeKey);
      if (it == edgeInds.end()) {
        it = edgeInds.insert(it, {edgeKey, std::vector<size_t>()});
      }

      // Add this halfedge to the entry
      it->second.push_back(iHe);
    }
  }

  // Second walk through, setting twins
  for (size_t iF = 0; iF < nFacesTriangulation(); iF++) {
    for (size_t j = 0; j < 3; j++) {
      size_t iV = triangleVertexInds.data[3 * iF + j];
      size_t iVNext = triangleVertexInds.data[3 * iF + ((j + 1) % 3)];
      size_t iHe = 3 * iF + j;

      std::pair<size_t, size_t> edgeKey(std::min(iV, iVNext), std::max(iV, iVNext));
      std::vector<size_t>& edgeHalfedges = edgeInds.find(edgeKey)->second;

      // Pick the first halfedge we find which is not this one
      size_t myTwin = INVALID_IND;
      for (size_t t : edgeHalfedges) {
        if (t != iHe) {
          myTwin = t;
          break;
        }
      }

      twinHalfedge[iHe] = myTwin;
    }
  }
}

void SurfaceMesh::draw() {
  if (!isEnabled()) {
    return;
  }

  render::engine->setBackfaceCull(backFacePolicy.get() == BackFacePolicy::Cull);

  // If no quantity is drawing the surface, we should draw it
  if (dominantQuantity == nullptr) {

    if (program == nullptr) {
      prepare();

      // do this now to reduce lag when picking later, etc
      // FIXME
      // preparePick();
    }

    // Set uniforms
    setStructureUniforms(*program);
    setSurfaceMeshUniforms(*program);
    program->setUniform("u_baseColor", getSurfaceColor());
    render::engine->setMaterialUniforms(*program, getMaterial());

    program->draw();
  }

  // Draw the quantities
  for (auto& x : quantities) {
    x.second->draw();
  }

  render::engine->setBackfaceCull(); // return to default setting

  for (auto& x : floatingQuantities) {
    x.second->draw();
  }
}

void SurfaceMesh::drawDelayed() {
  if (!isEnabled()) {
    return;
  }

  render::engine->setBackfaceCull(backFacePolicy.get() == BackFacePolicy::Cull);

  for (auto& x : quantities) {
    x.second->drawDelayed();
  }

  render::engine->setBackfaceCull(); // return to default setting

  for (auto& x : floatingQuantities) {
    x.second->drawDelayed();
  }
}

void SurfaceMesh::drawPick() {
  if (!isEnabled()) {
    return;
  }

  if (pickProgram == nullptr) {
    preparePick();
  }

  render::engine->setBackfaceCull(backFacePolicy.get() == BackFacePolicy::Cull);

  // Set uniforms
  setStructureUniforms(*pickProgram);

  pickProgram->draw();

  render::engine->setBackfaceCull(); // return to default setting
}

void SurfaceMesh::prepare() {
  // clang-format off
  program = render::engine->requestShader( "MESH", 
      render::engine->addMaterialRules(getMaterial(),
        addSurfaceMeshRules({"SHADE_BASECOLOR"})
      )
  );
  // clang-format on

  // Populate draw buffers
  setMeshGeometryAttributes(*program);
  render::engine->setMaterial(*program, getMaterial());
}

void SurfaceMesh::preparePick() {


  bool simplePick = !(edgesHaveBeenUsed || halfedgesHaveBeenUsed || cornersHaveBeenUsed);

  if (simplePick) {
    pickProgram =
        render::engine->requestShader("MESH", addSurfaceMeshRules({"MESH_PROPAGATE_PICK_SIMPLE"}, true, false),
                                      render::ShaderReplacementDefaults::Pick);
  } else {
    pickProgram = render::engine->requestShader("MESH", addSurfaceMeshRules({"MESH_PROPAGATE_PICK"}, true, false),
                                                render::ShaderReplacementDefaults::Pick);
  }

  // Populate draw buffers
  setMeshGeometryAttributes(*pickProgram);
  setMeshPickAttributes(*pickProgram);
}

void SurfaceMesh::setMeshGeometryAttributes(render::ShaderProgram& p) {
  if (p.hasAttribute("a_vertexPositions")) {
    p.setAttribute("a_vertexPositions", vertexPositions.getIndexedRenderAttributeBuffer(triangleVertexInds));
  }
  if (p.hasAttribute("a_vertexNormals")) {

    if (getShadeStyle() == MeshShadeStyle::Smooth) {
      p.setAttribute("a_vertexNormals", vertexNormals.getIndexedRenderAttributeBuffer(triangleVertexInds));
    } else {
      // these aren't actually used in in the automatically-generated case, but the shader is set up in a lazy way so
      // it is still needed
      p.setAttribute("a_vertexNormals", faceNormals.getIndexedRenderAttributeBuffer(triangleFaceInds));
    }
  }
  if (p.hasAttribute("a_normal")) {
    p.setAttribute("a_normal", faceNormals.getIndexedRenderAttributeBuffer(triangleFaceInds));
  }
  if (p.hasAttribute("a_barycoord")) {
    p.setAttribute("a_barycoord", baryCoord.getRenderAttributeBuffer());
  }
  if (p.hasAttribute("a_edgeIsReal")) {
    p.setAttribute("a_edgeIsReal", edgeIsReal.getRenderAttributeBuffer());
  }
  if (wantsCullPosition()) {
    p.setAttribute("a_cullPos", faceCenters.getIndexedRenderAttributeBuffer(triangleFaceInds));
  }

  if (transparencyQuantityName != "") {
    SurfaceScalarQuantity& transparencyQ = resolveTransparencyQuantity();
    p.setAttribute("a_valueAlpha", transparencyQ.getAttributeBuffer());
  }
}

void SurfaceMesh::setMeshPickAttributes(render::ShaderProgram& p) {

  // TODO in principle all of the data this shader needs is already available on the GPU via the [...]Inds attribute
  // buffers. We could move the encoding / offsetting logic to happen in a shader with some uniforms, and avoid any
  // CPU-side processing. Maybe the solution is to directly render ints?

  // make sure we have the relevant indexing data
  bool simplePick = !(edgesHaveBeenUsed || halfedgesHaveBeenUsed || cornersHaveBeenUsed);
  triangleVertexInds.ensureHostBufferPopulated();
  triangleFaceInds.ensureHostBufferPopulated();
  if (edgesHaveBeenUsed) triangleAllEdgeInds.ensureHostBufferPopulated();
  if (halfedgesHaveBeenUsed) triangleAllHalfedgeInds.ensureHostBufferPopulated();
  if (cornersHaveBeenUsed) triangleCornerInds.ensureHostBufferPopulated();

  // nEdges() requires computing number of edges, which is expensive and might not even be implemented for polygonal
  // meshes. This way we only call it if actually needed, and use 0 otherwise.
  size_t nEdgesSafe = edgesHaveBeenUsed ? nEdges() : 0;

  // Get element indices
  size_t totalPickElements = nVertices() + nFaces() + nEdgesSafe + nHalfedges() + nCorners();

  // In "local" indices, indexing elements only within this mesh, used for reading later
  facePickIndStart = nVertices();
  edgePickIndStart = facePickIndStart + nFaces();
  halfedgePickIndStart = edgePickIndStart + nEdgesSafe;
  cornerPickIndStart = halfedgePickIndStart + nHalfedges();

  // In "global" indices, indexing all elements in the scene, used to fill buffers for drawing here
  size_t pickStart = pick::requestPickBufferRange(this, totalPickElements);
  size_t vertexGlobalPickIndStart = pickStart;
  size_t faceGlobalPickIndStart = pickStart + facePickIndStart;
  size_t edgeGlobalPickIndStart = pickStart + edgePickIndStart;
  size_t halfedgeGlobalPickIndStart = pickStart + halfedgePickIndStart;
  size_t cornerGlobalPickIndStart = pickStart + cornerPickIndStart;

  // == Fill buffers
  std::vector<std::array<glm::vec3, 3>> vertexColors, halfedgeColors, cornerColors;
  std::vector<glm::vec3> faceColor;

  // Reserve space
  vertexColors.reserve(3 * nFacesTriangulation());
  faceColor.reserve(3 * nFacesTriangulation());
  if (!simplePick) {
    halfedgeColors.reserve(3 * nFacesTriangulation());
    cornerColors.reserve(3 * nFacesTriangulation());
  }


  // Build all quantities in each face
  size_t iFTri = 0;
  for (size_t iF = 0; iF < nFaces(); iF++) {
    size_t D = faceIndsStart[iF + 1] - faceIndsStart[iF];

    glm::vec3 fColor = pick::indToVec(iF + faceGlobalPickIndStart);

    for (size_t j = 1; (j + 1) < D; j++) {

      // == Build face & vertex index data

      // clang-format off
      std::array<glm::vec3, 3> vColor = {
        pick::indToVec(triangleVertexInds.data[3*iFTri + 0] + vertexGlobalPickIndStart),
        pick::indToVec(triangleVertexInds.data[3*iFTri + 1] + vertexGlobalPickIndStart),
        pick::indToVec(triangleVertexInds.data[3*iFTri + 2] + vertexGlobalPickIndStart),
      };
      // clang-format on

      for (int j = 0; j < 3; j++) {
        faceColor.push_back(fColor);
        vertexColors.push_back(vColor);
      }

      // Second half does halfedges/edges/corners, not used for simple mode
      if (simplePick) {
        iFTri++;
        continue;
      }


      // Fill the halfedge buffer with edge or halfedge data, depending on which are in use
      // In the pick function we will use the halfedge to look up the edge if needed
      // (this is an optimization to use one less array of values, because we hit implementation limits in the shader)

      // == Build edge index data, if needed


      if (edgesHaveBeenUsed || halfedgesHaveBeenUsed) {

        const std::vector<uint32_t>& eDataVec =
            (edgesHaveBeenUsed && !halfedgesHaveBeenUsed) ? triangleAllEdgeInds.data : triangleAllHalfedgeInds.data;
        size_t offset =
            (edgesHaveBeenUsed && !halfedgesHaveBeenUsed) ? edgeGlobalPickIndStart : halfedgeGlobalPickIndStart;

        // clang-format off
        std::array<glm::vec3, 3> eColor = { 
          fColor, 
          pick::indToVec(eDataVec[9*iFTri + 1] + offset), 
          fColor
        };
        // clang-format on
        if (j == 1) eColor[0] = pick::indToVec(eDataVec[9 * iFTri + 0] + offset);
        if (j + 2 == D) eColor[2] = pick::indToVec(eDataVec[9 * iFTri + 2] + offset);

        for (int j = 0; j < 3; j++) halfedgeColors.push_back(eColor);
      } else {
        for (int j = 0; j < 3; j++) halfedgeColors.push_back({fColor, fColor, fColor});
      }

      // == Build corner index data, if needed

      if (cornersHaveBeenUsed) {
        // clang-format off
        std::array<glm::vec3, 3> cColor = { 
          pick::indToVec(triangleCornerInds.data[3*iFTri + 0] + cornerGlobalPickIndStart), 
          pick::indToVec(triangleCornerInds.data[3*iFTri + 1] + cornerGlobalPickIndStart), 
          pick::indToVec(triangleCornerInds.data[3*iFTri + 2] + cornerGlobalPickIndStart), 
        };
        // clang-format on
        for (int j = 0; j < 3; j++) cornerColors.push_back(cColor);
      } else {
        for (int j = 0; j < 3; j++) cornerColors.push_back({vColor[0], vColor[1], vColor[2]});
      }

      iFTri++;
    }
  }

  // == Store data in buffers

  std::shared_ptr<render::AttributeBuffer> vertexColorsBuff =
      render::engine->generateAttributeBuffer(RenderDataType::Vector3Float, 3);
  vertexColorsBuff->setData(vertexColors);
  pickProgram->setAttribute("a_vertexColors", vertexColorsBuff);

  std::shared_ptr<render::AttributeBuffer> faceColorsBuff =