Volumetric grids and scalar isosurfaces

include/polyscope/marchingcubes/Vectors.h

@@ -0,0 +1,38 @@
+#pragma once
+
+#ifndef VECTORS_H
+#define VECTORS_H
+
+namespace polyscope {
+namespace marchingcubes {
+// File Name: Vectors.h
+// Last Modified: 7/8/2000
+// Author: Raghavendra Chandrashekara
+// Email: rc99@doc.ic.ac.uk, rchandrashekara@hotmail.com
+//
+// Description: This file contains some useful structures.
+
+// Edited and ported to polyscope by Christopher Yu.
+
+typedef float POINT3D[3];
+typedef float VECTOR3D[3];
+
+struct POINT3DXYZ {
+  float x, y, z;
+  friend POINT3DXYZ operator+(const POINT3DXYZ& pt3dPoint1, const POINT3DXYZ& pt3dPoint2);
+  friend POINT3DXYZ operator-(const POINT3DXYZ& pt3dPoint1, const POINT3DXYZ& pt3dPoint2);
+  friend POINT3DXYZ operator*(const POINT3DXYZ& pt3dPoint, float fScale);
+  friend POINT3DXYZ operator*(float fScale, const POINT3DXYZ& pt3dPoint);
+  friend POINT3DXYZ operator/(const POINT3DXYZ& pt3dPoint, float fScale);
+  friend POINT3DXYZ& operator*=(POINT3DXYZ& pt3dPoint, float fScale);
+  friend POINT3DXYZ& operator/=(POINT3DXYZ& pt3dPoint, float fScale);
+  friend POINT3DXYZ& operator+=(POINT3DXYZ& pt3dPoint1, const POINT3DXYZ& pt3dPoint2);
+  friend POINT3DXYZ& operator-=(POINT3DXYZ& pt3dPoint1, const POINT3DXYZ& pt3dPoint2);
+};
+
+typedef POINT3DXYZ VECTOR3DXYZ;
+
+} // namespace marchingcubes
+} // namespace polyscope
+
+#endif // VECTORS_H

include/polyscope/marchingcubes/mesh_implicit_surface.h

@@ -0,0 +1,69 @@
+#pragma once
+
+#include <array>
+#include <glm/glm.hpp>
+#include <vector>
+
+#include "CIsoSurface.h"
+
+namespace polyscope {
+namespace marchingcubes {
+
+template <typename Implicit, typename Data>
+void SampleFunctionToGrid(const Implicit &funct, size_t numCornersPerSide, glm::vec3 center, double sideLength, std::vector<Data> &field) {
+  double diameter = sideLength;
+  double cellSize = diameter / (numCornersPerSide - 1);
+  double radius = diameter / 2;
+
+  glm::vec3 lowerCorner = center - glm::vec3{radius, radius, radius};
+
+  int nSlice = numCornersPerSide * numCornersPerSide;
+  int nRow = numCornersPerSide;
+
+  for (size_t x = 0; x < numCornersPerSide; x++) {
+    for (size_t y = 0; y < numCornersPerSide; y++) {
+      for (size_t z = 0; z < numCornersPerSide; z++) {
+        glm::vec3 samplePt = lowerCorner + glm::vec3{(double)x, (double)y, (double)z} * (float)cellSize;
+        Data value = funct.ValueAt(samplePt);
+        field[nSlice * z + nRow * y + x] = value;
+      }
+    }
+  }
+}
+
+inline void MeshImplicitGrid(std::vector<double> &field, double isoLevel, size_t numCornersPerSide, glm::vec3 center, double sideLength,
+                             std::vector<glm::vec3> &nodes, std::vector<std::array<size_t, 3>> &triangles) {
+  CIsoSurface<double>* iso = new CIsoSurface<double>();
+  int numCells = numCornersPerSide - 1;
+  double diameter = sideLength;
+  double cellSize = diameter / numCells;
+  double radius = diameter / 2;
+  glm::vec3 lowerCorner = center - glm::vec3{radius, radius, radius};
+
+  iso->GenerateSurface(&field[0], isoLevel, numCells, numCells, numCells, cellSize, cellSize, cellSize);
+
+  int nVerts = iso->m_nVertices;
+
+  for (int i = 0; i < nVerts; i++) {
+    double x = iso->m_ppt3dVertices[i][0];
+    double y = iso->m_ppt3dVertices[i][1];
+    double z = iso->m_ppt3dVertices[i][2];
+
+    glm::vec3 p = lowerCorner + glm::vec3{x, y, z};
+    nodes.push_back(p);
+  }
+
+  int nTris = iso->m_nTriangles;
+
+  for (int i = 0; i < nTris; i++) {
+    int i1 = iso->m_piTriangleIndices[3 * i];
+    int i2 = iso->m_piTriangleIndices[3 * i + 1];
+    int i3 = iso->m_piTriangleIndices[3 * i + 2];
+
+    triangles.push_back({(size_t)i1, (size_t)i2, (size_t)i3});
+  }
+  delete iso;
+}
+
+} // namespace marchingcubes
+} // namespace polyscope

include/polyscope/volumetric_grid.h

@@ -0,0 +1,141 @@
+// Copyright 2017-2019, Nicholas Sharp and the Polyscope contributors. http://polyscope.run.
+#pragma once
+
+#include "polyscope/affine_remapper.h"
+#include "polyscope/color_management.h"
+#include "polyscope/gl/gl_utils.h"
+#include "polyscope/polyscope.h"
+#include "polyscope/standardize_data_array.h"
+#include "polyscope/structure.h"
+
+#include "polyscope/volumetric_grid_quantity.h"
+
+#include <vector>
+
+#include "marchingcubes/mesh_implicit_surface.h"
+
+namespace polyscope {
+
+class VolumetricGrid;
+class VolumetricGridScalarIsosurface;
+class VolumetricGridScalarQuantity;
+class VolumetricGridVectorQuantity;
+
+template <> // Specialize the quantity type
+struct QuantityTypeHelper<VolumetricGrid> {
+  typedef VolumetricGridQuantity type;
+};
+
+
+class VolumetricGrid : public QuantityStructure<VolumetricGrid> {
+public:
+  // Construct a new curve network structure
+  VolumetricGrid(std::string name, size_t nValuesPerSide, glm::vec3 center, double sideLen);
+
+  // === Overloads
+
+  // Build the imgui display
+  virtual void buildCustomUI() override;
+  virtual void buildPickUI(size_t localPickID) override;
+  // Render the the structure on screen
+  virtual void draw() override;
+  // Render for picking
+  virtual void drawPick() override;
+  // A characteristic length for the structure
+  virtual double lengthScale() override;
+  // Axis-aligned bounding box for the structure
+  virtual std::tuple<glm::vec3, glm::vec3> boundingBox() override;
+  virtual std::string typeName() override;
+
+  // Field data
+  size_t nCornersPerSide;
+  glm::vec3 gridCenter;
+  double sideLength;
+
+  inline size_t nValues() const { return nCornersPerSide * nCornersPerSide * nCornersPerSide; }
+
+  inline glm::vec3 positionOfIndex(size_t i) {
+    size_t nPerSlice = nCornersPerSide * nCornersPerSide;
+    size_t z = i / nPerSlice;
+    size_t i_in_slice = i % nPerSlice;
+    size_t y = i_in_slice / nCornersPerSide;
+    size_t x = i_in_slice % nCornersPerSide;
+
+    double cellSize = sideLength / (nCornersPerSide - 1);
+    double radius = sideLength / 2;
+    glm::vec3 lowerCorner = gridCenter - glm::vec3{radius, radius, radius};
+    return lowerCorner + glm::vec3{x * cellSize, y * cellSize, z * cellSize};
+  }
+
+  // Misc data
+  static const std::string structureTypeName;
+
+  template <class T>
+  VolumetricGridScalarIsosurface* addGridIsosurfaceQuantity(std::string name, double isoLevel, const T& values) {
+    validateSize(values, nValues(), "grid isosurface quantity " + name);
+    return addIsosurfaceQuantityImpl(name, isoLevel, standardizeArray<double, T>(values));
+  }
+
+  template <class Funct>
+  VolumetricGridScalarIsosurface* addGridIsosurfaceQuantityFromFunction(std::string name, double isoLevel, const Funct& funct) {
+    size_t totalValues = nCornersPerSide * nCornersPerSide * nCornersPerSide;
+    std::vector<double> field(totalValues);
+    marchingcubes::SampleFunctionToGrid(funct, nCornersPerSide, gridCenter, sideLength, field);
+    return addGridIsosurfaceQuantity(name, isoLevel, field);
+  }
+
+  template <class T>
+  VolumetricGridScalarQuantity* addGridScalarQuantity(std::string name, const T& values, DataType dataType_) {
+    validateSize(values, nValues(), "grid scalar quantity " + name);
+    return addScalarQuantityImpl(name, standardizeArray<double, T>(values), dataType_);
+  }
+
+  template <class Funct>
+  VolumetricGridScalarQuantity* addGridScalarQuantityFromFunction(std::string name, const Funct& funct, DataType dataType_) {
+    size_t totalValues = nCornersPerSide * nCornersPerSide * nCornersPerSide;
+    std::vector<double> field(totalValues);
+    marchingcubes::SampleFunctionToGrid(funct, nCornersPerSide, gridCenter, sideLength, field);
+    return addGridScalarQuantity(name, field, dataType_);
+  }
+
+  template <class T>
+  VolumetricGridVectorQuantity* addGridVectorQuantity(std::string name, const T& vecValues, VectorType dataType_) {
+    validateSize(vecValues, nValues(), "grid vector quantity " + name);
+    return addVectorQuantityImpl(name, standardizeArray<glm::vec3, T>(vecValues), dataType_);
+  }
+
+  template <class Funct>
+  VolumetricGridVectorQuantity* addGridVectorQuantityFromFunction(std::string name, const Funct& funct, VectorType dataType_) {
+    size_t totalValues = nCornersPerSide * nCornersPerSide * nCornersPerSide;
+    std::vector<glm::vec3> field(totalValues);
+    marchingcubes::SampleFunctionToGrid(funct, nCornersPerSide, gridCenter, sideLength, field);
+    return addGridVectorQuantity(name, field, dataType_);
+  }
+
+
+  glm::vec3 getColor();
+
+private:
+  PersistentValue<glm::vec3> color;
+  VolumetricGrid* setColor(glm::vec3 newVal);
+  VolumetricGridScalarIsosurface* addIsosurfaceQuantityImpl(std::string name, double isoLevel,
+                                                            const std::vector<double>& data);
+  VolumetricGridScalarQuantity* addScalarQuantityImpl(std::string name, const std::vector<double>& data,
+                                                      DataType dataType_);
+  VolumetricGridVectorQuantity* addVectorQuantityImpl(std::string name, const std::vector<glm::vec3>& data,
+                                                      VectorType dataType_);
+};
+
+
+VolumetricGrid* registerVolumetricGrid(std::string name, size_t nValuesPerSide, glm::vec3 center, double sideLen);
+
+template <typename Implicit>
+VolumetricGrid* registerIsosurfaceFromFunction(std::string name, const Implicit& funct, size_t nValuesPerSide,
+                                               glm::vec3 center, double sideLen, bool meshImmediately = true) {
+
+  VolumetricGrid* outputSurface = registerVolumetricGrid(name, nValuesPerSide, center, sideLen);
+  outputSurface->addGridIsosurfaceQuantityFromFunction("isosurface", 0, funct);
+  return outputSurface;
+}
+
+} // namespace polyscope

include/polyscope/volumetric_grid_quantity.h

@@ -0,0 +1,27 @@
+// Copyright 2017-2019, Nicholas Sharp and the Polyscope contributors. http://polyscope.run.
+#pragma once
+
+#include "polyscope/quantity.h"
+#include "polyscope/structure.h"
+
+namespace polyscope {
+
+// Forward declare
+class VolumetricGrid;
+
+// Extend Quantity<CurveNetwork> 
+class VolumetricGridQuantity : public Quantity<VolumetricGrid> {
+public:
+  VolumetricGridQuantity(std::string name, VolumetricGrid& parentStructure, bool dominates = false);
+  virtual ~VolumetricGridQuantity() {};
+
+  // Build GUI info an element
+  virtual void buildNodeInfoGUI(size_t vInd);
+  virtual void buildEdgeInfoGUI(size_t fInd);
+
+  // Invalidate geometric data
+  virtual void geometryChanged() = 0;
+};
+
+
+} // namespace polyscope

include/polyscope/volumetric_grid_scalar_isosurface.h

@@ -0,0 +1,37 @@
+#pragma once
+
+#include "polyscope/affine_remapper.h"
+#include "polyscope/gl/color_maps.h"
+#include "polyscope/histogram.h"
+#include "polyscope/volumetric_grid.h"
+
+namespace polyscope {
+
+class VolumetricGridScalarIsosurface : public VolumetricGridQuantity {
+public:
+  VolumetricGridScalarIsosurface(std::string name, VolumetricGrid& grid_, const std::vector<double> &values_, double levelSet_);
+
+  virtual void draw() override;
+  virtual void buildCustomUI() override;
+  virtual std::string niceName() override;
+  virtual void geometryChanged() override;
+
+  void recomputeNormals();
+  void meshIsosurface();
+
+protected:
+  std::vector<double> values;
+  double levelSet;
+  double newLevelSet;
+  double increment;
+  void prepare();
+  void prepareTriangleIndices();
+  void setProgramUniforms(gl::GLProgram& program);
+  void fillGeometryBuffersSmooth(gl::GLProgram& p);
+  std::unique_ptr<gl::GLProgram> meshProgram;
+  std::vector<glm::vec3> nodes;
+  std::vector<glm::vec3> normals;
+  std::vector<std::array<size_t, 3>> triangles;
+};
+
+} // namespace polyscope

include/polyscope/volumetric_grid_scalar_quantity.h

@@ -0,0 +1,36 @@
+#pragma once
+
+#include "polyscope/affine_remapper.h"
+#include "polyscope/gl/color_maps.h"
+#include "polyscope/histogram.h"
+#include "polyscope/volumetric_grid.h"
+
+namespace polyscope {
+
+class VolumetricGridScalarQuantity : public VolumetricGridQuantity {
+public:
+  VolumetricGridScalarQuantity(std::string name, VolumetricGrid& grid_, const std::vector<double>& values_, DataType dataType_);
+
+  virtual void draw() override;
+  virtual void buildCustomUI() override;
+  virtual std::string niceName() override;
+  virtual void geometryChanged() override;
+
+protected:
+  void createProgram();
+  void fillGeometryBuffersSmooth(gl::GLProgram& p);
+  void fillPositions();
+  void setPointCloudUniforms();
+  void resetMapRange();
+
+  const DataType dataType;
+  std::vector<double> values;
+  std::vector<glm::vec3> positions;
+  PersistentValue<gl::ColorMapID> cMap;
+  std::unique_ptr<gl::GLProgram> pointsProgram;
+  float pointRadius;
+  std::pair<float, float> vizRange;
+  std::pair<double, double> dataRange;
+};
+
+} // namespace polyscope

include/polyscope/volumetric_grid_vector_quantity.h

@@ -0,0 +1,31 @@
+#pragma once
+
+#include "polyscope/affine_remapper.h"
+#include "polyscope/gl/color_maps.h"
+#include "polyscope/histogram.h"
+#include "polyscope/volumetric_grid.h"
+
+namespace polyscope {
+class VolumetricGridVectorQuantity : public VolumetricGridQuantity {
+public:
+  VolumetricGridVectorQuantity(std::string name, VolumetricGrid& grid_, const std::vector<glm::vec3>& values_,
+                               VectorType vectorType_);
+
+  virtual void draw() override;
+  virtual void buildCustomUI() override;
+  virtual std::string niceName() override;
+  virtual void geometryChanged() override;
+
+protected:
+  VectorType vectorType;
+  std::vector<glm::vec3> vectorValues;
+  std::vector<glm::vec3> positions;
+  std::unique_ptr<gl::GLProgram> vectorsProgram;
+  AffineRemapper<glm::vec3> mapper;
+  float vectorRadius;
+  float vectorLengthMult;
+
+  void fillPositions();
+  void createProgram();
+};
+} // namespace polyscope