diff --git a/include/polyscope/marchingcubes/CIsoSurface.h b/include/polyscope/marchingcubes/CIsoSurface.h new file mode 100644 index 00000000..ac3a899b --- /dev/null +++ b/include/polyscope/marchingcubes/CIsoSurface.h @@ -0,0 +1,844 @@ +#pragma once + +#ifndef CISOSURFACE_H +#define CISOSURFACE_H + +// File Name: CIsoSurface.h +// Last Modified: 5/8/2000 +// Author: Raghavendra Chandrashekara (basesd on source code +// provided by Paul Bourke and Cory Gene Bloyd) +// Email: rc99@doc.ic.ac.uk, rchandrashekara@hotmail.com +// +// Description: This is the interface file for the CIsoSurface class. +// CIsoSurface can be used to construct an isosurface from a scalar +// field. + +// Edited and ported to polyscope by Christopher Yu. + +#include "Vectors.h" +#include +#include +#include + +namespace polyscope { +namespace marchingcubes { + +struct POINT3DID { + unsigned int newID; + float x, y, z; +}; + +typedef std::map ID2POINT3DID; + +struct TRIANGLE { + unsigned int pointID[3]; +}; + +typedef std::vector TRIANGLEVECTOR; + +template +class CIsoSurface { +public: + // Constructor and destructor. + CIsoSurface(); + ~CIsoSurface(); + + // Generates the isosurface from the scalar field contained in the + // buffer ptScalarField[]. + void GenerateSurface(const T* ptScalarField, T tIsoLevel, unsigned int nCellsX, unsigned int nCellsY, + unsigned int nCellsZ, float fCellLengthX, float fCellLengthY, float fCellLengthZ); + + // Returns true if a valid surface has been generated. + bool IsSurfaceValid(); + + // Deletes the isosurface. + void DeleteSurface(); + + // Returns the length, width, and height of the volume in which the + // isosurface in enclosed in. Returns -1 if the surface is not + // valid. + int GetVolumeLengths(float& fVolLengthX, float& fVolLengthY, float& fVolLengthZ); + + // The number of vertices which make up the isosurface. + unsigned int m_nVertices; + + // The vertices which make up the isosurface. + POINT3D* m_ppt3dVertices; + + // The number of triangles which make up the isosurface. + unsigned int m_nTriangles; + + // The indices of the vertices which make up the triangles. + unsigned int* m_piTriangleIndices; + + // The number of normals. + unsigned int m_nNormals; + + // The normals. + VECTOR3D* m_pvec3dNormals; + + // List of POINT3Ds which form the isosurface. + ID2POINT3DID m_i2pt3idVertices; + + // List of TRIANGLES which form the triangulation of the isosurface. + TRIANGLEVECTOR m_trivecTriangles; + + // Returns the edge ID. + unsigned int GetEdgeID(unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo); + + // Returns the vertex ID. + unsigned int GetVertexID(unsigned int nX, unsigned int nY, unsigned int nZ); + + // Calculates the intersection point of the isosurface with an + // edge. + POINT3DID CalculateIntersection(unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo); + + // Interpolates between two grid points to produce the point at which + // the isosurface intersects an edge. + POINT3DID Interpolate(float fX1, float fY1, float fZ1, float fX2, float fY2, float fZ2, T tVal1, T tVal2); + + // Renames vertices and triangles so that they can be accessed more + // efficiently. + void RenameVerticesAndTriangles(); + + // Calculates the normals. + void CalculateNormals(); + + // No. of cells in x, y, and z directions. + unsigned int m_nCellsX, m_nCellsY, m_nCellsZ; + + // Cell length in x, y, and z directions. + float m_fCellLengthX, m_fCellLengthY, m_fCellLengthZ; + + // The buffer holding the scalar field. + const T* m_ptScalarField; + + // The isosurface value. + T m_tIsoLevel; + + // Indicates whether a valid surface is present. + bool m_bValidSurface; + + // Lookup tables used in the construction of the isosurface. + static const unsigned int m_edgeTable[256]; + static const int m_triTable[256][16]; +}; + +template +const unsigned int CIsoSurface::m_edgeTable[256] = { + 0x0, 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c, 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00, + 0x190, 0x99, 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c, 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90, + 0x230, 0x339, 0x33, 0x13a, 0x636, 0x73f, 0x435, 0x53c, 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30, + 0x3a0, 0x2a9, 0x1a3, 0xaa, 0x7a6, 0x6af, 0x5a5, 0x4ac, 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0, + 0x460, 0x569, 0x663, 0x76a, 0x66, 0x16f, 0x265, 0x36c, 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60, + 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff, 0x3f5, 0x2fc, 0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0, + 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55, 0x15c, 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950, + 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc, 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0, + 0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc, 0xcc, 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0, + 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c, 0x15c, 0x55, 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650, + 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc, 0x2fc, 0x3f5, 0xff, 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0, + 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c, 0x36c, 0x265, 0x16f, 0x66, 0x76a, 0x663, 0x569, 0x460, + 0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac, 0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa, 0x1a3, 0x2a9, 0x3a0, + 0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c, 0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33, 0x339, 0x230, + 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c, 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99, 0x190, + 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c, 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0}; + +template +const int CIsoSurface::m_triTable[256][16] = {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1}, + {3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1}, + {3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1}, + {3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1}, + {9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1}, + {1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1}, + {9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1}, + {2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1}, + {8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1}, + {9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1}, + {4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1}, + {3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1}, + {1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1}, + {4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1}, + {4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1}, + {9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1}, + {1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1}, + {5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1}, + {2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1}, + {9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1}, + {0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1}, + {2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1}, + {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1}, + {4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1}, + {5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1}, + {5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1}, + {9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1}, + {0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1}, + {1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1}, + {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1}, + {8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1}, + {2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1}, + {7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1}, + {9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1}, + {2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1}, + {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1}, + {9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1}, + {5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1}, + {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1}, + {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1}, + {1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1}, + {9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1}, + {5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1}, + {2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1}, + {0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1}, + {5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1}, + {6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1}, + {0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1}, + {3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1}, + {6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1}, + {5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1}, + {1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1}, + {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1}, + {6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1}, + {1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1}, + {8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1}, + {7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1}, + {3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1}, + {5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1}, + {0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1}, + {9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1}, + {8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1}, + {5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1}, + {0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1}, + {6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1}, + {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1}, + {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1}, + {8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1}, + {1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1}, + {3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1}, + {0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1}, + {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1}, + {0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1}, + {3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1}, + {6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1}, + {9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1}, + {8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1}, + {3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1}, + {6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1}, + {0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1}, + {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1}, + {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1}, + {1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1}, + {2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1}, + {7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1}, + {7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1}, + {2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1}, + {1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1}, + {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1}, + {8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1}, + {0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1}, + {7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1}, + {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1}, + {2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1}, + {6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1}, + {7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1}, + {2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1}, + {1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1}, + {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1}, + {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1}, + {0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1}, + {7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1}, + {6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1}, + {8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1}, + {9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1}, + {6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1}, + {1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1}, + {4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1}, + {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1}, + {8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1}, + {0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1}, + {1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1}, + {8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1}, + {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1}, + {4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1}, + {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1}, + {5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1}, + {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1}, + {9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1}, + {6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1}, + {7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1}, + {3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1}, + {7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1}, + {9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1}, + {3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1}, + {6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1}, + {9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1}, + {1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1}, + {4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1}, + {7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1}, + {6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1}, + {3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1}, + {0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1}, + {6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1}, + {1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1}, + {0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1}, + {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1}, + {6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1}, + {5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1}, + {9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1}, + {1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1}, + {1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1}, + {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1}, + {0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1}, + {5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1}, + {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1}, + {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1}, + {0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1}, + {9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1}, + {7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1}, + {2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1}, + {8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1}, + {9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1}, + {9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1}, + {1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1}, + {9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1}, + {9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1}, + {5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1}, + {0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1}, + {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1}, + {2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1}, + {0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1}, + {0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1}, + {9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1}, + {5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1}, + {3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1}, + {5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1}, + {8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1}, + {0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1}, + {9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1}, + {0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1}, + {1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1}, + {3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1}, + {4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1}, + {9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1}, + {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1}, + {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1}, + {2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1}, + {9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1}, + {3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1}, + {1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1}, + {4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1}, + {4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1}, + {0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1}, + {3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1}, + {3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1}, + {0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1}, + {9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1}, + {1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, + {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}}; + +template +CIsoSurface::CIsoSurface() { + m_fCellLengthX = 0; + m_fCellLengthY = 0; + m_fCellLengthZ = 0; + m_nCellsX = 0; + m_nCellsY = 0; + m_nCellsZ = 0; + m_nTriangles = 0; + m_nNormals = 0; + m_nVertices = 0; + m_ppt3dVertices = NULL; + m_piTriangleIndices = NULL; + m_pvec3dNormals = NULL; + m_ptScalarField = NULL; + m_tIsoLevel = 0; + m_bValidSurface = false; +} + +template +CIsoSurface::~CIsoSurface() { + DeleteSurface(); +} + +template +void CIsoSurface::GenerateSurface(const T* ptScalarField, T tIsoLevel, unsigned int nCellsX, unsigned int nCellsY, + unsigned int nCellsZ, float fCellLengthX, float fCellLengthY, float fCellLengthZ) { + if (m_bValidSurface) DeleteSurface(); + + m_tIsoLevel = tIsoLevel; + m_nCellsX = nCellsX; + m_nCellsY = nCellsY; + m_nCellsZ = nCellsZ; + m_fCellLengthX = fCellLengthX; + m_fCellLengthY = fCellLengthY; + m_fCellLengthZ = fCellLengthZ; + m_ptScalarField = ptScalarField; + + unsigned int nPointsInXDirection = (m_nCellsX + 1); + unsigned int nPointsInSlice = nPointsInXDirection * (m_nCellsY + 1); + + // Generate isosurface. + for (unsigned int z = 0; z < m_nCellsZ; z++) + for (unsigned int y = 0; y < m_nCellsY; y++) + for (unsigned int x = 0; x < m_nCellsX; x++) { + // Calculate table lookup index from those + // vertices which are below the isolevel. + unsigned int tableIndex = 0; + if (m_ptScalarField[z * nPointsInSlice + y * nPointsInXDirection + x] < m_tIsoLevel) tableIndex |= 1; + if (m_ptScalarField[z * nPointsInSlice + (y + 1) * nPointsInXDirection + x] < m_tIsoLevel) tableIndex |= 2; + if (m_ptScalarField[z * nPointsInSlice + (y + 1) * nPointsInXDirection + (x + 1)] < m_tIsoLevel) + tableIndex |= 4; + if (m_ptScalarField[z * nPointsInSlice + y * nPointsInXDirection + (x + 1)] < m_tIsoLevel) tableIndex |= 8; + if (m_ptScalarField[(z + 1) * nPointsInSlice + y * nPointsInXDirection + x] < m_tIsoLevel) tableIndex |= 16; + if (m_ptScalarField[(z + 1) * nPointsInSlice + (y + 1) * nPointsInXDirection + x] < m_tIsoLevel) + tableIndex |= 32; + if (m_ptScalarField[(z + 1) * nPointsInSlice + (y + 1) * nPointsInXDirection + (x + 1)] < m_tIsoLevel) + tableIndex |= 64; + if (m_ptScalarField[(z + 1) * nPointsInSlice + y * nPointsInXDirection + (x + 1)] < m_tIsoLevel) + tableIndex |= 128; + + // Now create a triangulation of the isosurface in this + // cell. + if (m_edgeTable[tableIndex] != 0) { + if (m_edgeTable[tableIndex] & 8) { + POINT3DID pt = CalculateIntersection(x, y, z, 3); + unsigned int id = GetEdgeID(x, y, z, 3); + m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt)); + } + if (m_edgeTable[tableIndex] & 1) { + POINT3DID pt = CalculateIntersection(x, y, z, 0); + unsigned int id = GetEdgeID(x, y, z, 0); + m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt)); + } + if (m_edgeTable[tableIndex] & 256) { + POINT3DID pt = CalculateIntersection(x, y, z, 8); + unsigned int id = GetEdgeID(x, y, z, 8); + m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt)); + } + + if (x == m_nCellsX - 1) { + if (m_edgeTable[tableIndex] & 4) { + POINT3DID pt = CalculateIntersection(x, y, z, 2); + unsigned int id = GetEdgeID(x, y, z, 2); + m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt)); + } + if (m_edgeTable[tableIndex] & 2048) { + POINT3DID pt = CalculateIntersection(x, y, z, 11); + unsigned int id = GetEdgeID(x, y, z, 11); + m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt)); + } + } + if (y == m_nCellsY - 1) { + if (m_edgeTable[tableIndex] & 2) { + POINT3DID pt = CalculateIntersection(x, y, z, 1); + unsigned int id = GetEdgeID(x, y, z, 1); + m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt)); + } + if (m_edgeTable[tableIndex] & 512) { + POINT3DID pt = CalculateIntersection(x, y, z, 9); + unsigned int id = GetEdgeID(x, y, z, 9); + m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt)); + } + } + if (z == m_nCellsZ - 1) { + if (m_edgeTable[tableIndex] & 16) { + POINT3DID pt = CalculateIntersection(x, y, z, 4); + unsigned int id = GetEdgeID(x, y, z, 4); + m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt)); + } + if (m_edgeTable[tableIndex] & 128) { + POINT3DID pt = CalculateIntersection(x, y, z, 7); + unsigned int id = GetEdgeID(x, y, z, 7); + m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt)); + } + } + if ((x == m_nCellsX - 1) && (y == m_nCellsY - 1)) + if (m_edgeTable[tableIndex] & 1024) { + POINT3DID pt = CalculateIntersection(x, y, z, 10); + unsigned int id = GetEdgeID(x, y, z, 10); + m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt)); + } + if ((x == m_nCellsX - 1) && (z == m_nCellsZ - 1)) + if (m_edgeTable[tableIndex] & 64) { + POINT3DID pt = CalculateIntersection(x, y, z, 6); + unsigned int id = GetEdgeID(x, y, z, 6); + m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt)); + } + if ((y == m_nCellsY - 1) && (z == m_nCellsZ - 1)) + if (m_edgeTable[tableIndex] & 32) { + POINT3DID pt = CalculateIntersection(x, y, z, 5); + unsigned int id = GetEdgeID(x, y, z, 5); + m_i2pt3idVertices.insert(ID2POINT3DID::value_type(id, pt)); + } + + for (unsigned int i = 0; m_triTable[tableIndex][i] != -1; i += 3) { + TRIANGLE triangle; + unsigned int pointID0, pointID1, pointID2; + pointID0 = GetEdgeID(x, y, z, m_triTable[tableIndex][i]); + pointID1 = GetEdgeID(x, y, z, m_triTable[tableIndex][i + 1]); + pointID2 = GetEdgeID(x, y, z, m_triTable[tableIndex][i + 2]); + triangle.pointID[0] = pointID0; + triangle.pointID[1] = pointID1; + triangle.pointID[2] = pointID2; + m_trivecTriangles.push_back(triangle); + } + } + } + + RenameVerticesAndTriangles(); + CalculateNormals(); + m_bValidSurface = true; +} + +template +bool CIsoSurface::IsSurfaceValid() { + return m_bValidSurface; +} + +template +void CIsoSurface::DeleteSurface() { + m_fCellLengthX = 0; + m_fCellLengthY = 0; + m_fCellLengthZ = 0; + m_nCellsX = 0; + m_nCellsY = 0; + m_nCellsZ = 0; + m_nTriangles = 0; + m_nNormals = 0; + m_nVertices = 0; + if (m_ppt3dVertices != NULL) { + delete[] m_ppt3dVertices; + m_ppt3dVertices = NULL; + } + if (m_piTriangleIndices != NULL) { + delete[] m_piTriangleIndices; + m_piTriangleIndices = NULL; + } + if (m_pvec3dNormals != NULL) { + delete[] m_pvec3dNormals; + m_pvec3dNormals = NULL; + } + m_ptScalarField = NULL; + m_tIsoLevel = 0; + m_bValidSurface = false; +} + +template +int CIsoSurface::GetVolumeLengths(float& fVolLengthX, float& fVolLengthY, float& fVolLengthZ) { + if (IsSurfaceValid()) { + fVolLengthX = m_fCellLengthX * m_nCellsX; + fVolLengthY = m_fCellLengthY * m_nCellsY; + fVolLengthZ = m_fCellLengthZ * m_nCellsZ; + return 1; + } else + return -1; +} + +template +unsigned int CIsoSurface::GetEdgeID(unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nEdgeNo) { + switch (nEdgeNo) { + case 0: + return GetVertexID(nX, nY, nZ) + 1; + case 1: + return GetVertexID(nX, nY + 1, nZ); + case 2: + return GetVertexID(nX + 1, nY, nZ) + 1; + case 3: + return GetVertexID(nX, nY, nZ); + case 4: + return GetVertexID(nX, nY, nZ + 1) + 1; + case 5: + return GetVertexID(nX, nY + 1, nZ + 1); + case 6: + return GetVertexID(nX + 1, nY, nZ + 1) + 1; + case 7: + return GetVertexID(nX, nY, nZ + 1); + case 8: + return GetVertexID(nX, nY, nZ) + 2; + case 9: + return GetVertexID(nX, nY + 1, nZ) + 2; + case 10: + return GetVertexID(nX + 1, nY + 1, nZ) + 2; + case 11: + return GetVertexID(nX + 1, nY, nZ) + 2; + default: + // Invalid edge no. + return -1; + } +} + +template +unsigned int CIsoSurface::GetVertexID(unsigned int nX, unsigned int nY, unsigned int nZ) { + return 3 * (nZ * (m_nCellsY + 1) * (m_nCellsX + 1) + nY * (m_nCellsX + 1) + nX); +} + +template +POINT3DID CIsoSurface::CalculateIntersection(unsigned int nX, unsigned int nY, unsigned int nZ, + unsigned int nEdgeNo) { + float x1, y1, z1, x2, y2, z2; + unsigned int v1x = nX, v1y = nY, v1z = nZ; + unsigned int v2x = nX, v2y = nY, v2z = nZ; + + switch (nEdgeNo) { + case 0: + v2y += 1; + break; + case 1: + v1y += 1; + v2x += 1; + v2y += 1; + break; + case 2: + v1x += 1; + v1y += 1; + v2x += 1; + break; + case 3: + v1x += 1; + break; + case 4: + v1z += 1; + v2y += 1; + v2z += 1; + break; + case 5: + v1y += 1; + v1z += 1; + v2x += 1; + v2y += 1; + v2z += 1; + break; + case 6: + v1x += 1; + v1y += 1; + v1z += 1; + v2x += 1; + v2z += 1; + break; + case 7: + v1x += 1; + v1z += 1; + v2z += 1; + break; + case 8: + v2z += 1; + break; + case 9: + v1y += 1; + v2y += 1; + v2z += 1; + break; + case 10: + v1x += 1; + v1y += 1; + v2x += 1; + v2y += 1; + v2z += 1; + break; + case 11: + v1x += 1; + v2x += 1; + v2z += 1; + break; + } + + x1 = v1x * m_fCellLengthX; + y1 = v1y * m_fCellLengthY; + z1 = v1z * m_fCellLengthZ; + x2 = v2x * m_fCellLengthX; + y2 = v2y * m_fCellLengthY; + z2 = v2z * m_fCellLengthZ; + + unsigned int nPointsInXDirection = (m_nCellsX + 1); + unsigned int nPointsInSlice = nPointsInXDirection * (m_nCellsY + 1); + T val1 = m_ptScalarField[v1z * nPointsInSlice + v1y * nPointsInXDirection + v1x]; + T val2 = m_ptScalarField[v2z * nPointsInSlice + v2y * nPointsInXDirection + v2x]; + POINT3DID intersection = Interpolate(x1, y1, z1, x2, y2, z2, val1, val2); + + return intersection; +} + +template +POINT3DID CIsoSurface::Interpolate(float fX1, float fY1, float fZ1, float fX2, float fY2, float fZ2, T tVal1, + T tVal2) { + POINT3DID interpolation; + float mu; + + mu = float((m_tIsoLevel - tVal1)) / (tVal2 - tVal1); + interpolation.x = fX1 + mu * (fX2 - fX1); + interpolation.y = fY1 + mu * (fY2 - fY1); + interpolation.z = fZ1 + mu * (fZ2 - fZ1); + + return interpolation; +} + +template +void CIsoSurface::RenameVerticesAndTriangles() { + unsigned int nextID = 0; + ID2POINT3DID::iterator mapIterator = m_i2pt3idVertices.begin(); + TRIANGLEVECTOR::iterator vecIterator = m_trivecTriangles.begin(); + + // Rename vertices. + while (mapIterator != m_i2pt3idVertices.end()) { + (*mapIterator).second.newID = nextID; + nextID++; + mapIterator++; + } + + // Now rename triangles. + while (vecIterator != m_trivecTriangles.end()) { + for (unsigned int i = 0; i < 3; i++) { + unsigned int newID = m_i2pt3idVertices[(*vecIterator).pointID[i]].newID; + (*vecIterator).pointID[i] = newID; + } + vecIterator++; + } + + // Copy all the vertices and triangles into two arrays so that they + // can be efficiently accessed. + // Copy vertices. + mapIterator = m_i2pt3idVertices.begin(); + m_nVertices = m_i2pt3idVertices.size(); + m_ppt3dVertices = new POINT3D[m_nVertices]; + for (unsigned int i = 0; i < m_nVertices; i++, mapIterator++) { + m_ppt3dVertices[i][0] = (*mapIterator).second.x; + m_ppt3dVertices[i][1] = (*mapIterator).second.y; + m_ppt3dVertices[i][2] = (*mapIterator).second.z; + } + // Copy vertex indices which make triangles. + vecIterator = m_trivecTriangles.begin(); + m_nTriangles = m_trivecTriangles.size(); + m_piTriangleIndices = new unsigned int[m_nTriangles * 3]; + for (unsigned int i = 0; i < m_nTriangles; i++, vecIterator++) { + m_piTriangleIndices[i * 3] = (*vecIterator).pointID[0]; + m_piTriangleIndices[i * 3 + 1] = (*vecIterator).pointID[1]; + m_piTriangleIndices[i * 3 + 2] = (*vecIterator).pointID[2]; + } + + m_i2pt3idVertices.clear(); + m_trivecTriangles.clear(); +} + +template +void CIsoSurface::CalculateNormals() { + m_nNormals = m_nVertices; + m_pvec3dNormals = new VECTOR3D[m_nNormals]; + + // Set all normals to 0. + for (unsigned int i = 0; i < m_nNormals; i++) { + m_pvec3dNormals[i][0] = 0; + m_pvec3dNormals[i][1] = 0; + m_pvec3dNormals[i][2] = 0; + } + + // Calculate normals. + for (unsigned int i = 0; i < m_nTriangles; i++) { + VECTOR3D vec1, vec2, normal; + unsigned int id0, id1, id2; + id0 = m_piTriangleIndices[i * 3]; + id1 = m_piTriangleIndices[i * 3 + 1]; + id2 = m_piTriangleIndices[i * 3 + 2]; + vec1[0] = m_ppt3dVertices[id1][0] - m_ppt3dVertices[id0][0]; + vec1[1] = m_ppt3dVertices[id1][1] - m_ppt3dVertices[id0][1]; + vec1[2] = m_ppt3dVertices[id1][2] - m_ppt3dVertices[id0][2]; + vec2[0] = m_ppt3dVertices[id2][0] - m_ppt3dVertices[id0][0]; + vec2[1] = m_ppt3dVertices[id2][1] - m_ppt3dVertices[id0][1]; + vec2[2] = m_ppt3dVertices[id2][2] - m_ppt3dVertices[id0][2]; + normal[0] = vec1[2] * vec2[1] - vec1[1] * vec2[2]; + normal[1] = vec1[0] * vec2[2] - vec1[2] * vec2[0]; + normal[2] = vec1[1] * vec2[0] - vec1[0] * vec2[1]; + m_pvec3dNormals[id0][0] += normal[0]; + m_pvec3dNormals[id0][1] += normal[1]; + m_pvec3dNormals[id0][2] += normal[2]; + m_pvec3dNormals[id1][0] += normal[0]; + m_pvec3dNormals[id1][1] += normal[1]; + m_pvec3dNormals[id1][2] += normal[2]; + m_pvec3dNormals[id2][0] += normal[0]; + m_pvec3dNormals[id2][1] += normal[1]; + m_pvec3dNormals[id2][2] += normal[2]; + } + + // Normalize normals. + for (unsigned int i = 0; i < m_nNormals; i++) { + float length = sqrt(m_pvec3dNormals[i][0] * m_pvec3dNormals[i][0] + m_pvec3dNormals[i][1] * m_pvec3dNormals[i][1] + + m_pvec3dNormals[i][2] * m_pvec3dNormals[i][2]); + m_pvec3dNormals[i][0] /= length; + m_pvec3dNormals[i][1] /= length; + m_pvec3dNormals[i][2] /= length; + } +} + +template class CIsoSurface; +template class CIsoSurface; +template class CIsoSurface; + +} // namespace marchingcubes +} // namespace polyscope + + +#endif // CISOSURFACE_H diff --git a/include/polyscope/marchingcubes/Vectors.h b/include/polyscope/marchingcubes/Vectors.h new file mode 100644 index 00000000..00478910 --- /dev/null +++ b/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 diff --git a/include/polyscope/marchingcubes/mesh_implicit_surface.h b/include/polyscope/marchingcubes/mesh_implicit_surface.h new file mode 100644 index 00000000..35d07e2a --- /dev/null +++ b/include/polyscope/marchingcubes/mesh_implicit_surface.h @@ -0,0 +1,69 @@ +#pragma once + +#include +#include +#include + +#include "CIsoSurface.h" + +namespace polyscope { +namespace marchingcubes { + +template +void SampleFunctionToGrid(const Implicit &funct, size_t numCornersPerSide, glm::vec3 center, double sideLength, std::vector &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 &field, double isoLevel, size_t numCornersPerSide, glm::vec3 center, double sideLength, + std::vector &nodes, std::vector> &triangles) { + CIsoSurface* iso = new CIsoSurface(); + 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 \ No newline at end of file diff --git a/include/polyscope/volumetric_grid.h b/include/polyscope/volumetric_grid.h new file mode 100644 index 00000000..cc673324 --- /dev/null +++ b/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 + +#include "marchingcubes/mesh_implicit_surface.h" + +namespace polyscope { + +class VolumetricGrid; +class VolumetricGridScalarIsosurface; +class VolumetricGridScalarQuantity; +class VolumetricGridVectorQuantity; + +template <> // Specialize the quantity type +struct QuantityTypeHelper { + typedef VolumetricGridQuantity type; +}; + + +class VolumetricGrid : public QuantityStructure { +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 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 + VolumetricGridScalarIsosurface* addGridIsosurfaceQuantity(std::string name, double isoLevel, const T& values) { + validateSize(values, nValues(), "grid isosurface quantity " + name); + return addIsosurfaceQuantityImpl(name, isoLevel, standardizeArray(values)); + } + + template + VolumetricGridScalarIsosurface* addGridIsosurfaceQuantityFromFunction(std::string name, double isoLevel, const Funct& funct) { + size_t totalValues = nCornersPerSide * nCornersPerSide * nCornersPerSide; + std::vector field(totalValues); + marchingcubes::SampleFunctionToGrid(funct, nCornersPerSide, gridCenter, sideLength, field); + return addGridIsosurfaceQuantity(name, isoLevel, field); + } + + template + VolumetricGridScalarQuantity* addGridScalarQuantity(std::string name, const T& values, DataType dataType_) { + validateSize(values, nValues(), "grid scalar quantity " + name); + return addScalarQuantityImpl(name, standardizeArray(values), dataType_); + } + + template + VolumetricGridScalarQuantity* addGridScalarQuantityFromFunction(std::string name, const Funct& funct, DataType dataType_) { + size_t totalValues = nCornersPerSide * nCornersPerSide * nCornersPerSide; + std::vector field(totalValues); + marchingcubes::SampleFunctionToGrid(funct, nCornersPerSide, gridCenter, sideLength, field); + return addGridScalarQuantity(name, field, dataType_); + } + + template + VolumetricGridVectorQuantity* addGridVectorQuantity(std::string name, const T& vecValues, VectorType dataType_) { + validateSize(vecValues, nValues(), "grid vector quantity " + name); + return addVectorQuantityImpl(name, standardizeArray(vecValues), dataType_); + } + + template + VolumetricGridVectorQuantity* addGridVectorQuantityFromFunction(std::string name, const Funct& funct, VectorType dataType_) { + size_t totalValues = nCornersPerSide * nCornersPerSide * nCornersPerSide; + std::vector field(totalValues); + marchingcubes::SampleFunctionToGrid(funct, nCornersPerSide, gridCenter, sideLength, field); + return addGridVectorQuantity(name, field, dataType_); + } + + + glm::vec3 getColor(); + +private: + PersistentValue color; + VolumetricGrid* setColor(glm::vec3 newVal); + VolumetricGridScalarIsosurface* addIsosurfaceQuantityImpl(std::string name, double isoLevel, + const std::vector& data); + VolumetricGridScalarQuantity* addScalarQuantityImpl(std::string name, const std::vector& data, + DataType dataType_); + VolumetricGridVectorQuantity* addVectorQuantityImpl(std::string name, const std::vector& data, + VectorType dataType_); +}; + + +VolumetricGrid* registerVolumetricGrid(std::string name, size_t nValuesPerSide, glm::vec3 center, double sideLen); + +template +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 diff --git a/include/polyscope/volumetric_grid_quantity.h b/include/polyscope/volumetric_grid_quantity.h new file mode 100644 index 00000000..19c9de2b --- /dev/null +++ b/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 +class VolumetricGridQuantity : public Quantity { +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 diff --git a/include/polyscope/volumetric_grid_scalar_isosurface.h b/include/polyscope/volumetric_grid_scalar_isosurface.h new file mode 100644 index 00000000..cc256cab --- /dev/null +++ b/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 &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 values; + double levelSet; + double newLevelSet; + double increment; + void prepare(); + void prepareTriangleIndices(); + void setProgramUniforms(gl::GLProgram& program); + void fillGeometryBuffersSmooth(gl::GLProgram& p); + std::unique_ptr meshProgram; + std::vector nodes; + std::vector normals; + std::vector> triangles; +}; + +} // namespace polyscope diff --git a/include/polyscope/volumetric_grid_scalar_quantity.h b/include/polyscope/volumetric_grid_scalar_quantity.h new file mode 100644 index 00000000..2369ebf4 --- /dev/null +++ b/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& 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 values; + std::vector positions; + PersistentValue cMap; + std::unique_ptr pointsProgram; + float pointRadius; + std::pair vizRange; + std::pair dataRange; +}; + +} // namespace polyscope diff --git a/include/polyscope/volumetric_grid_vector_quantity.h b/include/polyscope/volumetric_grid_vector_quantity.h new file mode 100644 index 00000000..03e6702f --- /dev/null +++ b/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& 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 vectorValues; + std::vector positions; + std::unique_ptr vectorsProgram; + AffineRemapper mapper; + float vectorRadius; + float vectorLengthMult; + + void fillPositions(); + void createProgram(); +}; +} // namespace polyscope \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a55303c4..1f439b45 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -175,12 +175,20 @@ SET(SRCS volume_mesh_color_quantity.cpp volume_mesh_scalar_quantity.cpp volume_mesh_vector_quantity.cpp + + # Volume grid + volumetric_grid.cpp + volumetric_grid_scalar_quantity.cpp + volumetric_grid_vector_quantity.cpp + volumetric_grid_scalar_isosurface.cpp # Rendering utilities vector_artist.cpp trace_vector_field.cpp ribbon_artist.cpp image_scalar_artist.cpp + marchingcubes/Vectors.cpp + implicit_surface.cpp ## Embedded binary data render/bindata/bindata_font_lato_regular.cpp diff --git a/src/marchingcubes/CIsoSurface.cpp b/src/marchingcubes/CIsoSurface.cpp new file mode 100644 index 00000000..c9f63ad8 --- /dev/null +++ b/src/marchingcubes/CIsoSurface.cpp @@ -0,0 +1,13 @@ +// File Name: CIsoSurface.cpp +// Last Modified: 5/8/2000 +// Author: Raghavendra Chandrashekara (based on source code provided +// by Paul Bourke and Cory Gene Bloyd) +// Email: rc99@doc.ic.ac.uk, rchandrashekara@hotmail.com +// +// Description: This is the implementation file for the CIsoSurface class. + +// Edited and ported to polyscope by Christopher Yu. + +#include +#include "polyscope/marchingcubes/CIsoSurface.h" + diff --git a/src/marchingcubes/Vectors.cpp b/src/marchingcubes/Vectors.cpp new file mode 100644 index 00000000..bb255b7b --- /dev/null +++ b/src/marchingcubes/Vectors.cpp @@ -0,0 +1,99 @@ +// File Name: Vectors.h +// Last Modified: 9/8/2000 +// Author: Raghavendra Chandrashekara +// Email: rc99@doc.ic.ac.uk, rchandrashekara@hotmail.com +// +// Description: This is the implementation file for POINT3DXYZ class. + +// Edited and ported to polyscope by Christopher Yu. + +#include "polyscope/marchingcubes/Vectors.h" + +namespace polyscope { +namespace marchingcubes { + + +POINT3DXYZ operator+(const POINT3DXYZ& pt3dPoint1, const POINT3DXYZ& pt3dPoint2) { + POINT3DXYZ result; + + result.x = pt3dPoint1.x + pt3dPoint2.x; + result.y = pt3dPoint1.y + pt3dPoint2.y; + result.z = pt3dPoint1.z + pt3dPoint2.z; + + return result; +} + +POINT3DXYZ operator-(const POINT3DXYZ& pt3dPoint1, const POINT3DXYZ& pt3dPoint2) { + POINT3DXYZ result; + + result.x = pt3dPoint1.x - pt3dPoint2.x; + result.y = pt3dPoint1.y - pt3dPoint2.y; + result.z = pt3dPoint1.z - pt3dPoint2.z; + + return result; +} + +POINT3DXYZ operator*(const POINT3DXYZ& pt3dPoint, float fScale) { + POINT3DXYZ result; + + result.x = pt3dPoint.x * fScale; + result.y = pt3dPoint.y * fScale; + result.z = pt3dPoint.z * fScale; + + return result; +} + +POINT3DXYZ operator*(float fScale, const POINT3DXYZ& pt3dPoint) { + POINT3DXYZ result; + + result.x = pt3dPoint.x * fScale; + result.y = pt3dPoint.y * fScale; + result.z = pt3dPoint.z * fScale; + + return result; +} + +POINT3DXYZ operator/(const POINT3DXYZ& pt3dPoint, float fScale) { + POINT3DXYZ result; + + result.x = pt3dPoint.x / fScale; + result.y = pt3dPoint.y / fScale; + result.z = pt3dPoint.z / fScale; + + return result; +} + +POINT3DXYZ& operator*=(POINT3DXYZ& pt3dPoint, float fScale) { + pt3dPoint.x *= fScale; + pt3dPoint.y *= fScale; + pt3dPoint.z *= fScale; + + return pt3dPoint; +} + +POINT3DXYZ& operator/=(POINT3DXYZ& pt3dPoint, float fScale) { + pt3dPoint.x /= fScale; + pt3dPoint.y /= fScale; + pt3dPoint.z /= fScale; + + return pt3dPoint; +} + +POINT3DXYZ& operator+=(POINT3DXYZ& pt3dPoint1, const POINT3DXYZ& pt3dPoint2) { + pt3dPoint1.x += pt3dPoint2.x; + pt3dPoint1.y += pt3dPoint2.y; + pt3dPoint1.z += pt3dPoint2.z; + + return pt3dPoint1; +} + +POINT3DXYZ& operator-=(POINT3DXYZ& pt3dPoint1, const POINT3DXYZ& pt3dPoint2) { + pt3dPoint1.x -= pt3dPoint2.x; + pt3dPoint1.y -= pt3dPoint2.y; + pt3dPoint1.z -= pt3dPoint2.z; + + return pt3dPoint1; +} + +} // namespace marchingcubes +} // namespace polyscope diff --git a/src/volumetric_grid.cpp b/src/volumetric_grid.cpp new file mode 100644 index 00000000..6c49439c --- /dev/null +++ b/src/volumetric_grid.cpp @@ -0,0 +1,118 @@ +// Copyright 2017-2020, Christopher Yu, Nicholas Sharp and the +// Polyscope contributors. http://polyscope.run. + +#include "polyscope/volumetric_grid.h" +#include "polyscope/volumetric_grid_scalar_isosurface.h" +#include "polyscope/volumetric_grid_scalar_quantity.h" +#include "polyscope/volumetric_grid_vector_quantity.h" + +#include "polyscope/combining_hash_functions.h" +#include "polyscope/gl/gl_utils.h" +#include "polyscope/gl/materials/materials.h" +#include "polyscope/gl/shaders.h" +#include "polyscope/gl/shaders/surface_shaders.h" +#include "polyscope/gl/shaders/wireframe_shaders.h" +#include "polyscope/pick.h" +#include "polyscope/polyscope.h" + +#include "polyscope/surface_mesh.h" + +#include "imgui.h" + +namespace polyscope { +// Initialize statics +const std::string VolumetricGrid::structureTypeName = "Implicit Surface"; + +VolumetricGrid::VolumetricGrid(std::string name, size_t nValuesPerSide, glm::vec3 center, double sideLen) + : QuantityStructure(name, typeName()), color(uniquePrefix() + "#color", getNextUniqueColor()) { + nCornersPerSide = nValuesPerSide; + sideLength = sideLen; + gridCenter = center; +} + +VolumetricGrid* VolumetricGrid::setColor(glm::vec3 newVal) { + color = newVal; + polyscope::requestRedraw(); + return this; +} + +glm::vec3 VolumetricGrid::getColor() { return color.get(); } + +void VolumetricGrid::buildCustomUI() { + ImGui::Text("samples: %lld (%lld per side)", static_cast(nValues()), + static_cast(nCornersPerSide)); + ImGui::Text("center: (%.3f, %.3f, %.3f)", gridCenter.x, gridCenter.y, gridCenter.z); + ImGui::Text("grid side length: %.4f", sideLength); + if (ImGui::ColorEdit3("Color", &color.get()[0], ImGuiColorEditFlags_NoInputs)) { + setColor(getColor()); + } +} + +void VolumetricGrid::buildPickUI(size_t localPickID) { + // For now do nothing +} + +void VolumetricGrid::draw() { + // For now, do nothing for the actual grid + + // Draw the quantities + for (auto& x : quantities) { + x.second->draw(); + } +} + +void VolumetricGrid::drawPick() { + // For now do nothing +} + +double VolumetricGrid::lengthScale() { return sideLength; } + +std::tuple VolumetricGrid::boundingBox() { + double r = sideLength / 2; + glm::vec3 minCorner = gridCenter - glm::vec3{r, r, r}; + glm::vec3 maxCorner = gridCenter + glm::vec3{r, r, r}; + + return std::make_tuple(minCorner, maxCorner); +} + +std::string VolumetricGrid::typeName() { return structureTypeName; } + +VolumetricGridQuantity::VolumetricGridQuantity(std::string name_, VolumetricGrid& curveNetwork_, bool dominates_) + : Quantity(name_, curveNetwork_, dominates_) {} + + +void VolumetricGridQuantity::buildNodeInfoGUI(size_t nodeInd) {} +void VolumetricGridQuantity::buildEdgeInfoGUI(size_t edgeInd) {} + +VolumetricGridScalarIsosurface* VolumetricGrid::addIsosurfaceQuantityImpl(std::string name, double isoLevel, + const std::vector& data) { + VolumetricGridScalarIsosurface* q = new VolumetricGridScalarIsosurface(name, *this, data, isoLevel); + addQuantity(q); + q->setEnabled(true); + return q; +} + +VolumetricGridScalarQuantity* VolumetricGrid::addScalarQuantityImpl(std::string name, const std::vector& data, + DataType dataType_) { + VolumetricGridScalarQuantity* q = new VolumetricGridScalarQuantity(name, *this, data, dataType_); + addQuantity(q); + return q; +} + +VolumetricGridVectorQuantity* VolumetricGrid::addVectorQuantityImpl(std::string name, const std::vector& data, + VectorType dataType_) { + VolumetricGridVectorQuantity* q = new VolumetricGridVectorQuantity(name, *this, data, dataType_); + addQuantity(q); + return q; +} + +VolumetricGrid* registerVolumetricGrid(std::string name, size_t nValuesPerSide, glm::vec3 center, double sideLen) { + VolumetricGrid* s = new VolumetricGrid(name, nValuesPerSide, center, sideLen); + bool success = registerStructure(s); + if (!success) { + safeDelete(s); + } + return s; +} + +} // namespace polyscope \ No newline at end of file diff --git a/src/volumetric_grid_scalar_isosurface.cpp b/src/volumetric_grid_scalar_isosurface.cpp new file mode 100644 index 00000000..a9e63db3 --- /dev/null +++ b/src/volumetric_grid_scalar_isosurface.cpp @@ -0,0 +1,173 @@ +#include "polyscope/volumetric_grid_scalar_isosurface.h" +#include "polyscope/marchingcubes/mesh_implicit_surface.h" + +#include "polyscope/gl/materials/materials.h" +#include "polyscope/gl/shaders.h" +#include "polyscope/gl/shaders/surface_shaders.h" + +namespace polyscope { + +VolumetricGridScalarIsosurface::VolumetricGridScalarIsosurface(std::string name, VolumetricGrid& grid_, + const std::vector& values_, double levelSet_) + : VolumetricGridQuantity(name, grid_, true), values(std::move(values_)) { + levelSet = levelSet_; + meshIsosurface(); + increment = 0.01; + newLevelSet = levelSet; +} + +void VolumetricGridScalarIsosurface::recomputeNormals() { + normals.clear(); + normals.resize(nodes.size()); + + for (size_t i = 0; i < normals.size(); i++) { + normals[i] = glm::vec3{0., 0., 0.}; + } + + for (auto& tri : triangles) { + glm::vec3 v0 = nodes[tri[0]]; + glm::vec3 v1 = nodes[tri[1]]; + glm::vec3 v2 = nodes[tri[2]]; + + glm::vec3 e01 = v1 - v0; + glm::vec3 e02 = v2 - v0; + glm::vec3 areaWtNormal = glm::cross(e01, e02); + + normals[tri[0]] += areaWtNormal; + normals[tri[1]] += areaWtNormal; + normals[tri[2]] += areaWtNormal; + } + + for (size_t i = 0; i < normals.size(); i++) { + normals[i] = glm::normalize(normals[i]); + } +} + +void VolumetricGridScalarIsosurface::meshIsosurface() { + nodes.clear(); + triangles.clear(); + + size_t nCornersPerSide = parent.nCornersPerSide; + glm::vec3 gridCenter = parent.gridCenter; + double sideLength = parent.sideLength; + marchingcubes::MeshImplicitGrid(values, levelSet, nCornersPerSide, gridCenter, sideLength, nodes, triangles); + + std::cout << "Meshed isosurface, producing " << nodes.size() << " nodes, " << triangles.size() << " triangles" + << std::endl; + recomputeNormals(); + + prepare(); +} + +void VolumetricGridScalarIsosurface::draw() { + if (!isEnabled()) return; + + if (meshProgram == nullptr) { + prepare(); + } + + // Set uniforms + parent.setTransformUniforms(*meshProgram); + setProgramUniforms(*meshProgram); + meshProgram->draw(); +} + +void VolumetricGridScalarIsosurface::buildCustomUI() { + ImGui::Text("Current isolevel: %.4f", levelSet); + ImGui::PushItemWidth(100); + ImGui::InputDouble(" ", &newLevelSet); + + ImGui::SameLine(); + + if (ImGui::Button("Set isolevel")) { + levelSet = newLevelSet; + meshIsosurface(); + } + + ImGui::Spacing(); + + ImGui::Text("Increment level:"); + ImGui::SameLine(); + if (ImGui::Button("-", ImVec2{32, 20})) { + levelSet -= increment; + newLevelSet = levelSet; + meshIsosurface(); + } + ImGui::SameLine(); + if (ImGui::Button("+", ImVec2{32, 20})) { + levelSet += increment; + newLevelSet = levelSet; + meshIsosurface(); + } + + ImGui::InputDouble("Step size", &increment); + ImGui::PopItemWidth(); +} + +void VolumetricGridScalarIsosurface::setProgramUniforms(gl::GLProgram& program) { + program.setUniform("u_basecolor", parent.getColor()); +} + +std::string VolumetricGridScalarIsosurface::niceName() { return name + " (isosurface)"; } + +void VolumetricGridScalarIsosurface::geometryChanged() { + meshProgram.reset(); +} + +void VolumetricGridScalarIsosurface::prepareTriangleIndices() { + // Need to upload triangle indices + // Indicies are supplied as size_t but need to be cast to unsigned int + std::vector> index_triangles(triangles.size()); + for (size_t i = 0; i < index_triangles.size(); i++) { + size_t ui_max = (size_t)INT_MAX * 2; + if (triangles[i][0] > ui_max || triangles[i][1] > ui_max || triangles[i][2] > ui_max) { + std::cerr << "Number of triangles exceeds unsigned int max; exiting" << std::endl; + exit(1); + } + index_triangles[i] = {static_cast(triangles[i][0]), static_cast(triangles[i][1]), + static_cast(triangles[i][2])}; + } + + meshProgram->setIndex(index_triangles); +} + +void VolumetricGridScalarIsosurface::prepare() { + // TODO: figure out what shaders to use + meshProgram.reset(new gl::GLProgram(&gl::PLAIN_SURFACE_VERT_SHADER, &gl::PLAIN_SURFACE_FRAG_SHADER, + gl::DrawMode::IndexedTriangles)); + + prepareTriangleIndices(); + + // Populate draw buffers + fillGeometryBuffersSmooth(*meshProgram); + + setMaterialForProgram(*meshProgram, "wax"); +} + + +void VolumetricGridScalarIsosurface::fillGeometryBuffersSmooth(gl::GLProgram& p) { + std::vector bcoord; + size_t nFaces = triangles.size(); + bool wantsBary = p.hasAttribute("a_barycoord"); + + if (wantsBary) { + bcoord.reserve(3 * nFaces); + } + + for (size_t iF = 0; iF < nFaces; iF++) { + if (wantsBary) { + bcoord.push_back(glm::vec3{1., 0., 0.}); + bcoord.push_back(glm::vec3{0., 1., 0.}); + bcoord.push_back(glm::vec3{0., 0., 1.}); + } + } + + // Store data in buffers + p.setAttribute("a_position", nodes); + p.setAttribute("a_normal", normals); + if (wantsBary) { + p.setAttribute("a_barycoord", bcoord); + } +} + +} // namespace polyscope diff --git a/src/volumetric_grid_scalar_quantity.cpp b/src/volumetric_grid_scalar_quantity.cpp new file mode 100644 index 00000000..7f850a1b --- /dev/null +++ b/src/volumetric_grid_scalar_quantity.cpp @@ -0,0 +1,93 @@ +#include "polyscope/volumetric_grid_scalar_quantity.h" + +#include "polyscope/gl/materials/materials.h" +#include "polyscope/gl/shaders.h" +#include "polyscope/gl/shaders/sphere_shaders.h" + + +namespace polyscope { + +VolumetricGridScalarQuantity::VolumetricGridScalarQuantity(std::string name, VolumetricGrid& grid_, + const std::vector& values_, DataType dataType_) + : VolumetricGridQuantity(name, grid_, true), dataType(dataType_), values(std::move(values_)), + cMap(uniquePrefix() + "#cmap", defaultColorMap(dataType)) { + + fillPositions(); + dataRange = robustMinMax(values, 1e-5); + resetMapRange(); + pointRadius = 0.5 * parent.sideLength / (parent.nCornersPerSide - 1); +} + + +void VolumetricGridScalarQuantity::buildCustomUI() { + ImGui::SliderFloat("Point radius", &pointRadius, 0, 1); +} + +std::string VolumetricGridScalarQuantity::niceName() { return name + " (scalar)"; } + +void VolumetricGridScalarQuantity::geometryChanged() { + pointsProgram.reset(); +} + +void VolumetricGridScalarQuantity::fillPositions() { + positions.clear(); + // Fill the positions vector with each grid corner position + size_t nValues = parent.nValues(); + positions.resize(nValues); + for (size_t i = 0; i < nValues; i++) { + positions[i] = parent.positionOfIndex(i); + } +} + +void VolumetricGridScalarQuantity::draw() { + if (!isEnabled()) return; + + if (pointsProgram == nullptr) { + createProgram(); + } + parent.setTransformUniforms(*pointsProgram); + setPointCloudUniforms(); + pointsProgram->setUniform("u_rangeLow", vizRange.first); + pointsProgram->setUniform("u_rangeHigh", vizRange.second); + + pointsProgram->draw(); +} + +void VolumetricGridScalarQuantity::resetMapRange() { + switch (dataType) { + case DataType::STANDARD: + vizRange = dataRange; + break; + case DataType::SYMMETRIC: { + double absRange = std::max(std::abs(dataRange.first), std::abs(dataRange.second)); + vizRange = std::make_pair(-absRange, absRange); + } break; + case DataType::MAGNITUDE: + vizRange = std::make_pair(0., dataRange.second); + break; + } + + requestRedraw(); +} + +void VolumetricGridScalarQuantity::setPointCloudUniforms() { + pointsProgram->setUniform("u_pointRadius", pointRadius); + + glm::vec3 lookDir, upDir, rightDir; + view::getCameraFrame(lookDir, upDir, rightDir); + pointsProgram->setUniform("u_camZ", lookDir); + pointsProgram->setUniform("u_camUp", upDir); + pointsProgram->setUniform("u_camRight", rightDir); +} + +void VolumetricGridScalarQuantity::createProgram() { + pointsProgram.reset(new gl::GLProgram(&gl::SPHERE_VALUE_VERT_SHADER, &gl::SPHERE_VALUE_BILLBOARD_GEOM_SHADER, + &gl::SPHERE_VALUE_BILLBOARD_FRAG_SHADER, gl::DrawMode::Points)); + gl::setMaterialForProgram(*pointsProgram, "wax"); + + pointsProgram->setAttribute("a_position", positions); + pointsProgram->setAttribute("a_value", values); + pointsProgram->setTextureFromColormap("t_colormap", gl::getColorMap(cMap.get())); +} + +} // namespace polyscope \ No newline at end of file diff --git a/src/volumetric_grid_vector_quantity.cpp b/src/volumetric_grid_vector_quantity.cpp new file mode 100644 index 00000000..5704092e --- /dev/null +++ b/src/volumetric_grid_vector_quantity.cpp @@ -0,0 +1,86 @@ +#include "polyscope/volumetric_grid_vector_quantity.h" + +#include "polyscope/gl/materials/materials.h" +#include "polyscope/gl/shaders.h" +#include "polyscope/gl/shaders/vector_shaders.h" + +namespace polyscope { + +VolumetricGridVectorQuantity::VolumetricGridVectorQuantity(std::string name, VolumetricGrid& grid_, + const std::vector& values_, + VectorType vectorType_) + : VolumetricGridQuantity(name, grid_, true), vectorValues(std::move(values_)) { + vectorType = vectorType_; + vectorRadius = (vectorType == VectorType::AMBIENT) ? 1.0f : 0.02f; + vectorLengthMult = 0.2f; + + fillPositions(); + + // Create a mapper (default mapper is identity) + if (vectorType == VectorType::AMBIENT) { + mapper.setMinMax(vectorValues); + } else { + mapper = AffineRemapper(vectorValues, DataType::MAGNITUDE); + } +} + +void VolumetricGridVectorQuantity::fillPositions() { + positions.clear(); + // Fill the positions vector with each grid corner position + size_t nValues = parent.nValues(); + positions.resize(nValues); + for (size_t i = 0; i < nValues; i++) { + positions[i] = parent.positionOfIndex(i); + } +} + +void VolumetricGridVectorQuantity::buildCustomUI() { + ImGui::SliderFloat("Length", &vectorLengthMult, 0.0, 0.5, "%.5f", 3.); + ImGui::SliderFloat("Radius", &vectorRadius, 0.0, 0.5, "%.5f", 3.); + + ImGui::TextUnformatted(mapper.printBounds().c_str()); +} + +std::string VolumetricGridVectorQuantity::niceName() { return name + " (vector)"; } + +void VolumetricGridVectorQuantity::geometryChanged() { vectorsProgram.reset(); } + +void VolumetricGridVectorQuantity::draw() { + if (!isEnabled()) return; + + if (vectorsProgram == nullptr) { + createProgram(); + } + + // Set uniforms + parent.setTransformUniforms(*vectorsProgram); + + vectorsProgram->setUniform("u_radius", vectorRadius); + vectorsProgram->setUniform("u_color", parent.getColor()); + + if (vectorType == VectorType::AMBIENT) { + vectorsProgram->setUniform("u_lengthMult", 1.0); + } else { + vectorsProgram->setUniform("u_lengthMult", vectorLengthMult); + } + + vectorsProgram->draw(); +} + +void VolumetricGridVectorQuantity::createProgram() { + vectorsProgram.reset(new gl::GLProgram(&gl::PASSTHRU_VECTOR_VERT_SHADER, &gl::VECTOR_GEOM_SHADER, + &gl::SHINY_VECTOR_FRAG_SHADER, gl::DrawMode::Points)); + + // Fill buffers + std::vector mappedVectors; + for (glm::vec3 v : vectorValues) { + mappedVectors.push_back(mapper.map(v)); + } + + vectorsProgram->setAttribute("a_vector", mappedVectors); + vectorsProgram->setAttribute("a_position", positions); + + setMaterialForProgram(*vectorsProgram, "wax"); +} + +} // namespace polyscope \ No newline at end of file