diff options
Diffstat (limited to 'thirdparty/thekla_atlas/nvmesh/param')
16 files changed, 5992 insertions, 0 deletions
diff --git a/thirdparty/thekla_atlas/nvmesh/param/Atlas.cpp b/thirdparty/thekla_atlas/nvmesh/param/Atlas.cpp new file mode 100644 index 0000000000..98f92cef96 --- /dev/null +++ b/thirdparty/thekla_atlas/nvmesh/param/Atlas.cpp @@ -0,0 +1,1519 @@ +// Copyright NVIDIA Corporation 2006 -- Ignacio Castano <icastano@nvidia.com> + +#include "nvmesh.h" // pch + +#include "Atlas.h" +#include "Util.h" +#include "AtlasBuilder.h" +#include "AtlasPacker.h" +#include "SingleFaceMap.h" +#include "OrthogonalProjectionMap.h" +#include "LeastSquaresConformalMap.h" +#include "ParameterizationQuality.h" + +//#include "nvmesh/export/MeshExportOBJ.h" + +#include "nvmesh/halfedge/Mesh.h" +#include "nvmesh/halfedge/Face.h" +#include "nvmesh/halfedge/Vertex.h" + +#include "nvmesh/MeshBuilder.h" +#include "nvmesh/MeshTopology.h" +#include "nvmesh/param/Util.h" +#include "nvmesh/geometry/Measurements.h" + +#include "nvmath/Vector.inl" +#include "nvmath/Fitting.h" +#include "nvmath/Box.inl" +#include "nvmath/ProximityGrid.h" +#include "nvmath/Morton.h" + +#include "nvcore/StrLib.h" +#include "nvcore/Array.inl" +#include "nvcore/HashMap.inl" + +using namespace nv; + + +/// Ctor. +Atlas::Atlas() +{ + failed=false; +} + +// Dtor. +Atlas::~Atlas() +{ + deleteAll(m_meshChartsArray); +} + +uint Atlas::chartCount() const +{ + uint count = 0; + foreach(c, m_meshChartsArray) { + count += m_meshChartsArray[c]->chartCount(); + } + return count; +} + +const Chart * Atlas::chartAt(uint i) const +{ + foreach(c, m_meshChartsArray) { + uint count = m_meshChartsArray[c]->chartCount(); + + if (i < count) { + return m_meshChartsArray[c]->chartAt(i); + } + + i -= count; + } + + return NULL; +} + +Chart * Atlas::chartAt(uint i) +{ + foreach(c, m_meshChartsArray) { + uint count = m_meshChartsArray[c]->chartCount(); + + if (i < count) { + return m_meshChartsArray[c]->chartAt(i); + } + + i -= count; + } + + return NULL; +} + +// Extract the charts and add to this atlas. +void Atlas::addMeshCharts(MeshCharts * meshCharts) +{ + m_meshChartsArray.append(meshCharts); +} + +void Atlas::extractCharts(const HalfEdge::Mesh * mesh) +{ + MeshCharts * meshCharts = new MeshCharts(mesh); + meshCharts->extractCharts(); + addMeshCharts(meshCharts); +} + +void Atlas::computeCharts(const HalfEdge::Mesh * mesh, const SegmentationSettings & settings, const Array<uint> & unchartedMaterialArray) +{ + failed=false; + MeshCharts * meshCharts = new MeshCharts(mesh); + meshCharts->computeCharts(settings, unchartedMaterialArray); + addMeshCharts(meshCharts); +} + + + + +#if 0 + +/// Compute a seamless texture atlas. +bool Atlas::computeSeamlessTextureAtlas(bool groupFaces/*= true*/, bool scaleTiles/*= false*/, uint w/*= 1024*/, uint h/* = 1024*/) +{ + // Implement seamless texture atlas similar to what ZBrush does. See also: + // "Meshed Atlases for Real-Time Procedural Solid Texturing" + // http://graphics.cs.uiuc.edu/~jch/papers/rtpst.pdf + + // Other methods that we should experiment with: + // + // Seamless Texture Atlases: + // http://www.cs.jhu.edu/~bpurnomo/STA/index.html + // + // Rectangular Multi-Chart Geometry Images: + // http://graphics.cs.uiuc.edu/~jch/papers/rmcgi.pdf + // + // Discrete differential geometry also provide a way of constructing + // seamless quadrangulations as shown in: + // http://www.geometry.caltech.edu/pubs/TACD06.pdf + // + +#pragma message(NV_FILE_LINE "TODO: Implement seamless texture atlas.") + + if (groupFaces) + { + // @@ TODO. + } + else + { + // @@ Create one atlas per face. + } + + if (scaleTiles) + { + // @@ TODO + } + + /* + if (!isQuadMesh(m_mesh)) { + // Only handle quads for now. + return false; + } + + // Each face is a chart. + const uint faceCount = m_mesh->faceCount(); + m_chartArray.resize(faceCount); + + for(uint f = 0; f < faceCount; f++) { + m_chartArray[f].faceArray.clear(); + m_chartArray[f].faceArray.append(f); + } + + // Map each face to a separate square. + + // Determine face layout according to width and height. + float aspect = float(m_width) / float(m_height); + + uint i = 2; + uint total = (m_width / (i+1)) * (m_height / (i+1)); + while(total > faceCount) { + i *= 2; + total = (m_width / (i+1)) * (m_height / (i+1)); + } + + uint tileSize = i / 2; + + int x = 0; + int y = 0; + + m_result = new HalfEdge::Mesh(); + + // Once you have that it's just matter of traversing the faces. + for(uint f = 0; f < faceCount; f++) { + // Compute texture coordinates. + Vector2 tex[4]; + tex[0] = Vector2(float(x), float(y)); + tex[1] = Vector2(float(x+tileSize), float(y)); + tex[2] = Vector2(float(x+tileSize), float(y+tileSize)); + tex[3] = Vector2(float(x), float(y+tileSize)); + + Array<uint> indexArray(4); + + const HalfEdge::Face * face = m_mesh->faceAt(f); + + int i = 0; + for(HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance(), i++) { + const HalfEdge::Edge * edge = it.current(); + const HalfEdge::Vertex * vertex = edge->from(); + + HalfEdge::Vertex * newVertex = m_result->addVertex(vertex->id(), vertex->pos()); + + newVertex->setTex(Vector3(tex[i], 0)); + newVertex->setNor(vertex->nor()); + + indexArray.append(m_result->vertexCount() + 1); + } + + m_result->addFace(indexArray); + + // Move to the next tile. + x += tileSize + 1; + if (x + tileSize > m_width) { + x = 0; + y += tileSize + 1; + } + } + */ + + return false; +} + +#endif + + +void Atlas::parameterizeCharts() +{ + foreach(i, m_meshChartsArray) { + m_meshChartsArray[i]->parameterizeCharts(); + } +} + + +float Atlas::packCharts(int quality, float texelsPerUnit, bool blockAlign, bool conservative) +{ + AtlasPacker packer(this); + packer.packCharts(quality, texelsPerUnit, blockAlign, conservative); + if (hasFailed()) + return 0; + return packer.computeAtlasUtilization(); +} + + + + +/// Ctor. +MeshCharts::MeshCharts(const HalfEdge::Mesh * mesh) : m_mesh(mesh) +{ +} + +// Dtor. +MeshCharts::~MeshCharts() +{ + deleteAll(m_chartArray); +} + + +void MeshCharts::extractCharts() +{ + const uint faceCount = m_mesh->faceCount(); + + int first = 0; + Array<uint> queue(faceCount); + + BitArray bitFlags(faceCount); + bitFlags.clearAll(); + + for (uint f = 0; f < faceCount; f++) + { + if (bitFlags.bitAt(f) == false) + { + // Start new patch. Reset queue. + first = 0; + queue.clear(); + queue.append(f); + bitFlags.setBitAt(f); + + while (first != queue.count()) + { + const HalfEdge::Face * face = m_mesh->faceAt(queue[first]); + + // Visit face neighbors of queue[first] + for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + const HalfEdge::Edge * edge = it.current(); + nvDebugCheck(edge->pair != NULL); + + if (!edge->isBoundary() && /*!edge->isSeam()*/ + //!(edge->from()->tex() != edge->pair()->to()->tex() || edge->to()->tex() != edge->pair()->from()->tex())) + !(edge->from() != edge->pair->to() || edge->to() != edge->pair->from())) // Preserve existing seams (not just texture seams). + { + const HalfEdge::Face * neighborFace = edge->pair->face; + nvDebugCheck(neighborFace != NULL); + + if (bitFlags.bitAt(neighborFace->id) == false) + { + queue.append(neighborFace->id); + bitFlags.setBitAt(neighborFace->id); + } + } + } + + first++; + } + + Chart * chart = new Chart(); + chart->build(m_mesh, queue); + + m_chartArray.append(chart); + } + } +} + + +/* +LSCM: +- identify sharp features using local dihedral angles. +- identify seed faces farthest from sharp features. +- grow charts from these seeds. + +MCGIM: +- phase 1: chart growth + - grow all charts simultaneously using dijkstra search on the dual graph of the mesh. + - graph edges are weighted based on planarity metric. + - metric uses distance to global chart normal. + - terminate when all faces have been assigned. +- phase 2: seed computation: + - place new seed of the chart at the most interior face. + - most interior is evaluated using distance metric only. + +- method repeates the two phases, until the location of the seeds does not change. + - cycles are detected by recording all the previous seeds and chartification terminates. + +D-Charts: + +- Uniaxial conic metric: + - N_c = axis of the generalized cone that best fits the chart. (cone can a be cylinder or a plane). + - omega_c = angle between the face normals and the axis. + - Fitting error between chart C and tringle t: F(c,t) = (N_c*n_t - cos(omega_c))^2 + +- Compactness metrics: + - Roundness: + - C(c,t) = pi * D(S_c,t)^2 / A_c + - S_c = chart seed. + - D(S_c,t) = length of the shortest path inside the chart betwen S_c and t. + - A_c = chart area. + - Straightness: + - P(c,t) = l_out(c,t) / l_in(c,t) + - l_out(c,t) = lenght of the edges not shared between C and t. + - l_in(c,t) = lenght of the edges shared between C and t. + +- Combined metric: + - Cost(c,t) = F(c,t)^alpha + C(c,t)^beta + P(c,t)^gamma + - alpha = 1, beta = 0.7, gamma = 0.5 + + + + +Our basic approach: +- Just one iteration of k-means? +- Avoid dijkstra by greedily growing charts until a threshold is met. Increase threshold and repeat until no faces left. +- If distortion metric is too high, split chart, add two seeds. +- If chart size is low, try removing chart. + + +Postprocess: +- If topology is not disk: + - Fill holes, if new faces fit proxy. + - Find best cut, otherwise. +- After parameterization: + - If boundary self-intersects: + - cut chart along the closest two diametral boundary vertices, repeat parametrization. + - what if the overlap is on an appendix? How do we find that out and cut appropiately? + - emphasize roundness metrics to prevent those cases. + - If interior self-overlaps: preserve boundary parameterization and use mean-value map. + +*/ + + +SegmentationSettings::SegmentationSettings() +{ + // Charts have no area or boundary limits right now. + maxChartArea = NV_FLOAT_MAX; + maxBoundaryLength = NV_FLOAT_MAX; + + proxyFitMetricWeight = 1.0f; + roundnessMetricWeight = 0.1f; + straightnessMetricWeight = 0.25f; + normalSeamMetricWeight = 1.0f; + textureSeamMetricWeight = 0.1f; +} + + + +void MeshCharts::computeCharts(const SegmentationSettings & settings, const Array<uint> & unchartedMaterialArray) +{ + Chart * vertexMap = NULL; + + if (unchartedMaterialArray.count() != 0) { + vertexMap = new Chart(); + vertexMap->buildVertexMap(m_mesh, unchartedMaterialArray); + + if (vertexMap->faceCount() == 0) { + delete vertexMap; + vertexMap = NULL; + } + } + + + AtlasBuilder builder(m_mesh); + + if (vertexMap != NULL) { + // Mark faces that do not need to be charted. + builder.markUnchartedFaces(vertexMap->faceArray()); + + m_chartArray.append(vertexMap); + } + + if (builder.facesLeft != 0) { + + // Tweak these values: + const float maxThreshold = 2; + const uint growFaceCount = 32; + const uint maxIterations = 4; + + builder.settings = settings; + + //builder.settings.proxyFitMetricWeight *= 0.75; // relax proxy fit weight during initial seed placement. + //builder.settings.roundnessMetricWeight = 0; + //builder.settings.straightnessMetricWeight = 0; + + // This seems a reasonable estimate. + uint maxSeedCount = max(6U, builder.facesLeft); + + // Create initial charts greedely. + nvDebug("### Placing seeds\n"); + builder.placeSeeds(maxThreshold, maxSeedCount); + nvDebug("### Placed %d seeds (max = %d)\n", builder.chartCount(), maxSeedCount); + + builder.updateProxies(); + + builder.mergeCharts(); + + #if 1 + nvDebug("### Relocating seeds\n"); + builder.relocateSeeds(); + + nvDebug("### Reset charts\n"); + builder.resetCharts(); + + if (vertexMap != NULL) { + builder.markUnchartedFaces(vertexMap->faceArray()); + } + + builder.settings = settings; + + nvDebug("### Growing charts\n"); + + // Restart process growing charts in parallel. + uint iteration = 0; + while (true) + { + if (!builder.growCharts(maxThreshold, growFaceCount)) + { + nvDebug("### Can't grow anymore\n"); + + // If charts cannot grow more: fill holes, merge charts, relocate seeds and start new iteration. + + nvDebug("### Filling holes\n"); + builder.fillHoles(maxThreshold); + nvDebug("### Using %d charts now\n", builder.chartCount()); + + builder.updateProxies(); + + nvDebug("### Merging charts\n"); + builder.mergeCharts(); + nvDebug("### Using %d charts now\n", builder.chartCount()); + + nvDebug("### Reseeding\n"); + if (!builder.relocateSeeds()) + { + nvDebug("### Cannot relocate seeds anymore\n"); + + // Done! + break; + } + + if (iteration == maxIterations) + { + nvDebug("### Reached iteration limit\n"); + break; + } + iteration++; + + nvDebug("### Reset charts\n"); + builder.resetCharts(); + + if (vertexMap != NULL) { + builder.markUnchartedFaces(vertexMap->faceArray()); + } + + nvDebug("### Growing charts\n"); + } + }; + #endif + + // Make sure no holes are left! + nvDebugCheck(builder.facesLeft == 0); + + const uint chartCount = builder.chartArray.count(); + for (uint i = 0; i < chartCount; i++) + { + Chart * chart = new Chart(); + m_chartArray.append(chart); + + chart->build(m_mesh, builder.chartFaces(i)); + } + } + + + const uint chartCount = m_chartArray.count(); + + // Build face indices. + m_faceChart.resize(m_mesh->faceCount()); + m_faceIndex.resize(m_mesh->faceCount()); + + for (uint i = 0; i < chartCount; i++) + { + const Chart * chart = m_chartArray[i]; + + const uint faceCount = chart->faceCount(); + for (uint f = 0; f < faceCount; f++) + { + uint idx = chart->faceAt(f); + m_faceChart[idx] = i; + m_faceIndex[idx] = f; + } + } + + // Build an exclusive prefix sum of the chart vertex counts. + m_chartVertexCountPrefixSum.resize(chartCount); + + if (chartCount > 0) + { + m_chartVertexCountPrefixSum[0] = 0; + + for (uint i = 1; i < chartCount; i++) + { + const Chart * chart = m_chartArray[i-1]; + m_chartVertexCountPrefixSum[i] = m_chartVertexCountPrefixSum[i-1] + chart->vertexCount(); + } + + m_totalVertexCount = m_chartVertexCountPrefixSum[chartCount - 1] + m_chartArray[chartCount-1]->vertexCount(); + } + else + { + m_totalVertexCount = 0; + } +} + + +void MeshCharts::parameterizeCharts() +{ + ParameterizationQuality globalParameterizationQuality; + + // Parameterize the charts. + uint diskCount = 0; + const uint chartCount = m_chartArray.count(); + for (uint i = 0; i < chartCount; i++)\ + { + Chart * chart = m_chartArray[i]; + + bool isValid = false; + + if (chart->isVertexMapped()) { + continue; + } + + if (chart->isDisk()) + { + diskCount++; + + ParameterizationQuality chartParameterizationQuality; + + if (chart->faceCount() == 1) { + computeSingleFaceMap(chart->unifiedMesh()); + + chartParameterizationQuality = ParameterizationQuality(chart->unifiedMesh()); + } + else { + computeOrthogonalProjectionMap(chart->unifiedMesh()); + ParameterizationQuality orthogonalQuality(chart->unifiedMesh()); + + computeLeastSquaresConformalMap(chart->unifiedMesh()); + ParameterizationQuality lscmQuality(chart->unifiedMesh()); + + // If the orthogonal projection produces better results, just use that. + // @@ It may be dangerous to do this, because isValid() does not detect self-overlaps. + // @@ Another problem is that with very thin patches with nearly zero parametric area, the results of our metric are not accurate. + /*if (orthogonalQuality.isValid() && orthogonalQuality.rmsStretchMetric() < lscmQuality.rmsStretchMetric()) { + computeOrthogonalProjectionMap(chart->unifiedMesh()); + chartParameterizationQuality = orthogonalQuality; + } + else*/ { + chartParameterizationQuality = lscmQuality; + } + + // If conformal map failed, + + // @@ Experiment with other parameterization methods. + //computeCircularBoundaryMap(chart->unifiedMesh()); + //computeConformalMap(chart->unifiedMesh()); + //computeNaturalConformalMap(chart->unifiedMesh()); + //computeGuidanceGradientMap(chart->unifiedMesh()); + } + + //ParameterizationQuality chartParameterizationQuality(chart->unifiedMesh()); + + isValid = chartParameterizationQuality.isValid(); + + if (!isValid) + { + nvDebug("*** Invalid parameterization.\n"); +#if 0 + // Dump mesh to inspect problem: + static int pieceCount = 0; + + StringBuilder fileName; + fileName.format("invalid_chart_%d.obj", pieceCount++); + exportMesh(chart->unifiedMesh(), fileName.str()); +#endif + } + + // @@ Check that parameterization quality is above a certain threshold. + + // @@ Detect boundary self-intersections. + + globalParameterizationQuality += chartParameterizationQuality; + } + + if (!isValid) + { + //nvDebugBreak(); + // @@ Run the builder again, but only on this chart. + //AtlasBuilder builder(chart->chartMesh()); + } + + // Transfer parameterization from unified mesh to chart mesh. + chart->transferParameterization(); + + } + + nvDebug(" Parameterized %d/%d charts.\n", diskCount, chartCount); + nvDebug(" RMS stretch metric: %f\n", globalParameterizationQuality.rmsStretchMetric()); + nvDebug(" MAX stretch metric: %f\n", globalParameterizationQuality.maxStretchMetric()); + nvDebug(" RMS conformal metric: %f\n", globalParameterizationQuality.rmsConformalMetric()); + nvDebug(" RMS authalic metric: %f\n", globalParameterizationQuality.maxAuthalicMetric()); +} + + + +Chart::Chart() : m_chartMesh(NULL), m_unifiedMesh(NULL), m_isDisk(false), m_isVertexMapped(false) +{ +} + +void Chart::build(const HalfEdge::Mesh * originalMesh, const Array<uint> & faceArray) +{ + // Copy face indices. + m_faceArray = faceArray; + + const uint meshVertexCount = originalMesh->vertexCount(); + + m_chartMesh = new HalfEdge::Mesh(); + m_unifiedMesh = new HalfEdge::Mesh(); + + Array<uint> chartMeshIndices; + chartMeshIndices.resize(meshVertexCount, ~0); + + Array<uint> unifiedMeshIndices; + unifiedMeshIndices.resize(meshVertexCount, ~0); + + // Add vertices. + const uint faceCount = faceArray.count(); + for (uint f = 0; f < faceCount; f++) + { + const HalfEdge::Face * face = originalMesh->faceAt(faceArray[f]); + nvDebugCheck(face != NULL); + + for(HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + const HalfEdge::Vertex * vertex = it.current()->vertex; + const HalfEdge::Vertex * unifiedVertex = vertex->firstColocal(); + + if (unifiedMeshIndices[unifiedVertex->id] == ~0) + { + unifiedMeshIndices[unifiedVertex->id] = m_unifiedMesh->vertexCount(); + + nvDebugCheck(vertex->pos == unifiedVertex->pos); + m_unifiedMesh->addVertex(vertex->pos); + } + + if (chartMeshIndices[vertex->id] == ~0) + { + chartMeshIndices[vertex->id] = m_chartMesh->vertexCount(); + m_chartToOriginalMap.append(vertex->id); + m_chartToUnifiedMap.append(unifiedMeshIndices[unifiedVertex->id]); + + HalfEdge::Vertex * v = m_chartMesh->addVertex(vertex->pos); + v->nor = vertex->nor; + v->tex = vertex->tex; + } + } + } + + // This is ignoring the canonical map: + // - Is it really necessary to link colocals? + + m_chartMesh->linkColocals(); + //m_unifiedMesh->linkColocals(); // Not strictly necessary, no colocals in the unified mesh. # Wrong. + + // This check is not valid anymore, if the original mesh vertices were linked with a canonical map, then it might have + // some colocal vertices that were unlinked. So, the unified mesh might have some duplicate vertices, because firstColocal() + // is not guaranteed to return the same vertex for two colocal vertices. + //nvCheck(m_chartMesh->colocalVertexCount() == m_unifiedMesh->vertexCount()); + + // Is that OK? What happens in meshes were that happens? Does anything break? Apparently not... + + + + Array<uint> faceIndices(7); + + // Add faces. + for (uint f = 0; f < faceCount; f++) + { + const HalfEdge::Face * face = originalMesh->faceAt(faceArray[f]); + nvDebugCheck(face != NULL); + + faceIndices.clear(); + + for(HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + const HalfEdge::Vertex * vertex = it.current()->vertex; + nvDebugCheck(vertex != NULL); + + faceIndices.append(chartMeshIndices[vertex->id]); + } + + m_chartMesh->addFace(faceIndices); + + faceIndices.clear(); + + for(HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + const HalfEdge::Vertex * vertex = it.current()->vertex; + nvDebugCheck(vertex != NULL); + + vertex = vertex->firstColocal(); + + faceIndices.append(unifiedMeshIndices[vertex->id]); + } + + m_unifiedMesh->addFace(faceIndices); + } + + m_chartMesh->linkBoundary(); + m_unifiedMesh->linkBoundary(); + + //exportMesh(m_unifiedMesh.ptr(), "debug_input.obj"); + + if (m_unifiedMesh->splitBoundaryEdges()) { + m_unifiedMesh = unifyVertices(m_unifiedMesh.ptr()); + } + + //exportMesh(m_unifiedMesh.ptr(), "debug_split.obj"); + + // Closing the holes is not always the best solution and does not fix all the problems. + // We need to do some analysis of the holes and the genus to: + // - Find cuts that reduce genus. + // - Find cuts to connect holes. + // - Use minimal spanning trees or seamster. + if (!closeHoles()) { + /*static int pieceCount = 0; + StringBuilder fileName; + fileName.format("debug_hole_%d.obj", pieceCount++); + exportMesh(m_unifiedMesh.ptr(), fileName.str());*/ + } + + m_unifiedMesh = triangulate(m_unifiedMesh.ptr()); + + //exportMesh(m_unifiedMesh.ptr(), "debug_triangulated.obj"); + + + // Analyze chart topology. + MeshTopology topology(m_unifiedMesh.ptr()); + m_isDisk = topology.isDisk(); + + // This is sometimes failing, when triangulate fails to add a triangle, it generates a hole in the mesh. + //nvDebugCheck(m_isDisk); + + /*if (!m_isDisk) { + static int pieceCount = 0; + StringBuilder fileName; + fileName.format("debug_hole_%d.obj", pieceCount++); + exportMesh(m_unifiedMesh.ptr(), fileName.str()); + }*/ + + +#if 0 + if (!m_isDisk) { + nvDebugBreak(); + + static int pieceCount = 0; + + StringBuilder fileName; + fileName.format("debug_nodisk_%d.obj", pieceCount++); + exportMesh(m_chartMesh.ptr(), fileName.str()); + } +#endif + +} + + +void Chart::buildVertexMap(const HalfEdge::Mesh * originalMesh, const Array<uint> & unchartedMaterialArray) +{ + nvCheck(m_chartMesh == NULL && m_unifiedMesh == NULL); + + m_isVertexMapped = true; + + // Build face indices. + m_faceArray.clear(); + + const uint meshFaceCount = originalMesh->faceCount(); + for (uint f = 0; f < meshFaceCount; f++) { + const HalfEdge::Face * face = originalMesh->faceAt(f); + + if (unchartedMaterialArray.contains(face->material)) { + m_faceArray.append(f); + } + } + + const uint faceCount = m_faceArray.count(); + + if (faceCount == 0) { + return; + } + + + // @@ The chartMesh construction is basically the same as with regular charts, don't duplicate! + + const uint meshVertexCount = originalMesh->vertexCount(); + + m_chartMesh = new HalfEdge::Mesh(); + + Array<uint> chartMeshIndices; + chartMeshIndices.resize(meshVertexCount, ~0); + + // Vertex map mesh only has disconnected vertices. + for (uint f = 0; f < faceCount; f++) + { + const HalfEdge::Face * face = originalMesh->faceAt(m_faceArray[f]); + nvDebugCheck(face != NULL); + + for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + const HalfEdge::Vertex * vertex = it.current()->vertex; + + if (chartMeshIndices[vertex->id] == ~0) + { + chartMeshIndices[vertex->id] = m_chartMesh->vertexCount(); + m_chartToOriginalMap.append(vertex->id); + + HalfEdge::Vertex * v = m_chartMesh->addVertex(vertex->pos); + v->nor = vertex->nor; + v->tex = vertex->tex; // @@ Not necessary. + } + } + } + + // @@ Link colocals using the original mesh canonical map? Build canonical map on the fly? Do we need to link colocals at all for this? + //m_chartMesh->linkColocals(); + + Array<uint> faceIndices(7); + + // Add faces. + for (uint f = 0; f < faceCount; f++) + { + const HalfEdge::Face * face = originalMesh->faceAt(m_faceArray[f]); + nvDebugCheck(face != NULL); + + faceIndices.clear(); + + for(HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + const HalfEdge::Vertex * vertex = it.current()->vertex; + nvDebugCheck(vertex != NULL); + nvDebugCheck(chartMeshIndices[vertex->id] != ~0); + + faceIndices.append(chartMeshIndices[vertex->id]); + } + + HalfEdge::Face * new_face = m_chartMesh->addFace(faceIndices); + nvDebugCheck(new_face != NULL); + } + + m_chartMesh->linkBoundary(); + + + const uint chartVertexCount = m_chartMesh->vertexCount(); + + Box bounds; + bounds.clearBounds(); + + for (uint i = 0; i < chartVertexCount; i++) { + HalfEdge::Vertex * vertex = m_chartMesh->vertexAt(i); + bounds.addPointToBounds(vertex->pos); + } + + ProximityGrid grid; + grid.init(bounds, chartVertexCount); + + for (uint i = 0; i < chartVertexCount; i++) { + HalfEdge::Vertex * vertex = m_chartMesh->vertexAt(i); + grid.add(vertex->pos, i); + } + + +#if 0 + // Arrange vertices in a rectangle. + vertexMapWidth = ftoi_ceil(sqrtf(float(chartVertexCount))); + vertexMapHeight = (chartVertexCount + vertexMapWidth - 1) / vertexMapWidth; + nvDebugCheck(vertexMapWidth >= vertexMapHeight); + + int x = 0, y = 0; + for (uint i = 0; i < chartVertexCount; i++) { + HalfEdge::Vertex * vertex = m_chartMesh->vertexAt(i); + + vertex->tex.x = float(x); + vertex->tex.y = float(y); + + x++; + if (x == vertexMapWidth) { + x = 0; + y++; + nvCheck(y < vertexMapHeight); + } + } + +#elif 0 + // Arrange vertices in a rectangle, traversing grid in 3D morton order and laying them down in 2D morton order. + vertexMapWidth = ftoi_ceil(sqrtf(float(chartVertexCount))); + vertexMapHeight = (chartVertexCount + vertexMapWidth - 1) / vertexMapWidth; + nvDebugCheck(vertexMapWidth >= vertexMapHeight); + + int n = 0; + uint32 texelCode = 0; + + uint cellsVisited = 0; + + const uint32 cellCodeCount = grid.mortonCount(); + for (uint32 cellCode = 0; cellCode < cellCodeCount; cellCode++) { + int cell = grid.mortonIndex(cellCode); + if (cell < 0) continue; + + cellsVisited++; + + const Array<uint> & indexArray = grid.cellArray[cell].indexArray; + + foreach(i, indexArray) { + uint idx = indexArray[i]; + HalfEdge::Vertex * vertex = m_chartMesh->vertexAt(idx); + + //vertex->tex.x = float(n % rectangleWidth) + 0.5f; + //vertex->tex.y = float(n / rectangleWidth) + 0.5f; + + // Lay down the points in z order too. + uint x, y; + do { + x = decodeMorton2X(texelCode); + y = decodeMorton2Y(texelCode); + texelCode++; + } while (x >= U32(vertexMapWidth) || y >= U32(vertexMapHeight)); + + vertex->tex.x = float(x); + vertex->tex.y = float(y); + + n++; + } + } + + nvDebugCheck(cellsVisited == grid.cellArray.count()); + nvDebugCheck(n == chartVertexCount); + +#else + + uint texelCount = 0; + + const float positionThreshold = 0.01f; + const float normalThreshold = 0.01f; + + uint verticesVisited = 0; + uint cellsVisited = 0; + + Array<int> vertexIndexArray; + vertexIndexArray.resize(chartVertexCount, -1); // Init all indices to -1. + + // Traverse vertices in morton order. @@ It may be more interesting to sort them based on orientation. + const uint cellCodeCount = grid.mortonCount(); + for (uint cellCode = 0; cellCode < cellCodeCount; cellCode++) { + int cell = grid.mortonIndex(cellCode); + if (cell < 0) continue; + + cellsVisited++; + + const Array<uint> & indexArray = grid.cellArray[cell].indexArray; + + foreach(i, indexArray) { + uint idx = indexArray[i]; + HalfEdge::Vertex * vertex = m_chartMesh->vertexAt(idx); + + nvDebugCheck(vertexIndexArray[idx] == -1); + + Array<uint> neighbors; + grid.gather(vertex->pos, positionThreshold, /*ref*/neighbors); + + // Compare against all nearby vertices, cluster greedily. + foreach(j, neighbors) { + uint otherIdx = neighbors[j]; + + if (vertexIndexArray[otherIdx] != -1) { + HalfEdge::Vertex * otherVertex = m_chartMesh->vertexAt(otherIdx); + + if (distance(vertex->pos, otherVertex->pos) < positionThreshold && + distance(vertex->nor, otherVertex->nor) < normalThreshold) + { + vertexIndexArray[idx] = vertexIndexArray[otherIdx]; + break; + } + } + } + + // If index not assigned, assign new one. + if (vertexIndexArray[idx] == -1) { + vertexIndexArray[idx] = texelCount++; + } + + verticesVisited++; + } + } + + nvDebugCheck(cellsVisited == grid.cellArray.count()); + nvDebugCheck(verticesVisited == chartVertexCount); + + vertexMapWidth = ftoi_ceil(sqrtf(float(texelCount))); + vertexMapWidth = (vertexMapWidth + 3) & ~3; // Width aligned to 4. + vertexMapHeight = vertexMapWidth == 0 ? 0 : (texelCount + vertexMapWidth - 1) / vertexMapWidth; + //vertexMapHeight = (vertexMapHeight + 3) & ~3; // Height aligned to 4. + nvDebugCheck(vertexMapWidth >= vertexMapHeight); + + nvDebug("Reduced vertex count from %d to %d.\n", chartVertexCount, texelCount); + +#if 0 + // This lays down the clustered vertices linearly. + for (uint i = 0; i < chartVertexCount; i++) { + HalfEdge::Vertex * vertex = m_chartMesh->vertexAt(i); + + int idx = vertexIndexArray[i]; + + vertex->tex.x = float(idx % vertexMapWidth); + vertex->tex.y = float(idx / vertexMapWidth); + } +#else + // Lay down the clustered vertices in morton order. + + Array<uint> texelCodes; + texelCodes.resize(texelCount); + + // For each texel, assign one morton code. + uint texelCode = 0; + for (uint i = 0; i < texelCount; i++) { + uint x, y; + do { + x = decodeMorton2X(texelCode); + y = decodeMorton2Y(texelCode); + texelCode++; + } while (x >= U32(vertexMapWidth) || y >= U32(vertexMapHeight)); + + texelCodes[i] = texelCode - 1; + } + + for (uint i = 0; i < chartVertexCount; i++) { + HalfEdge::Vertex * vertex = m_chartMesh->vertexAt(i); + + int idx = vertexIndexArray[i]; + if (idx != -1) { + uint texelCode = texelCodes[idx]; + uint x = decodeMorton2X(texelCode); + uint y = decodeMorton2Y(texelCode); + + vertex->tex.x = float(x); + vertex->tex.y = float(y); + } + } + +#endif + +#endif + +} + + + +static void getBoundaryEdges(HalfEdge::Mesh * mesh, Array<HalfEdge::Edge *> & boundaryEdges) +{ + nvDebugCheck(mesh != NULL); + + const uint edgeCount = mesh->edgeCount(); + + BitArray bitFlags(edgeCount); + bitFlags.clearAll(); + + boundaryEdges.clear(); + + // Search for boundary edges. Mark all the edges that belong to the same boundary. + for (uint e = 0; e < edgeCount; e++) + { + HalfEdge::Edge * startEdge = mesh->edgeAt(e); + + if (startEdge != NULL && startEdge->isBoundary() && bitFlags.bitAt(e) == false) + { + nvDebugCheck(startEdge->face != NULL); + nvDebugCheck(startEdge->pair->face == NULL); + + startEdge = startEdge->pair; + + const HalfEdge::Edge * edge = startEdge; + do { + nvDebugCheck(edge->face == NULL); + nvDebugCheck(bitFlags.bitAt(edge->id/2) == false); + + bitFlags.setBitAt(edge->id / 2); + edge = edge->next; + } while(startEdge != edge); + + boundaryEdges.append(startEdge); + } + } +} + + +bool Chart::closeLoop(uint start, const Array<HalfEdge::Edge *> & loop) +{ + const uint vertexCount = loop.count() - start; + + nvDebugCheck(vertexCount >= 3); + if (vertexCount < 3) return false; + + nvDebugCheck(loop[start]->vertex->isColocal(loop[start+vertexCount-1]->to())); + + // If the hole is planar, then we add a single face that will be properly triangulated later. + // If the hole is not planar, we add a triangle fan with a vertex at the hole centroid. + // This is still a bit of a hack. There surely are better hole filling algorithms out there. + + Array<Vector3> points; + points.resize(vertexCount); + for (uint i = 0; i < vertexCount; i++) { + points[i] = loop[start+i]->vertex->pos; + } + + bool isPlanar = Fit::isPlanar(vertexCount, points.buffer()); + + if (isPlanar) { + // Add face and connect edges. + HalfEdge::Face * face = m_unifiedMesh->addFace(); + for (uint i = 0; i < vertexCount; i++) { + HalfEdge::Edge * edge = loop[start + i]; + + edge->face = face; + edge->setNext(loop[start + (i + 1) % vertexCount]); + } + face->edge = loop[start]; + + nvDebugCheck(face->isValid()); + } + else { + // If the polygon is not planar, we just cross our fingers, and hope this will work: + + // Compute boundary centroid: + Vector3 centroidPos(0); + + for (uint i = 0; i < vertexCount; i++) { + centroidPos += points[i]; + } + + centroidPos *= (1.0f / vertexCount); + + HalfEdge::Vertex * centroid = m_unifiedMesh->addVertex(centroidPos); + + // Add one pair of edges for each boundary vertex. + for (uint j = vertexCount-1, i = 0; i < vertexCount; j = i++) { + HalfEdge::Face * face = m_unifiedMesh->addFace(centroid->id, loop[start+j]->vertex->id, loop[start+i]->vertex->id); + nvDebugCheck(face != NULL); + } + } + + return true; +} + + +bool Chart::closeHoles() +{ + nvDebugCheck(!m_isVertexMapped); + + Array<HalfEdge::Edge *> boundaryEdges; + getBoundaryEdges(m_unifiedMesh.ptr(), boundaryEdges); + + uint boundaryCount = boundaryEdges.count(); + if (boundaryCount <= 1) + { + // Nothing to close. + return true; + } + + // Compute lengths and areas. + Array<float> boundaryLengths; + //Array<Vector3> boundaryCentroids; + + for (uint i = 0; i < boundaryCount; i++) + { + const HalfEdge::Edge * startEdge = boundaryEdges[i]; + nvCheck(startEdge->face == NULL); + + //float boundaryEdgeCount = 0; + float boundaryLength = 0.0f; + //Vector3 boundaryCentroid(zero); + + const HalfEdge::Edge * edge = startEdge; + do { + Vector3 t0 = edge->from()->pos; + Vector3 t1 = edge->to()->pos; + + //boundaryEdgeCount++; + boundaryLength += length(t1 - t0); + //boundaryCentroid += edge->vertex()->pos; + + edge = edge->next; + } while(edge != startEdge); + + boundaryLengths.append(boundaryLength); + //boundaryCentroids.append(boundaryCentroid / boundaryEdgeCount); + } + + + // Find disk boundary. + uint diskBoundary = 0; + float maxLength = boundaryLengths[0]; + + for (uint i = 1; i < boundaryCount; i++) + { + if (boundaryLengths[i] > maxLength) + { + maxLength = boundaryLengths[i]; + diskBoundary = i; + } + } + + + // Sew holes. + /*for (uint i = 0; i < boundaryCount; i++) + { + if (diskBoundary == i) + { + // Skip disk boundary. + continue; + } + + HalfEdge::Edge * startEdge = boundaryEdges[i]; + nvCheck(startEdge->face() == NULL); + + boundaryEdges[i] = m_unifiedMesh->sewBoundary(startEdge); + } + + exportMesh(m_unifiedMesh.ptr(), "debug_sewn.obj");*/ + + //bool hasNewHoles = false; + + // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + // @@ Close loop is wrong, after closing a loop, we do not only have to add the face, but make sure that every edge in he loop is pointing to the right place. + + // Close holes. + for (uint i = 0; i < boundaryCount; i++) + { + if (diskBoundary == i) + { + // Skip disk boundary. + continue; + } + + HalfEdge::Edge * startEdge = boundaryEdges[i]; + nvDebugCheck(startEdge != NULL); + nvDebugCheck(startEdge->face == NULL); + +#if 1 + Array<HalfEdge::Vertex *> vertexLoop; + Array<HalfEdge::Edge *> edgeLoop; + + HalfEdge::Edge * edge = startEdge; + do { + HalfEdge::Vertex * vertex = edge->next->vertex; // edge->to() + + uint i; + for (i = 0; i < vertexLoop.count(); i++) { + if (vertex->isColocal(vertexLoop[i])) { + break; + } + } + + bool isCrossing = (i != vertexLoop.count()); + + if (isCrossing) { + + HalfEdge::Edge * prev = edgeLoop[i]; // Previous edge before the loop. + HalfEdge::Edge * next = edge->next; // Next edge after the loop. + + nvDebugCheck(prev->to()->isColocal(next->from())); + + // Close loop. + edgeLoop.append(edge); + closeLoop(i+1, edgeLoop); + + // Link boundary loop. + prev->setNext(next); + vertex->setEdge(next); + + // Start over again. + vertexLoop.clear(); + edgeLoop.clear(); + + edge = startEdge; + vertex = edge->to(); + } + + vertexLoop.append(vertex); + edgeLoop.append(edge); + + edge = edge->next; + } while(edge != startEdge); + + closeLoop(0, edgeLoop); +#endif + + /* + + // Add face and connect boundary edges. + HalfEdge::Face * face = m_unifiedMesh->addFace(); + face->setEdge(startEdge); + + HalfEdge::Edge * edge = startEdge; + do { + edge->setFace(face); + + edge = edge->next(); + } while(edge != startEdge); + + */ + + + /* + uint edgeCount = 0; + HalfEdge::Edge * edge = startEdge; + do { + edgeCount++; + edge = edge->next(); + } while(edge != startEdge); + + + + // Count edges in this boundary. + uint edgeCount = 0; + HalfEdge::Edge * edge = startEdge; + do { + edgeCount++; + edge = edge->next(); + } while(edge != startEdge); + + // Trivial hole, fill with one triangle. This actually works for all convex boundaries with non colinear vertices. + if (edgeCount == 3) { + // Add face and connect boundary edges. + HalfEdge::Face * face = m_unifiedMesh->addFace(); + face->setEdge(startEdge); + + edge = startEdge; + do { + edge->setFace(face); + + edge = edge->next(); + } while(edge != startEdge); + + // @@ Implement the above using addFace, it should now work with existing edges, as long as their face pointers is zero. + + } + else { + // Ideally we should: + // - compute best fit plane of boundary vertices. + // - project boundary polygon onto plane. + // - triangulate boundary polygon. + // - add faces of the resulting triangulation. + + // I don't have a good triangulator available. A more simple solution that works in more (but not all) cases: + // - compute boundary centroid. + // - add vertex centroid. + // - connect centroid vertex with boundary vertices. + // - connect radial edges with boundary edges. + + // This should work for non-convex boundaries with colinear vertices as long as the kernel of the polygon is not empty. + + // Compute boundary centroid: + Vector3 centroid_pos(0); + Vector2 centroid_tex(0); + + HalfEdge::Edge * edge = startEdge; + do { + centroid_pos += edge->vertex()->pos; + centroid_tex += edge->vertex()->tex; + edge = edge->next(); + } while(edge != startEdge); + + centroid_pos *= (1.0f / edgeCount); + centroid_tex *= (1.0f / edgeCount); + + HalfEdge::Vertex * centroid = m_unifiedMesh->addVertex(centroid_pos); + centroid->tex = centroid_tex; + + // Add one pair of edges for each boundary vertex. + edge = startEdge; + do { + HalfEdge::Edge * next = edge->next(); + + nvCheck(edge->face() == NULL); + HalfEdge::Face * face = m_unifiedMesh->addFace(centroid->id(), edge->from()->id(), edge->to()->id()); + + if (face != NULL) { + nvCheck(edge->face() == face); + } + else { + hasNewHoles = true; + } + + edge = next; + } while(edge != startEdge); + } + */ + } + + /*nvDebugCheck(!hasNewHoles); + + if (hasNewHoles) { + // Link boundary again, in case closeHoles created new holes! + m_unifiedMesh->linkBoundary(); + }*/ + + // Because some algorithms do not expect sparse edge buffers. + //m_unifiedMesh->compactEdges(); + + // In case we messed up: + //m_unifiedMesh->linkBoundary(); + + getBoundaryEdges(m_unifiedMesh.ptr(), boundaryEdges); + + boundaryCount = boundaryEdges.count(); + nvDebugCheck(boundaryCount == 1); + + //exportMesh(m_unifiedMesh.ptr(), "debug_hole_filled.obj"); + + return boundaryCount == 1; +} + + +// Transfer parameterization from unified mesh to chart mesh. +void Chart::transferParameterization() { + nvDebugCheck(!m_isVertexMapped); + + uint vertexCount = m_chartMesh->vertexCount(); + for (uint v = 0; v < vertexCount; v++) { + HalfEdge::Vertex * vertex = m_chartMesh->vertexAt(v); + HalfEdge::Vertex * unifiedVertex = m_unifiedMesh->vertexAt(mapChartVertexToUnifiedVertex(v)); + vertex->tex = unifiedVertex->tex; + } +} + +float Chart::computeSurfaceArea() const { + return nv::computeSurfaceArea(m_chartMesh.ptr()) * scale; +} + +float Chart::computeParametricArea() const { + // This only makes sense in parameterized meshes. + nvDebugCheck(m_isDisk); + nvDebugCheck(!m_isVertexMapped); + + return nv::computeParametricArea(m_chartMesh.ptr()); +} + +Vector2 Chart::computeParametricBounds() const { + // This only makes sense in parameterized meshes. + nvDebugCheck(m_isDisk); + nvDebugCheck(!m_isVertexMapped); + + Box bounds; + bounds.clearBounds(); + + uint vertexCount = m_chartMesh->vertexCount(); + for (uint v = 0; v < vertexCount; v++) { + HalfEdge::Vertex * vertex = m_chartMesh->vertexAt(v); + bounds.addPointToBounds(Vector3(vertex->tex, 0)); + } + + return bounds.extents().xy(); +} diff --git a/thirdparty/thekla_atlas/nvmesh/param/Atlas.h b/thirdparty/thekla_atlas/nvmesh/param/Atlas.h new file mode 100644 index 0000000000..41cfaea9cb --- /dev/null +++ b/thirdparty/thekla_atlas/nvmesh/param/Atlas.h @@ -0,0 +1,186 @@ +// Copyright NVIDIA Corporation 2006 -- Ignacio Castano <icastano@nvidia.com> + +#pragma once +#ifndef NV_MESH_ATLAS_H +#define NV_MESH_ATLAS_H + +#include "nvcore/Array.h" +#include "nvcore/Ptr.h" +#include "nvmath/Vector.h" +#include "nvmesh/nvmesh.h" +#include "nvmesh/halfedge/Mesh.h" + + +namespace nv +{ + namespace HalfEdge { class Mesh; } + + class Chart; + class MeshCharts; + class VertexMap; + + struct SegmentationSettings + { + SegmentationSettings(); + + float maxChartArea; + float maxBoundaryLength; + + float proxyFitMetricWeight; + float roundnessMetricWeight; + float straightnessMetricWeight; + float normalSeamMetricWeight; + float textureSeamMetricWeight; + }; + + + /// An atlas is a set of charts. + class Atlas + { + public: + + Atlas(); + ~Atlas(); + + uint meshCount() const { return m_meshChartsArray.count(); } + const MeshCharts * meshAt(uint i) const { return m_meshChartsArray[i]; } + MeshCharts * meshAt(uint i) { return m_meshChartsArray[i]; } + + uint chartCount() const; + const Chart * chartAt(uint i) const; + Chart * chartAt(uint i); + + // Add mesh charts and takes ownership. + void addMeshCharts(MeshCharts * meshCharts); + + void extractCharts(const HalfEdge::Mesh * mesh); + void computeCharts(const HalfEdge::Mesh * mesh, const SegmentationSettings & settings, const Array<uint> & unchartedMaterialArray); + + + // Compute a trivial seamless texture similar to ZBrush. + //bool computeSeamlessTextureAtlas(bool groupFaces = true, bool scaleTiles = false, uint w = 1024, uint h = 1024); + + void parameterizeCharts(); + + // Pack charts in the smallest possible rectangle. + float packCharts(int quality, float texelArea, bool blockAlign, bool conservative); + void setFailed() { failed = true; } + bool hasFailed() const { return failed; } + + private: + + bool failed; + Array<MeshCharts *> m_meshChartsArray; + + }; + + + // Set of charts corresponding to a single mesh. + class MeshCharts + { + public: + MeshCharts(const HalfEdge::Mesh * mesh); + ~MeshCharts(); + + uint chartCount() const { return m_chartArray.count(); } + uint vertexCount () const { return m_totalVertexCount; } + + const Chart * chartAt(uint i) const { return m_chartArray[i]; } + Chart * chartAt(uint i) { return m_chartArray[i]; } + + void computeVertexMap(const Array<uint> & unchartedMaterialArray); + + // Extract the charts of the input mesh. + void extractCharts(); + + // Compute charts using a simple segmentation algorithm. + void computeCharts(const SegmentationSettings & settings, const Array<uint> & unchartedMaterialArray); + + void parameterizeCharts(); + + uint faceChartAt(uint i) const { return m_faceChart[i]; } + uint faceIndexWithinChartAt(uint i) const { return m_faceIndex[i]; } + + uint vertexCountBeforeChartAt(uint i) const { return m_chartVertexCountPrefixSum[i]; } + + private: + + const HalfEdge::Mesh * m_mesh; + + Array<Chart *> m_chartArray; + + Array<uint> m_chartVertexCountPrefixSum; + uint m_totalVertexCount; + + Array<uint> m_faceChart; // the chart of every face of the input mesh. + Array<uint> m_faceIndex; // the index within the chart for every face of the input mesh. + }; + + + /// A chart is a connected set of faces with a certain topology (usually a disk). + class Chart + { + public: + + Chart(); + + void build(const HalfEdge::Mesh * originalMesh, const Array<uint> & faceArray); + void buildVertexMap(const HalfEdge::Mesh * originalMesh, const Array<uint> & unchartedMaterialArray); + + bool closeHoles(); + + bool isDisk() const { return m_isDisk; } + bool isVertexMapped() const { return m_isVertexMapped; } + + uint vertexCount() const { return m_chartMesh->vertexCount(); } + uint colocalVertexCount() const { return m_unifiedMesh->vertexCount(); } + + uint faceCount() const { return m_faceArray.count(); } + uint faceAt(uint i) const { return m_faceArray[i]; } + + const HalfEdge::Mesh * chartMesh() const { return m_chartMesh.ptr(); } + HalfEdge::Mesh * chartMesh() { return m_chartMesh.ptr(); } + const HalfEdge::Mesh * unifiedMesh() const { return m_unifiedMesh.ptr(); } + HalfEdge::Mesh * unifiedMesh() { return m_unifiedMesh.ptr(); } + + //uint vertexIndex(uint i) const { return m_vertexIndexArray[i]; } + + uint mapChartVertexToOriginalVertex(uint i) const { return m_chartToOriginalMap[i]; } + uint mapChartVertexToUnifiedVertex(uint i) const { return m_chartToUnifiedMap[i]; } + + const Array<uint> & faceArray() const { return m_faceArray; } + + void transferParameterization(); + + float computeSurfaceArea() const; + float computeParametricArea() const; + Vector2 computeParametricBounds() const; + + + float scale = 1.0f; + uint vertexMapWidth; + uint vertexMapHeight; + + private: + + bool closeLoop(uint start, const Array<HalfEdge::Edge *> & loop); + + // Chart mesh. + AutoPtr<HalfEdge::Mesh> m_chartMesh; + AutoPtr<HalfEdge::Mesh> m_unifiedMesh; + + bool m_isDisk; + bool m_isVertexMapped; + + // List of faces of the original mesh that belong to this chart. + Array<uint> m_faceArray; + + // Map vertices of the chart mesh to vertices of the original mesh. + Array<uint> m_chartToOriginalMap; + + Array<uint> m_chartToUnifiedMap; + }; + +} // nv namespace + +#endif // NV_MESH_ATLAS_H diff --git a/thirdparty/thekla_atlas/nvmesh/param/AtlasBuilder.cpp b/thirdparty/thekla_atlas/nvmesh/param/AtlasBuilder.cpp new file mode 100644 index 0000000000..bd2140c2f3 --- /dev/null +++ b/thirdparty/thekla_atlas/nvmesh/param/AtlasBuilder.cpp @@ -0,0 +1,1320 @@ +// This code is in the public domain -- castano@gmail.com + +#include "nvmesh.h" // pch + +#include "AtlasBuilder.h" +#include "Util.h" + +#include "nvmesh/halfedge/Mesh.h" +#include "nvmesh/halfedge/Face.h" +#include "nvmesh/halfedge/Vertex.h" + +#include "nvmath/Matrix.inl" +#include "nvmath/Vector.inl" + +//#include "nvcore/IntroSort.h" +#include "nvcore/Array.inl" + +#include <algorithm> // std::sort + +#include <float.h> // FLT_MAX +#include <limits.h> // UINT_MAX + +using namespace nv; + +namespace +{ + + // Dummy implementation of a priority queue using sort at insertion. + // - Insertion is o(n) + // - Smallest element goes at the end, so that popping it is o(1). + // - Resorting is n*log(n) + // @@ Number of elements in the queue is usually small, and we'd have to rebalance often. I'm not sure it's worth implementing a heap. + // @@ Searcing at removal would remove the need for sorting when priorities change. + struct PriorityQueue + { + PriorityQueue(uint size = UINT_MAX) : maxSize(size) {} + + void push(float priority, uint face) { + uint i = 0; + const uint count = pairs.count(); + for (; i < count; i++) { + if (pairs[i].priority > priority) break; + } + + Pair p = { priority, face }; + pairs.insertAt(i, p); + + if (pairs.count() > maxSize) { + pairs.removeAt(0); + } + } + + // push face out of order, to be sorted later. + void push(uint face) { + Pair p = { 0.0f, face }; + pairs.append(p); + } + + uint pop() { + uint f = pairs.back().face; + pairs.pop_back(); + return f; + } + + void sort() { + //nv::sort(pairs); // @@ My intro sort appears to be much slower than it should! + std::sort(pairs.buffer(), pairs.buffer() + pairs.count()); + } + + void clear() { + pairs.clear(); + } + + uint count() const { return pairs.count(); } + + float firstPriority() const { return pairs.back().priority; } + + + const uint maxSize; + + struct Pair { + bool operator <(const Pair & p) const { return priority > p.priority; } // !! Sort in inverse priority order! + float priority; + uint face; + }; + + + Array<Pair> pairs; + }; + + static bool isNormalSeam(const HalfEdge::Edge * edge) { + return (edge->vertex->nor != edge->pair->next->vertex->nor || edge->next->vertex->nor != edge->pair->vertex->nor); + } + + static bool isTextureSeam(const HalfEdge::Edge * edge) { + return (edge->vertex->tex != edge->pair->next->vertex->tex || edge->next->vertex->tex != edge->pair->vertex->tex); + } + +} // namespace + + +struct nv::ChartBuildData +{ + ChartBuildData(int id) : id(id) { + planeNormal = Vector3(0); + centroid = Vector3(0); + coneAxis = Vector3(0); + coneAngle = 0; + area = 0; + boundaryLength = 0; + normalSum = Vector3(0); + centroidSum = Vector3(0); + } + + int id; + + // Proxy info: + Vector3 planeNormal; + Vector3 centroid; + Vector3 coneAxis; + float coneAngle; + + float area; + float boundaryLength; + Vector3 normalSum; + Vector3 centroidSum; + + Array<uint> seeds; // @@ These could be a pointers to the HalfEdge faces directly. + Array<uint> faces; + PriorityQueue candidates; +}; + + + +AtlasBuilder::AtlasBuilder(const HalfEdge::Mesh * m) : mesh(m), facesLeft(m->faceCount()) +{ + const uint faceCount = m->faceCount(); + faceChartArray.resize(faceCount, -1); + faceCandidateArray.resize(faceCount, -1); + + // @@ Floyd for the whole mesh is too slow. We could compute floyd progressively per patch as the patch grows. We need a better solution to compute most central faces. + //computeShortestPaths(); + + // Precompute edge lengths and face areas. + uint edgeCount = m->edgeCount(); + edgeLengths.resize(edgeCount); + + for (uint i = 0; i < edgeCount; i++) { + uint id = m->edgeAt(i)->id; + nvDebugCheck(id / 2 == i); + + edgeLengths[i] = m->edgeAt(i)->length(); + } + + faceAreas.resize(faceCount); + for (uint i = 0; i < faceCount; i++) { + faceAreas[i] = m->faceAt(i)->area(); + } +} + +AtlasBuilder::~AtlasBuilder() +{ + const uint chartCount = chartArray.count(); + for (uint i = 0; i < chartCount; i++) + { + delete chartArray[i]; + } +} + + +void AtlasBuilder::markUnchartedFaces(const Array<uint> & unchartedFaces) +{ + const uint unchartedFaceCount = unchartedFaces.count(); + for (uint i = 0; i < unchartedFaceCount; i++){ + uint f = unchartedFaces[i]; + faceChartArray[f] = -2; + //faceCandidateArray[f] = -2; // @@ ? + + removeCandidate(f); + } + + nvDebugCheck(facesLeft >= unchartedFaceCount); + facesLeft -= unchartedFaceCount; +} + + +void AtlasBuilder::computeShortestPaths() +{ + const uint faceCount = mesh->faceCount(); + shortestPaths.resize(faceCount*faceCount, FLT_MAX); + + // Fill edges: + for (uint i = 0; i < faceCount; i++) + { + shortestPaths[i*faceCount + i] = 0.0f; + + const HalfEdge::Face * face_i = mesh->faceAt(i); + Vector3 centroid_i = face_i->centroid(); + + for (HalfEdge::Face::ConstEdgeIterator it(face_i->edges()); !it.isDone(); it.advance()) + { + const HalfEdge::Edge * edge = it.current(); + + if (!edge->isBoundary()) + { + const HalfEdge::Face * face_j = edge->pair->face; + + uint j = face_j->id; + Vector3 centroid_j = face_j->centroid(); + + shortestPaths[i*faceCount + j] = shortestPaths[j*faceCount + i] = length(centroid_i - centroid_j); + } + } + } + + // Use Floyd-Warshall algorithm to compute all paths: + for (uint k = 0; k < faceCount; k++) + { + for (uint i = 0; i < faceCount; i++) + { + for (uint j = 0; j < faceCount; j++) + { + shortestPaths[i*faceCount + j] = min(shortestPaths[i*faceCount + j], shortestPaths[i*faceCount + k]+shortestPaths[k*faceCount + j]); + } + } + } +} + + +void AtlasBuilder::placeSeeds(float threshold, uint maxSeedCount) +{ + // Instead of using a predefiened number of seeds: + // - Add seeds one by one, growing chart until a certain treshold. + // - Undo charts and restart growing process. + + // @@ How can we give preference to faces far from sharp features as in the LSCM paper? + // - those points can be found using a simple flood filling algorithm. + // - how do we weight the probabilities? + + for (uint i = 0; i < maxSeedCount; i++) + { + if (facesLeft == 0) { + // No faces left, stop creating seeds. + break; + } + + createRandomChart(threshold); + } +} + + +void AtlasBuilder::createRandomChart(float threshold) +{ + ChartBuildData * chart = new ChartBuildData(chartArray.count()); + chartArray.append(chart); + + // Pick random face that is not used by any chart yet. + uint randomFaceIdx = rand.getRange(facesLeft - 1); + uint i = 0; + for (uint f = 0; f != randomFaceIdx; f++, i++) + { + while (faceChartArray[i] != -1) i++; + } + while (faceChartArray[i] != -1) i++; + + chart->seeds.append(i); + + addFaceToChart(chart, i, true); + + // Grow the chart as much as possible within the given threshold. + growChart(chart, threshold * 0.5f, facesLeft); + //growCharts(threshold - threshold * 0.75f / chartCount(), facesLeft); +} + +void AtlasBuilder::addFaceToChart(ChartBuildData * chart, uint f, bool recomputeProxy) +{ + // Add face to chart. + chart->faces.append(f); + + nvDebugCheck(faceChartArray[f] == -1); + faceChartArray[f] = chart->id; + + facesLeft--; + + // Update area and boundary length. + chart->area = evaluateChartArea(chart, f); + chart->boundaryLength = evaluateBoundaryLength(chart, f); + chart->normalSum = evaluateChartNormalSum(chart, f); + chart->centroidSum = evaluateChartCentroidSum(chart, f); + + if (recomputeProxy) { + // Update proxy and candidate's priorities. + updateProxy(chart); + } + + // Update candidates. + removeCandidate(f); + updateCandidates(chart, f); + updatePriorities(chart); +} + +// @@ Get N best candidates in one pass. +const AtlasBuilder::Candidate & AtlasBuilder::getBestCandidate() const +{ + uint best = 0; + float bestCandidateMetric = FLT_MAX; + + const uint candidateCount = candidateArray.count(); + nvCheck(candidateCount > 0); + + for (uint i = 0; i < candidateCount; i++) + { + const Candidate & candidate = candidateArray[i]; + + if (candidate.metric < bestCandidateMetric) { + bestCandidateMetric = candidate.metric; + best = i; + } + } + + return candidateArray[best]; +} + + +// Returns true if any of the charts can grow more. +bool AtlasBuilder::growCharts(float threshold, uint faceCount) +{ +#if 1 // Using one global list. + + faceCount = min(faceCount, facesLeft); + + for (uint i = 0; i < faceCount; i++) + { + const Candidate & candidate = getBestCandidate(); + + if (candidate.metric > threshold) { + return false; // Can't grow more. + } + + addFaceToChart(candidate.chart, candidate.face); + } + + return facesLeft != 0; // Can continue growing. + +#else // Using one list per chart. + bool canGrowMore = false; + + const uint chartCount = chartArray.count(); + for (uint i = 0; i < chartCount; i++) + { + if (growChart(chartArray[i], threshold, faceCount)) + { + canGrowMore = true; + } + } + + return canGrowMore; +#endif +} + +bool AtlasBuilder::growChart(ChartBuildData * chart, float threshold, uint faceCount) +{ + // Try to add faceCount faces within threshold to chart. + for (uint i = 0; i < faceCount; ) + { + if (chart->candidates.count() == 0 || chart->candidates.firstPriority() > threshold) + { + return false; + } + + uint f = chart->candidates.pop(); + if (faceChartArray[f] == -1) + { + addFaceToChart(chart, f); + i++; + } + } + + if (chart->candidates.count() == 0 || chart->candidates.firstPriority() > threshold) + { + return false; + } + + return true; +} + + +void AtlasBuilder::resetCharts() +{ + const uint faceCount = mesh->faceCount(); + for (uint i = 0; i < faceCount; i++) + { + faceChartArray[i] = -1; + faceCandidateArray[i] = -1; + } + + facesLeft = faceCount; + + candidateArray.clear(); + + const uint chartCount = chartArray.count(); + for (uint i = 0; i < chartCount; i++) + { + ChartBuildData * chart = chartArray[i]; + + const uint seed = chart->seeds.back(); + + chart->area = 0.0f; + chart->boundaryLength = 0.0f; + chart->normalSum = Vector3(0); + chart->centroidSum = Vector3(0); + + chart->faces.clear(); + chart->candidates.clear(); + + addFaceToChart(chart, seed); + } +} + + +void AtlasBuilder::updateCandidates(ChartBuildData * chart, uint f) +{ + const HalfEdge::Face * face = mesh->faceAt(f); + + // Traverse neighboring faces, add the ones that do not belong to any chart yet. + for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + const HalfEdge::Edge * edge = it.current()->pair; + + if (!edge->isBoundary()) + { + uint f = edge->face->id; + + if (faceChartArray[f] == -1) + { + chart->candidates.push(f); + } + } + } +} + + +void AtlasBuilder::updateProxies() +{ + const uint chartCount = chartArray.count(); + for (uint i = 0; i < chartCount; i++) + { + updateProxy(chartArray[i]); + } +} + + +namespace { + + float absoluteSum(Vector4::Arg v) + { + return fabs(v.x) + fabs(v.y) + fabs(v.z) + fabs(v.w); + } + + //#pragma message(NV_FILE_LINE "FIXME: Using the c=cos(teta) substitution, the equation system becomes linear and we can avoid the newton solver.") + + struct ConeFitting + { + ConeFitting(const HalfEdge::Mesh * m, float g, float tf, float tx) : mesh(m), gamma(g), tolf(tf), tolx(tx), F(0), D(0), H(0) { + } + + void addTerm(Vector3 N, float A) + { + const float c = cosf(X.w); + const float s = sinf(X.w); + const float tmp = dot(X.xyz(), N) - c; + + F += tmp * tmp; + + D.x += 2 * X.x * tmp; + D.y += 2 * X.y * tmp; + D.z += 2 * X.z * tmp; + D.w += 2 * s * tmp; + + H(0,0) = 2 * X.x * N.x + 2 * tmp; + H(0,1) = 2 * X.x * N.y; + H(0,2) = 2 * X.x * N.z; + H(0,3) = 2 * X.x * s; + + H(1,0) = 2 * X.y * N.x; + H(1,1) = 2 * X.y * N.y + 2 * tmp; + H(1,2) = 2 * X.y * N.z; + H(1,3) = 2 * X.y * s; + + H(2,0) = 2 * X.z * N.x; + H(2,1) = 2 * X.z * N.y; + H(2,2) = 2 * X.z * N.z + 2 * tmp; + H(2,3) = 2 * X.z * s; + + H(3,0) = 2 * s * N.x; + H(3,1) = 2 * s * N.y; + H(3,2) = 2 * s * N.z; + H(3,3) = 2 * s * s + 2 * c * tmp; + } + + Vector4 solve(ChartBuildData * chart, Vector4 start) + { + const uint faceCount = chart->faces.count(); + + X = start; + + Vector4 dX; + + do { + for (uint i = 0; i < faceCount; i++) + { + const HalfEdge::Face * face = mesh->faceAt(chart->faces[i]); + + addTerm(face->normal(), face->area()); + } + + Vector4 dX; + //solveKramer(H, D, &dX); + solveLU(H, D, &dX); + + // @@ Do a full newton step and reduce by half if F doesn't decrease. + X -= gamma * dX; + + // Constrain normal to be normalized. + X = Vector4(normalize(X.xyz()), X.w); + + } while(absoluteSum(D) > tolf || absoluteSum(dX) > tolx); + + return X; + } + + HalfEdge::Mesh const * const mesh; + const float gamma; + const float tolf; + const float tolx; + + Vector4 X; + + float F; + Vector4 D; + Matrix H; + }; + + // Unnormalized face normal assuming it's a triangle. + static Vector3 triangleNormal(const HalfEdge::Face * face) + { + Vector3 p0 = face->edge->vertex->pos; + Vector3 p1 = face->edge->next->vertex->pos; + Vector3 p2 = face->edge->next->next->vertex->pos; + + Vector3 e0 = p2 - p0; + Vector3 e1 = p1 - p0; + + return normalizeSafe(cross(e0, e1), Vector3(0), 0.0f); + } + + static Vector3 triangleNormalAreaScaled(const HalfEdge::Face * face) + { + Vector3 p0 = face->edge->vertex->pos; + Vector3 p1 = face->edge->next->vertex->pos; + Vector3 p2 = face->edge->next->next->vertex->pos; + + Vector3 e0 = p2 - p0; + Vector3 e1 = p1 - p0; + + return cross(e0, e1); + } + + // Average of the edge midpoints weighted by the edge length. + // I want a point inside the triangle, but closer to the cirumcenter. + static Vector3 triangleCenter(const HalfEdge::Face * face) + { + Vector3 p0 = face->edge->vertex->pos; + Vector3 p1 = face->edge->next->vertex->pos; + Vector3 p2 = face->edge->next->next->vertex->pos; + + float l0 = length(p1 - p0); + float l1 = length(p2 - p1); + float l2 = length(p0 - p2); + + Vector3 m0 = (p0 + p1) * l0 / (l0 + l1 + l2); + Vector3 m1 = (p1 + p2) * l1 / (l0 + l1 + l2); + Vector3 m2 = (p2 + p0) * l2 / (l0 + l1 + l2); + + return m0 + m1 + m2; + } + +} // namespace + +void AtlasBuilder::updateProxy(ChartBuildData * chart) +{ + //#pragma message(NV_FILE_LINE "TODO: Use best fit plane instead of average normal.") + + chart->planeNormal = normalizeSafe(chart->normalSum, Vector3(0), 0.0f); + chart->centroid = chart->centroidSum / float(chart->faces.count()); + + //#pragma message(NV_FILE_LINE "TODO: Experiment with conic fitting.") + + // F = (Nc*Nt - cos Oc)^2 = (x*Nt_x + y*Nt_y + z*Nt_z - cos w)^2 + // dF/dx = 2 * x * (x*Nt_x + y*Nt_y + z*Nt_z - cos w) + // dF/dy = 2 * y * (x*Nt_x + y*Nt_y + z*Nt_z - cos w) + // dF/dz = 2 * z * (x*Nt_x + y*Nt_y + z*Nt_z - cos w) + // dF/dw = 2 * sin w * (x*Nt_x + y*Nt_y + z*Nt_z - cos w) + + // JacobianMatrix({ + // 2 * x * (x*Nt_x + y*Nt_y + z*Nt_z - Cos(w)), + // 2 * y * (x*Nt_x + y*Nt_y + z*Nt_z - Cos(w)), + // 2 * z * (x*Nt_x + y*Nt_y + z*Nt_z - Cos(w)), + // 2 * Sin(w) * (x*Nt_x + y*Nt_y + z*Nt_z - Cos(w))}, {x,y,z,w}) + + // H[0,0] = 2 * x * Nt_x + 2 * (x*Nt_x + y*Nt_y + z*Nt_z - cos(w)); + // H[0,1] = 2 * x * Nt_y; + // H[0,2] = 2 * x * Nt_z; + // H[0,3] = 2 * x * sin(w); + + // H[1,0] = 2 * y * Nt_x; + // H[1,1] = 2 * y * Nt_y + 2 * (x*Nt_x + y*Nt_y + z*Nt_z - cos(w)); + // H[1,2] = 2 * y * Nt_z; + // H[1,3] = 2 * y * sin(w); + + // H[2,0] = 2 * z * Nt_x; + // H[2,1] = 2 * z * Nt_y; + // H[2,2] = 2 * z * Nt_z + 2 * (x*Nt_x + y*Nt_y + z*Nt_z - cos(w)); + // H[2,3] = 2 * z * sin(w); + + // H[3,0] = 2 * sin(w) * Nt_x; + // H[3,1] = 2 * sin(w) * Nt_y; + // H[3,2] = 2 * sin(w) * Nt_z; + // H[3,3] = 2 * sin(w) * sin(w) + 2 * cos(w) * (x*Nt_x + y*Nt_y + z*Nt_z - cos(w)); + + // @@ Cone fitting might be quite slow. + + /*ConeFitting coneFitting(mesh, 0.1f, 0.001f, 0.001f); + + Vector4 start = Vector4(chart->coneAxis, chart->coneAngle); + Vector4 solution = coneFitting.solve(chart, start); + + chart->coneAxis = solution.xyz(); + chart->coneAngle = solution.w;*/ +} + + + +bool AtlasBuilder::relocateSeeds() +{ + bool anySeedChanged = false; + + const uint chartCount = chartArray.count(); + for (uint i = 0; i < chartCount; i++) + { + if (relocateSeed(chartArray[i])) + { + anySeedChanged = true; + } + } + + return anySeedChanged; +} + + +bool AtlasBuilder::relocateSeed(ChartBuildData * chart) +{ + Vector3 centroid = computeChartCentroid(chart); + + const uint N = 10; // @@ Hardcoded to 10? + PriorityQueue bestTriangles(N); + + // Find the first N triangles that fit the proxy best. + const uint faceCount = chart->faces.count(); + for (uint i = 0; i < faceCount; i++) + { + float priority = evaluateProxyFitMetric(chart, chart->faces[i]); + bestTriangles.push(priority, chart->faces[i]); + } + + // Of those, choose the most central triangle. + uint mostCentral; + float maxDistance = -1; + + const uint bestCount = bestTriangles.count(); + for (uint i = 0; i < bestCount; i++) + { + const HalfEdge::Face * face = mesh->faceAt(bestTriangles.pairs[i].face); + Vector3 faceCentroid = triangleCenter(face); + + float distance = length(centroid - faceCentroid); + + /*#pragma message(NV_FILE_LINE "TODO: Implement evaluateDistanceToBoundary.") + float distance = evaluateDistanceToBoundary(chart, bestTriangles.pairs[i].face);*/ + + if (distance > maxDistance) + { + maxDistance = distance; + mostCentral = bestTriangles.pairs[i].face; + } + } + nvDebugCheck(maxDistance >= 0); + + // In order to prevent k-means cyles we record all the previously chosen seeds. + uint index; + if (chart->seeds.find(mostCentral, &index)) + { + // Move new seed to the end of the seed array. + uint last = chart->seeds.count() - 1; + swap(chart->seeds[index], chart->seeds[last]); + return false; + } + else + { + // Append new seed. + chart->seeds.append(mostCentral); + return true; + } +} + +void AtlasBuilder::removeCandidate(uint f) +{ + int c = faceCandidateArray[f]; + if (c != -1) { + faceCandidateArray[f] = -1; + + if (c == candidateArray.count() - 1) { + candidateArray.popBack(); + } + else { + candidateArray.replaceWithLast(c); + faceCandidateArray[candidateArray[c].face] = c; + } + } +} + +void AtlasBuilder::updateCandidate(ChartBuildData * chart, uint f, float metric) +{ + if (faceCandidateArray[f] == -1) { + const uint index = candidateArray.count(); + faceCandidateArray[f] = index; + candidateArray.resize(index + 1); + candidateArray[index].face = f; + candidateArray[index].chart = chart; + candidateArray[index].metric = metric; + } + else { + int c = faceCandidateArray[f]; + nvDebugCheck(c != -1); + + Candidate & candidate = candidateArray[c]; + nvDebugCheck(candidate.face == f); + + if (metric < candidate.metric || chart == candidate.chart) { + candidate.metric = metric; + candidate.chart = chart; + } + } + +} + + +void AtlasBuilder::updatePriorities(ChartBuildData * chart) +{ + // Re-evaluate candidate priorities. + uint candidateCount = chart->candidates.count(); + for (uint i = 0; i < candidateCount; i++) + { + chart->candidates.pairs[i].priority = evaluatePriority(chart, chart->candidates.pairs[i].face); + + if (faceChartArray[chart->candidates.pairs[i].face] == -1) + { + updateCandidate(chart, chart->candidates.pairs[i].face, chart->candidates.pairs[i].priority); + } + } + + // Sort candidates. + chart->candidates.sort(); +} + + +// Evaluate combined metric. +float AtlasBuilder::evaluatePriority(ChartBuildData * chart, uint face) +{ + // Estimate boundary length and area: + float newBoundaryLength = evaluateBoundaryLength(chart, face); + float newChartArea = evaluateChartArea(chart, face); + + float F = evaluateProxyFitMetric(chart, face); + float C = evaluateRoundnessMetric(chart, face, newBoundaryLength, newChartArea); + float P = evaluateStraightnessMetric(chart, face); + + // Penalize faces that cross seams, reward faces that close seams or reach boundaries. + float N = evaluateNormalSeamMetric(chart, face); + float T = evaluateTextureSeamMetric(chart, face); + + //float R = evaluateCompletenessMetric(chart, face); + + //float D = evaluateDihedralAngleMetric(chart, face); + // @@ Add a metric based on local dihedral angle. + + // @@ Tweaking the normal and texture seam metrics. + // - Cause more impedance. Never cross 90 degree edges. + // - + + float cost = float( + settings.proxyFitMetricWeight * F + + settings.roundnessMetricWeight * C + + settings.straightnessMetricWeight * P + + settings.normalSeamMetricWeight * N + + settings.textureSeamMetricWeight * T); + + /*cost = settings.proxyFitMetricWeight * powf(F, settings.proxyFitMetricExponent); + cost = max(cost, settings.roundnessMetricWeight * powf(C, settings.roundnessMetricExponent)); + cost = max(cost, settings.straightnessMetricWeight * pow(P, settings.straightnessMetricExponent)); + cost = max(cost, settings.normalSeamMetricWeight * N); + cost = max(cost, settings.textureSeamMetricWeight * T);*/ + + // Enforce limits strictly: + if (newChartArea > settings.maxChartArea) cost = FLT_MAX; + if (newBoundaryLength > settings.maxBoundaryLength) cost = FLT_MAX; + + // Make sure normal seams are fully respected: + if (settings.normalSeamMetricWeight >= 1000 && N != 0) cost = FLT_MAX; + + nvCheck(isFinite(cost)); + return cost; +} + + +// Returns a value in [0-1]. +float AtlasBuilder::evaluateProxyFitMetric(ChartBuildData * chart, uint f) +{ + const HalfEdge::Face * face = mesh->faceAt(f); + Vector3 faceNormal = triangleNormal(face); + //return square(dot(chart->coneAxis, faceNormal) - cosf(chart->coneAngle)); + + // Use plane fitting metric for now: + //return square(1 - dot(faceNormal, chart->planeNormal)); // @@ normal deviations should be weighted by face area + return 1 - dot(faceNormal, chart->planeNormal); // @@ normal deviations should be weighted by face area + + // Find distance to chart. + /*Vector3 faceCentroid = face->centroid(); + + float dist = 0; + int count = 0; + + for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + const HalfEdge::Edge * edge = it.current(); + + if (!edge->isBoundary()) { + const HalfEdge::Face * neighborFace = edge->pair()->face(); + if (faceChartArray[neighborFace->id()] == chart->id) { + dist += length(neighborFace->centroid() - faceCentroid); + count++; + } + } + } + + dist /= (count * count); + + return (1 - dot(faceNormal, chart->planeNormal)) * dist;*/ + + //return (1 - dot(faceNormal, chart->planeNormal)); +} + +float AtlasBuilder::evaluateDistanceToBoundary(ChartBuildData * chart, uint face) +{ +//#pragma message(NV_FILE_LINE "TODO: Evaluate distance to boundary metric.") + + // @@ This is needed for the seed relocation code. + // @@ This could provide a better roundness metric. + + return 0.0f; +} + +float AtlasBuilder::evaluateDistanceToSeed(ChartBuildData * chart, uint f) +{ + //const uint seed = chart->seeds.back(); + //const uint faceCount = mesh->faceCount(); + //return shortestPaths[seed * faceCount + f]; + + const HalfEdge::Face * seed = mesh->faceAt(chart->seeds.back()); + const HalfEdge::Face * face = mesh->faceAt(f); + return length(triangleCenter(seed) - triangleCenter(face)); +} + + +float AtlasBuilder::evaluateRoundnessMetric(ChartBuildData * chart, uint face, float newBoundaryLength, float newChartArea) +{ + // @@ D-charts use distance to seed. + // C(c,t) = pi * D(S_c,t)^2 / A_c + //return PI * square(evaluateDistanceToSeed(chart, face)) / chart->area; + //return PI * square(evaluateDistanceToSeed(chart, face)) / chart->area; + //return 2 * PI * evaluateDistanceToSeed(chart, face) / chart->boundaryLength; + + // Garland's Hierarchical Face Clustering paper uses ratio between boundary and area, which is easier to compute and might work as well: + // roundness = D^2/4*pi*A -> circle = 1, non circle greater than 1 + + //return square(newBoundaryLength) / (newChartArea * 4 * PI); + float roundness = square(chart->boundaryLength) / chart->area; + float newRoundness = square(newBoundaryLength) / newChartArea; + if (newRoundness > roundness) { + return square(newBoundaryLength) / (newChartArea * 4 * PI); + } + else { + // Offer no impedance to faces that improve roundness. + return 0; + } + + //return square(newBoundaryLength) / (4 * PI * newChartArea); + //return clamp(1 - (4 * PI * newChartArea) / square(newBoundaryLength), 0.0f, 1.0f); + + // Use the ratio between the new roundness vs. the previous roundness. + // - If we use the absolute metric, when the initial face is very long, then it's hard to make any progress. + //return (square(newBoundaryLength) * chart->area) / (square(chart->boundaryLength) * newChartArea); + //return (4 * PI * newChartArea) / square(newBoundaryLength) - (4 * PI * chart->area) / square(chart->boundaryLength); + + //if (square(newBoundaryLength) * chart->area) / (square(chart->boundaryLength) * newChartArea); + +} + +float AtlasBuilder::evaluateStraightnessMetric(ChartBuildData * chart, uint f) +{ + float l_out = 0.0f; + float l_in = 0.0f; + + const HalfEdge::Face * face = mesh->faceAt(f); + for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + const HalfEdge::Edge * edge = it.current(); + + //float l = edge->length(); + float l = edgeLengths[edge->id/2]; + + if (edge->isBoundary()) + { + l_out += l; + } + else + { + uint neighborFaceId = edge->pair->face->id; + if (faceChartArray[neighborFaceId] != chart->id) { + l_out += l; + } + else { + l_in += l; + } + } + } + nvDebugCheck(l_in != 0.0f); // Candidate face must be adjacent to chart. @@ This is not true if the input mesh has zero-length edges. + + //return l_out / l_in; + float ratio = (l_out - l_in) / (l_out + l_in); + //if (ratio < 0) ratio *= 10; // Encourage closing gaps. + return min(ratio, 0.0f); // Only use the straightness metric to close gaps. + //return ratio; +} + + +float AtlasBuilder::evaluateNormalSeamMetric(ChartBuildData * chart, uint f) +{ + float seamFactor = 0.0f; + float totalLength = 0.0f; + + const HalfEdge::Face * face = mesh->faceAt(f); + for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + const HalfEdge::Edge * edge = it.current(); + + if (edge->isBoundary()) { + continue; + } + + const uint neighborFaceId = edge->pair->face->id; + if (faceChartArray[neighborFaceId] != chart->id) { + continue; + } + + //float l = edge->length(); + float l = edgeLengths[edge->id/2]; + + totalLength += l; + + if (!edge->isSeam()) { + continue; + } + + // Make sure it's a normal seam. + if (isNormalSeam(edge)) + { + float d0 = clamp(dot(edge->vertex->nor, edge->pair->next->vertex->nor), 0.0f, 1.0f); + float d1 = clamp(dot(edge->next->vertex->nor, edge->pair->vertex->nor), 0.0f, 1.0f); + //float a0 = clamp(acosf(d0) / (PI/2), 0.0f, 1.0f); + //float a1 = clamp(acosf(d1) / (PI/2), 0.0f, 1.0f); + //l *= (a0 + a1) * 0.5f; + + l *= 1 - (d0 + d1) * 0.5f; + + seamFactor += l; + } + } + + if (seamFactor == 0) return 0.0f; + return seamFactor / totalLength; +} + + +float AtlasBuilder::evaluateTextureSeamMetric(ChartBuildData * chart, uint f) +{ + float seamLength = 0.0f; + //float newSeamLength = 0.0f; + //float oldSeamLength = 0.0f; + float totalLength = 0.0f; + + const HalfEdge::Face * face = mesh->faceAt(f); + for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + const HalfEdge::Edge * edge = it.current(); + + /*float l = edge->length(); + totalLength += l; + + if (edge->isBoundary() || !edge->isSeam()) { + continue; + } + + // Make sure it's a texture seam. + if (isTextureSeam(edge)) + { + uint neighborFaceId = edge->pair()->face()->id(); + if (faceChartArray[neighborFaceId] != chart->id) { + newSeamLength += l; + } + else { + oldSeamLength += l; + } + }*/ + + if (edge->isBoundary()) { + continue; + } + + const uint neighborFaceId = edge->pair->face->id; + if (faceChartArray[neighborFaceId] != chart->id) { + continue; + } + + //float l = edge->length(); + float l = edgeLengths[edge->id/2]; + totalLength += l; + + if (!edge->isSeam()) { + continue; + } + + // Make sure it's a texture seam. + if (isTextureSeam(edge)) + { + seamLength += l; + } + } + + if (seamLength == 0.0f) { + return 0.0f; // Avoid division by zero. + } + + return seamLength / totalLength; +} + + +float AtlasBuilder::evaluateSeamMetric(ChartBuildData * chart, uint f) +{ + float newSeamLength = 0.0f; + float oldSeamLength = 0.0f; + float totalLength = 0.0f; + + const HalfEdge::Face * face = mesh->faceAt(f); + for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + const HalfEdge::Edge * edge = it.current(); + + //float l = edge->length(); + float l = edgeLengths[edge->id/2]; + + if (edge->isBoundary()) + { + newSeamLength += l; + } + else + { + if (edge->isSeam()) + { + uint neighborFaceId = edge->pair->face->id; + if (faceChartArray[neighborFaceId] != chart->id) { + newSeamLength += l; + } + else { + oldSeamLength += l; + } + } + } + + totalLength += l; + } + + return (newSeamLength - oldSeamLength) / totalLength; +} + + +float AtlasBuilder::evaluateChartArea(ChartBuildData * chart, uint f) +{ + const HalfEdge::Face * face = mesh->faceAt(f); + //return chart->area + face->area(); + return chart->area + faceAreas[face->id]; +} + + +float AtlasBuilder::evaluateBoundaryLength(ChartBuildData * chart, uint f) +{ + float boundaryLength = chart->boundaryLength; + + // Add new edges, subtract edges shared with the chart. + const HalfEdge::Face * face = mesh->faceAt(f); + for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + const HalfEdge::Edge * edge = it.current(); + //float edgeLength = edge->length(); + float edgeLength = edgeLengths[edge->id/2]; + + if (edge->isBoundary()) + { + boundaryLength += edgeLength; + } + else + { + uint neighborFaceId = edge->pair->face->id; + if (faceChartArray[neighborFaceId] != chart->id) { + boundaryLength += edgeLength; + } + else { + boundaryLength -= edgeLength; + } + } + } + //nvDebugCheck(boundaryLength >= 0); + + return max(0.0f, boundaryLength); // @@ Hack! +} + +Vector3 AtlasBuilder::evaluateChartNormalSum(ChartBuildData * chart, uint f) +{ + const HalfEdge::Face * face = mesh->faceAt(f); + return chart->normalSum + triangleNormalAreaScaled(face); +} + +Vector3 AtlasBuilder::evaluateChartCentroidSum(ChartBuildData * chart, uint f) +{ + const HalfEdge::Face * face = mesh->faceAt(f); + return chart->centroidSum + face->centroid(); +} + + +Vector3 AtlasBuilder::computeChartCentroid(const ChartBuildData * chart) +{ + Vector3 centroid(0); + + const uint faceCount = chart->faces.count(); + for (uint i = 0; i < faceCount; i++) + { + const HalfEdge::Face * face = mesh->faceAt(chart->faces[i]); + centroid += triangleCenter(face); + } + + return centroid / float(faceCount); +} + + +void AtlasBuilder::fillHoles(float threshold) +{ + while (facesLeft > 0) + { + createRandomChart(threshold); + } +} + + +void AtlasBuilder::mergeChart(ChartBuildData * owner, ChartBuildData * chart, float sharedBoundaryLength) +{ + const uint faceCount = chart->faces.count(); + for (uint i = 0; i < faceCount; i++) + { + uint f = chart->faces[i]; + + nvDebugCheck(faceChartArray[f] == chart->id); + faceChartArray[f] = owner->id; + + owner->faces.append(f); + } + + // Update adjacencies? + + owner->area += chart->area; + owner->boundaryLength += chart->boundaryLength - sharedBoundaryLength; + + owner->normalSum += chart->normalSum; + owner->centroidSum += chart->centroidSum; + + updateProxy(owner); +} + +void AtlasBuilder::mergeCharts() +{ + Array<float> sharedBoundaryLengths; + + const uint chartCount = chartArray.count(); + for (int c = chartCount-1; c >= 0; c--) + { + sharedBoundaryLengths.clear(); + sharedBoundaryLengths.resize(chartCount, 0.0f); + + ChartBuildData * chart = chartArray[c]; + + float externalBoundary = 0.0f; + + const uint faceCount = chart->faces.count(); + for (uint i = 0; i < faceCount; i++) + { + uint f = chart->faces[i]; + const HalfEdge::Face * face = mesh->faceAt(f); + + for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + const HalfEdge::Edge * edge = it.current(); + + //float l = edge->length(); + float l = edgeLengths[edge->id/2]; + + if (edge->isBoundary()) { + externalBoundary += l; + } + else { + uint neighborFace = edge->pair->face->id; + uint neighborChart = faceChartArray[neighborFace]; + + if (neighborChart != c) { + if ((edge->isSeam() && (isNormalSeam(edge) || isTextureSeam(edge))) || neighborChart == -2) { + externalBoundary += l; + } + else { + sharedBoundaryLengths[neighborChart] += l; + } + } + } + } + } + + for (int cc = chartCount-1; cc >= 0; cc--) + { + if (cc == c) + continue; + + ChartBuildData * chart2 = chartArray[cc]; + if (chart2 == NULL) + continue; + + if (sharedBoundaryLengths[cc] > 0.8 * max(0.0f, chart->boundaryLength - externalBoundary)) { + + // Try to avoid degenerate configurations. + if (chart2->boundaryLength > sharedBoundaryLengths[cc]) + { + if (dot(chart2->planeNormal, chart->planeNormal) > -0.25) { + mergeChart(chart2, chart, sharedBoundaryLengths[cc]); + delete chart; + chartArray[c] = NULL; + break; + } + } + } + + if (sharedBoundaryLengths[cc] > 0.20 * max(0.0f, chart->boundaryLength - externalBoundary)) { + + // Compare proxies. + if (dot(chart2->planeNormal, chart->planeNormal) > 0) { + mergeChart(chart2, chart, sharedBoundaryLengths[cc]); + delete chart; + chartArray[c] = NULL; + break; + } + } + } + } + + // Remove deleted charts. + for (int c = 0; c < I32(chartArray.count()); /*do not increment if removed*/) + { + if (chartArray[c] == NULL) { + chartArray.removeAt(c); + + // Update faceChartArray. + const uint faceCount = faceChartArray.count(); + for (uint i = 0; i < faceCount; i++) { + nvDebugCheck (faceChartArray[i] != -1); + nvDebugCheck (faceChartArray[i] != c); + nvDebugCheck (faceChartArray[i] <= I32(chartArray.count())); + + if (faceChartArray[i] > c) { + faceChartArray[i]--; + } + } + } + else { + chartArray[c]->id = c; + c++; + } + } +} + + + +const Array<uint> & AtlasBuilder::chartFaces(uint i) const +{ + return chartArray[i]->faces; +} diff --git a/thirdparty/thekla_atlas/nvmesh/param/AtlasBuilder.h b/thirdparty/thekla_atlas/nvmesh/param/AtlasBuilder.h new file mode 100644 index 0000000000..f25c724f7e --- /dev/null +++ b/thirdparty/thekla_atlas/nvmesh/param/AtlasBuilder.h @@ -0,0 +1,111 @@ +// This code is in the public domain -- castano@gmail.com + +#pragma once +#ifndef NV_MESH_ATLASBUILDER_H +#define NV_MESH_ATLASBUILDER_H + +#include "Atlas.h" + +#include "nvmath/Vector.h" +#include "nvmath/Random.h" +#include "nvmesh/nvmesh.h" + +#include "nvcore/Array.h" +#include "nvcore/BitArray.h" + + + +namespace nv +{ + namespace HalfEdge { class Mesh; } + + struct ChartBuildData; + + struct AtlasBuilder + { + AtlasBuilder(const HalfEdge::Mesh * m); + ~AtlasBuilder(); + + void markUnchartedFaces(const Array<uint> & unchartedFaces); + + void computeShortestPaths(); + + void placeSeeds(float threshold, uint maxSeedCount); + void createRandomChart(float threshold); + + void addFaceToChart(ChartBuildData * chart, uint f, bool recomputeProxy=false); + + bool growCharts(float threshold, uint faceCount); + bool growChart(ChartBuildData * chart, float threshold, uint faceCount); + + void resetCharts(); + + void updateCandidates(ChartBuildData * chart, uint face); + + void updateProxies(); + void updateProxy(ChartBuildData * chart); + + bool relocateSeeds(); + bool relocateSeed(ChartBuildData * chart); + + void updatePriorities(ChartBuildData * chart); + + float evaluatePriority(ChartBuildData * chart, uint face); + float evaluateProxyFitMetric(ChartBuildData * chart, uint face); + float evaluateDistanceToBoundary(ChartBuildData * chart, uint face); + float evaluateDistanceToSeed(ChartBuildData * chart, uint face); + float evaluateRoundnessMetric(ChartBuildData * chart, uint face, float newBoundaryLength, float newChartArea); + float evaluateStraightnessMetric(ChartBuildData * chart, uint face); + + float evaluateNormalSeamMetric(ChartBuildData * chart, uint f); + float evaluateTextureSeamMetric(ChartBuildData * chart, uint f); + float evaluateSeamMetric(ChartBuildData * chart, uint f); + + float evaluateChartArea(ChartBuildData * chart, uint f); + float evaluateBoundaryLength(ChartBuildData * chart, uint f); + Vector3 evaluateChartNormalSum(ChartBuildData * chart, uint f); + Vector3 evaluateChartCentroidSum(ChartBuildData * chart, uint f); + + Vector3 computeChartCentroid(const ChartBuildData * chart); + + + void fillHoles(float threshold); + void mergeCharts(); + + // @@ Cleanup. + struct Candidate { + uint face; + ChartBuildData * chart; + float metric; + }; + + const Candidate & getBestCandidate() const; + void removeCandidate(uint f); + void updateCandidate(ChartBuildData * chart, uint f, float metric); + + void mergeChart(ChartBuildData * owner, ChartBuildData * chart, float sharedBoundaryLength); + + + uint chartCount() const { return chartArray.count(); } + const Array<uint> & chartFaces(uint i) const; + + const HalfEdge::Mesh * mesh; + uint facesLeft; + Array<int> faceChartArray; + Array<ChartBuildData *> chartArray; + Array<float> shortestPaths; + + Array<float> edgeLengths; + Array<float> faceAreas; + + Array<Candidate> candidateArray; // + Array<uint> faceCandidateArray; // Map face index to candidate index. + + MTRand rand; + + SegmentationSettings settings; + }; + +} // nv namespace + +#endif // NV_MESH_ATLASBUILDER_H diff --git a/thirdparty/thekla_atlas/nvmesh/param/AtlasPacker.cpp b/thirdparty/thekla_atlas/nvmesh/param/AtlasPacker.cpp new file mode 100644 index 0000000000..5ce452cb9e --- /dev/null +++ b/thirdparty/thekla_atlas/nvmesh/param/AtlasPacker.cpp @@ -0,0 +1,1387 @@ +// This code is in the public domain -- castano@gmail.com + +#include "nvmesh.h" // pch + +#include "AtlasPacker.h" +#include "nvmesh/halfedge/Vertex.h" +#include "nvmesh/halfedge/Face.h" +#include "nvmesh/param/Atlas.h" +#include "nvmesh/param/Util.h" +#include "nvmesh/raster/Raster.h" + +#include "nvmath/Vector.inl" +#include "nvmath/ConvexHull.h" +#include "nvmath/Color.h" +#include "nvmath/ftoi.h" + +#include "nvcore/StrLib.h" // debug +#include "nvcore/StdStream.h" // fileOpen + +#include <float.h> // FLT_MAX +#include <limits.h> // UINT_MAX + +using namespace nv; + +#define DEBUG_OUTPUT 0 + +#if DEBUG_OUTPUT + +#include "nvimage/ImageIO.h" + +namespace +{ + const uint TGA_TYPE_GREY = 3; + const uint TGA_TYPE_RGB = 2; + const uint TGA_ORIGIN_UPPER = 0x20; + +#pragma pack(push, 1) + struct TgaHeader { + uint8 id_length; + uint8 colormap_type; + uint8 image_type; + uint16 colormap_index; + uint16 colormap_length; + uint8 colormap_size; + uint16 x_origin; + uint16 y_origin; + uint16 width; + uint16 height; + uint8 pixel_size; + uint8 flags; + + enum { Size = 18 }; //const static int SIZE = 18; + }; +#pragma pack(pop) + + static void outputDebugBitmap(const char * fileName, const BitMap & bitmap, int w, int h) + { + FILE * fp = fileOpen(fileName, "wb"); + if (fp == NULL) return; + + nvStaticCheck(sizeof(TgaHeader) == TgaHeader::Size); + TgaHeader tga; + tga.id_length = 0; + tga.colormap_type = 0; + tga.image_type = TGA_TYPE_GREY; + + tga.colormap_index = 0; + tga.colormap_length = 0; + tga.colormap_size = 0; + + tga.x_origin = 0; + tga.y_origin = 0; + tga.width = w; + tga.height = h; + tga.pixel_size = 8; + tga.flags = TGA_ORIGIN_UPPER; + + fwrite(&tga, sizeof(TgaHeader), 1, fp); + + for (int j = 0; j < h; j++) { + for (int i = 0; i < w; i++) { + uint8 color = bitmap.bitAt(i, j) ? 0xFF : 0x0; + fwrite(&color, 1, 1, fp); + } + } + + fclose(fp); + } + + static void outputDebugImage(const char * fileName, const Image & bitmap, int w, int h) + { + FILE * fp = fileOpen(fileName, "wb"); + if (fp == NULL) return; + + nvStaticCheck(sizeof(TgaHeader) == TgaHeader::Size); + TgaHeader tga; + tga.id_length = 0; + tga.colormap_type = 0; + tga.image_type = TGA_TYPE_RGB; + + tga.colormap_index = 0; + tga.colormap_length = 0; + tga.colormap_size = 0; + + tga.x_origin = 0; + tga.y_origin = 0; + tga.width = w; + tga.height = h; + tga.pixel_size = 24; + tga.flags = TGA_ORIGIN_UPPER; + + fwrite(&tga, sizeof(TgaHeader), 1, fp); + + for (int j = 0; j < h; j++) { + for (int i = 0; i < w; i++) { + Color32 color = bitmap.pixel(i, j); + fwrite(&color.r, 1, 1, fp); + fwrite(&color.g, 1, 1, fp); + fwrite(&color.b, 1, 1, fp); + } + } + + fclose(fp); + } +} + +#endif // DEBUG_OUTPUT + +inline int align(int x, int a) { + //return a * ((x + a - 1) / a); + //return (x + a - 1) & -a; + return (x + a - 1) & ~(a - 1); +} + +inline bool isAligned(int x, int a) { + return (x & (a - 1)) == 0; +} + + + +AtlasPacker::AtlasPacker(Atlas * atlas) : m_atlas(atlas), m_bitmap(256, 256) +{ + m_width = 0; + m_height = 0; + + m_debug_bitmap.allocate(256, 256); + m_debug_bitmap.fill(Color32(0,0,0,0)); +} + +AtlasPacker::~AtlasPacker() +{ +} + +// This should compute convex hull and use rotating calipers to find the best box. Currently it uses a brute force method. +static bool computeBoundingBox(Chart * chart, Vector2 * majorAxis, Vector2 * minorAxis, Vector2 * minCorner, Vector2 * maxCorner) +{ + // Compute list of boundary points. + Array<Vector2> points(16); + + HalfEdge::Mesh * mesh = chart->chartMesh(); + const uint vertexCount = mesh->vertexCount(); + + for (uint i = 0; i < vertexCount; i++) { + HalfEdge::Vertex * vertex = mesh->vertexAt(i); + if (vertex->isBoundary()) { + points.append(vertex->tex); + } + } + + // This is not valid anymore. The chart mesh may have multiple boundaries! + /*const HalfEdge::Vertex * vertex = findBoundaryVertex(chart->chartMesh()); + + // Traverse boundary. + const HalfEdge::Edge * const firstEdge = vertex->edge(); + const HalfEdge::Edge * edge = firstEdge; + do { + vertex = edge->vertex(); + + nvDebugCheck (vertex->isBoundary()); + points.append(vertex->tex); + + edge = edge->next(); + } while (edge != firstEdge);*/ + +#if 1 + Array<Vector2> hull; + if (points.size()==0) { + return false; + } + + convexHull(points, hull, 0.00001f); + + // @@ Ideally I should use rotating calipers to find the best box. Using brute force for now. + + float best_area = FLT_MAX; + Vector2 best_min; + Vector2 best_max; + Vector2 best_axis; + + const uint hullCount = hull.count(); + for (uint i = 0, j = hullCount-1; i < hullCount; j = i, i++) { + + if (equal(hull[i], hull[j])) { + continue; + } + + Vector2 axis = normalize(hull[i] - hull[j], 0.0f); + nvDebugCheck(isFinite(axis)); + + // Compute bounding box. + Vector2 box_min(FLT_MAX, FLT_MAX); + Vector2 box_max(-FLT_MAX, -FLT_MAX); + + for (uint v = 0; v < hullCount; v++) { + + Vector2 point = hull[v]; + + float x = dot(axis, point); + if (x < box_min.x) box_min.x = x; + if (x > box_max.x) box_max.x = x; + + float y = dot(Vector2(-axis.y, axis.x), point); + if (y < box_min.y) box_min.y = y; + if (y > box_max.y) box_max.y = y; + } + + // Compute box area. + float area = (box_max.x - box_min.x) * (box_max.y - box_min.y); + + if (area < best_area) { + best_area = area; + best_min = box_min; + best_max = box_max; + best_axis = axis; + } + } + + // Make sure the box contains all the input points since the convex hull is not 100% accurate. + /*const uint pointCount = points.count(); + for (uint v = 0; v < pointCount; v++) { + + Vector2 point = points[v]; + + float x = dot(best_axis, point); + if (x < best_min.x) best_min.x = x; + + float y = dot(Vector2(-best_axis.y, best_axis.x), point); + if (y < best_min.y) best_min.y = y; + }*/ + + // Consider all points, not only boundary points, in case the input chart is malformed. + for (uint i = 0; i < vertexCount; i++) { + HalfEdge::Vertex * vertex = mesh->vertexAt(i); + Vector2 point = vertex->tex; + + float x = dot(best_axis, point); + if (x < best_min.x) best_min.x = x; + if (x > best_max.x) best_max.x = x; + + float y = dot(Vector2(-best_axis.y, best_axis.x), point); + if (y < best_min.y) best_min.y = y; + if (y > best_max.y) best_max.y = y; + } + + *majorAxis = best_axis; + *minorAxis = Vector2(-best_axis.y, best_axis.x); + *minCorner = best_min; + *maxCorner = best_max; + +#else + // Approximate implementation: try 16 different directions and keep the best. + + const uint N = 16; + Vector2 axis[N]; + + float minAngle = 0; + float maxAngle = PI / 2; + + int best; + Vector2 mins[N]; + Vector2 maxs[N]; + + const int iterationCount = 1; + for (int j = 0; j < iterationCount; j++) + { + // Init predefined directions. + for (int i = 0; i < N; i++) + { + float angle = lerp(minAngle, maxAngle, float(i)/N); + axis[i].set(cosf(angle), sinf(angle)); + } + + // Compute box for each direction. + for (int i = 0; i < N; i++) + { + mins[i].set(FLT_MAX, FLT_MAX); + maxs[i].set(-FLT_MAX, -FLT_MAX); + } + + for (uint p = 0; p < points.count(); p++) + { + Vector2 point = points[p]; + + for (int i = 0; i < N; i++) + { + float x = dot(axis[i], point); + if (x < mins[i].x) mins[i].x = x; + if (x > maxs[i].x) maxs[i].x = x; + + float y = dot(Vector2(-axis[i].y, axis[i].x), point); + if (y < mins[i].y) mins[i].y = y; + if (y > maxs[i].y) maxs[i].y = y; + } + } + + // Find box with minimum area. + best = -1; + int second_best = -1; + float best_area = FLT_MAX; + float second_best_area = FLT_MAX; + + for (int i = 0; i < N; i++) + { + float area = (maxs[i].x - mins[i].x) * (maxs[i].y - mins[i].y); + + if (area < best_area) + { + second_best_area = best_area; + second_best = best; + + best_area = area; + best = i; + } + else if (area < second_best_area) + { + second_best_area = area; + second_best = i; + } + } + nvDebugCheck(best != -1); + nvDebugCheck(second_best != -1); + nvDebugCheck(best != second_best); + + if (j != iterationCount-1) + { + // Handle wrap-around during the first iteration. + if (j == 0) { + if (best == 0 && second_best == N-1) best = N; + if (best == N-1 && second_best == 0) second_best = N; + } + + if (best < second_best) swap(best, second_best); + + // Update angles. + float deltaAngle = (maxAngle - minAngle) / N; + maxAngle = minAngle + (best - 0.5f) * deltaAngle; + minAngle = minAngle + (second_best + 0.5f) * deltaAngle; + } + } + + // Compute major and minor axis, and origin. + *majorAxis = axis[best]; + *minorAxis = Vector2(-axis[best].y, axis[best].x); + *origin = mins[best]; + + // @@ If the parameterization is invalid, we could have an interior vertex outside the boundary. + // @@ In that case the returned bounding box would be incorrect. Compute updated bounds here. + /*for (uint p = 0; p < points.count(); p++) + { + Vector2 point = points[p]; + + for (int i = 0; i < N; i++) + { + float x = dot(*majorAxis, point); + float y = dot(*minorAxis, point); + } + }*/ +#endif + + return true; +} + + +void AtlasPacker::packCharts(int quality, float texelsPerUnit, bool blockAligned, bool conservative) +{ + const uint chartCount = m_atlas->chartCount(); + if (chartCount == 0) return; + + Array<float> chartOrderArray; + chartOrderArray.resize(chartCount); + + Array<Vector2> chartExtents; + chartExtents.resize(chartCount); + + float meshArea = 0; + for (uint c = 0; c < chartCount; c++) + { + Chart * chart = m_atlas->chartAt(c); + + if (!chart->isVertexMapped() && !chart->isDisk()) { + chartOrderArray[c] = 0; + + // Skip non-disks. + continue; + } + + Vector2 extents(0.0f); + + if (chart->isVertexMapped()) { + // Let's assume vertex maps are arranged in a rectangle. + //HalfEdge::Mesh * mesh = chart->chartMesh(); + + // Arrange vertices in a rectangle. + extents.x = float(chart->vertexMapWidth); + extents.y = float(chart->vertexMapHeight); + } + else { + // Compute surface area to sort charts. + float chartArea = chart->computeSurfaceArea(); + meshArea += chartArea; + //chartOrderArray[c] = chartArea; + + // Compute chart scale + float parametricArea = fabs(chart->computeParametricArea()); // @@ There doesn't seem to be anything preventing parametric area to be negative. + if (parametricArea < NV_EPSILON) { + // When the parametric area is too small we use a rough approximation to prevent divisions by very small numbers. + Vector2 bounds = chart->computeParametricBounds(); + parametricArea = bounds.x * bounds.y; + } + float scale = (chartArea / parametricArea) * texelsPerUnit; + if (parametricArea == 0) // < NV_EPSILON) + { + scale = 0; + } + nvCheck(isFinite(scale)); + + // Compute bounding box of chart. + Vector2 majorAxis, minorAxis, origin, end; + if (!computeBoundingBox(chart, &majorAxis, &minorAxis, &origin, &end)) { + m_atlas->setFailed(); + return; + } + + nvCheck(isFinite(majorAxis) && isFinite(minorAxis) && isFinite(origin)); + + // Sort charts by perimeter. @@ This is sometimes producing somewhat unexpected results. Is this right? + //chartOrderArray[c] = ((end.x - origin.x) + (end.y - origin.y)) * scale; + + // Translate, rotate and scale vertices. Compute extents. + HalfEdge::Mesh * mesh = chart->chartMesh(); + const uint vertexCount = mesh->vertexCount(); + for (uint i = 0; i < vertexCount; i++) + { + HalfEdge::Vertex * vertex = mesh->vertexAt(i); + + //Vector2 t = vertex->tex - origin; + Vector2 tmp; + tmp.x = dot(vertex->tex, majorAxis); + tmp.y = dot(vertex->tex, minorAxis); + tmp -= origin; + tmp *= scale; + if (tmp.x < 0 || tmp.y < 0) { + nvDebug("tmp: %f %f\n", tmp.x, tmp.y); + nvDebug("scale: %f\n", scale); + nvDebug("origin: %f %f\n", origin.x, origin.y); + nvDebug("majorAxis: %f %f\n", majorAxis.x, majorAxis.y); + nvDebug("minorAxis: %f %f\n", minorAxis.x, minorAxis.y); + nvDebugBreak(); + } + //nvCheck(tmp.x >= 0 && tmp.y >= 0); + + vertex->tex = tmp; + + nvCheck(isFinite(vertex->tex.x) && isFinite(vertex->tex.y)); + + extents = max(extents, tmp); + } + nvDebugCheck(extents.x >= 0 && extents.y >= 0); + + // Limit chart size. + if (extents.x > 1024 || extents.y > 1024) { + float limit = max(extents.x, extents.y); + + scale = 1024 / (limit + 1); + + for (uint i = 0; i < vertexCount; i++) + { + HalfEdge::Vertex * vertex = mesh->vertexAt(i); + vertex->tex *= scale; + } + + extents *= scale; + + nvDebugCheck(extents.x <= 1024 && extents.y <= 1024); + } + + + // Scale the charts to use the entire texel area available. So, if the width is 0.1 we could scale it to 1 without increasing the lightmap usage and making a better + // use of it. In many cases this also improves the look of the seams, since vertices on the chart boundaries have more chances of being aligned with the texel centers. + + float scale_x = 1.0f; + float scale_y = 1.0f; + + float divide_x = 1.0f; + float divide_y = 1.0f; + + if (extents.x > 0) { + int cw = ftoi_ceil(extents.x); + + if (blockAligned) { + // Align all chart extents to 4x4 blocks, but taking padding into account. + if (conservative) { + cw = align(cw + 2, 4) - 2; + } + else { + cw = align(cw + 1, 4) - 1; + } + } + + scale_x = (float(cw) - NV_EPSILON); + divide_x = extents.x; + extents.x = float(cw); + } + + if (extents.y > 0) { + int ch = ftoi_ceil(extents.y); + + if (blockAligned) { + // Align all chart extents to 4x4 blocks, but taking padding into account. + if (conservative) { + ch = align(ch + 2, 4) - 2; + } + else { + ch = align(ch + 1, 4) - 1; + } + } + + scale_y = (float(ch) - NV_EPSILON); + divide_y = extents.y; + extents.y = float(ch); + } + + for (uint v = 0; v < vertexCount; v++) { + HalfEdge::Vertex * vertex = mesh->vertexAt(v); + + vertex->tex.x /= divide_x; + vertex->tex.y /= divide_y; + vertex->tex.x *= scale_x; + vertex->tex.y *= scale_y; + + nvCheck(isFinite(vertex->tex.x) && isFinite(vertex->tex.y)); + } + } + + chartExtents[c] = extents; + + // Sort charts by perimeter. + chartOrderArray[c] = extents.x + extents.y; + } + + // @@ We can try to improve compression of small charts by sorting them by proximity like we do with vertex samples. + // @@ How to do that? One idea: compute chart centroid, insert into grid, compute morton index of the cell, sort based on morton index. + // @@ We would sort by morton index, first, then quantize the chart sizes, so that all small charts have the same size, and sort by size preserving the morton order. + + //nvDebug("Sorting charts.\n"); + + // Sort charts by area. + m_radix.sort(chartOrderArray); + const uint32 * ranks = m_radix.ranks(); + + // Estimate size of the map based on the mesh surface area and given texel scale. + float texelCount = meshArea * square(texelsPerUnit) / 0.75f; // Assume 75% utilization. + if (texelCount < 1) texelCount = 1; + uint approximateExtent = nextPowerOfTwo(uint(sqrtf(texelCount))); + + //nvDebug("Init bitmap.\n"); + + // @@ Pack all charts smaller than a texel into a compact rectangle. + // @@ Start considering only 1x1 charts. Extend to 1xn charts later. + + /*for (uint i = 0; i < chartCount; i++) + { + uint c = ranks[chartCount - i - 1]; // largest chart first + + Chart * chart = m_atlas->chartAt(c); + + if (!chart->isDisk()) continue; + + if (iceil(chartExtents[c].x) == 1 && iceil(chartExtents[c].x) == 1) { + // @@ Add to + } + }*/ + + + + // Init bit map. + m_bitmap.clearAll(); + if (approximateExtent > m_bitmap.width()) { + m_bitmap.resize(approximateExtent, approximateExtent, false); + m_debug_bitmap.resize(approximateExtent, approximateExtent); + m_debug_bitmap.fill(Color32(0,0,0,0)); + } + + + int w = 0; + int h = 0; + +#if 1 + // Add sorted charts to bitmap. + for (uint i = 0; i < chartCount; i++) + { + uint c = ranks[chartCount - i - 1]; // largest chart first + + Chart * chart = m_atlas->chartAt(c); + + if (!chart->isVertexMapped() && !chart->isDisk()) continue; + + //float scale_x = 1; + //float scale_y = 1; + + BitMap chart_bitmap; + + if (chart->isVertexMapped()) { + // Init all bits to 1. + chart_bitmap.resize(ftoi_ceil(chartExtents[c].x), ftoi_ceil(chartExtents[c].y), /*initValue=*/true); + + // @@ Another alternative would be to try to map each vertex to a different texel trying to fill all the available unused texels. + } + else { + // @@ Add special cases for dot and line charts. @@ Lightmap rasterizer also needs to handle these special cases. + // @@ We could also have a special case for chart quads. If the quad surface <= 4 texels, align vertices with texel centers and do not add padding. May be very useful for foliage. + + // @@ In general we could reduce the padding of all charts by one texel by using a rasterizer that takes into account the 2-texel footprint of the tent bilinear filter. For example, + // if we have a chart that is less than 1 texel wide currently we add one texel to the left and one texel to the right creating a 3-texel-wide bitmap. However, if we know that the + // chart is only 1 texel wide we could align it so that it only touches the footprint of two texels: + + // | | <- Touches texels 0, 1 and 2. + // | | <- Only touches texels 0 and 1. + // \ \ / \ / / + // \ X X / + // \ / \ / \ / + // V V V + // 0 1 2 + + if (conservative) { + // Init all bits to 0. + chart_bitmap.resize(ftoi_ceil(chartExtents[c].x) + 2, ftoi_ceil(chartExtents[c].y) + 2, /*initValue=*/false); // + 2 to add padding on both sides. + + // Rasterize chart and dilate. + drawChartBitmapDilate(chart, &chart_bitmap, /*padding=*/1); + } + else { + // Init all bits to 0. + chart_bitmap.resize(ftoi_ceil(chartExtents[c].x) + 1, ftoi_ceil(chartExtents[c].y) + 1, /*initValue=*/false); // Add half a texels on each side. + + // Rasterize chart and dilate. + drawChartBitmap(chart, &chart_bitmap, Vector2(1), Vector2(0.5)); + } + } + + int best_x, best_y; + int best_cw, best_ch; // Includes padding now. + int best_r; + findChartLocation(quality, &chart_bitmap, chartExtents[c], w, h, &best_x, &best_y, &best_cw, &best_ch, &best_r); + + /*if (w < best_x + best_cw || h < best_y + best_ch) + { + nvDebug("Resize extents to (%d, %d).\n", best_x + best_cw, best_y + best_ch); + }*/ + + // Update parametric extents. + w = max(w, best_x + best_cw); + h = max(h, best_y + best_ch); + + w = align(w, 4); + h = align(h, 4); + + // Resize bitmap if necessary. + if (uint(w) > m_bitmap.width() || uint(h) > m_bitmap.height()) + { + //nvDebug("Resize bitmap (%d, %d).\n", nextPowerOfTwo(w), nextPowerOfTwo(h)); + m_bitmap.resize(nextPowerOfTwo(U32(w)), nextPowerOfTwo(U32(h)), false); + m_debug_bitmap.resize(nextPowerOfTwo(U32(w)), nextPowerOfTwo(U32(h))); + } + + //nvDebug("Add chart at (%d, %d).\n", best_x, best_y); + + addChart(&chart_bitmap, w, h, best_x, best_y, best_r, /*debugOutput=*/NULL); + + // IC: Output chart again to debug bitmap. + if (chart->isVertexMapped()) { + addChart(&chart_bitmap, w, h, best_x, best_y, best_r, &m_debug_bitmap); + } + else { + addChart(chart, w, h, best_x, best_y, best_r, &m_debug_bitmap); + } + + //float best_angle = 2 * PI * best_r; + + // Translate and rotate chart texture coordinates. + HalfEdge::Mesh * mesh = chart->chartMesh(); + const uint vertexCount = mesh->vertexCount(); + for (uint v = 0; v < vertexCount; v++) + { + HalfEdge::Vertex * vertex = mesh->vertexAt(v); + + Vector2 t = vertex->tex; + if (best_r) swap(t.x, t.y); + //vertex->tex.x = best_x + t.x * cosf(best_angle) - t.y * sinf(best_angle); + //vertex->tex.y = best_y + t.x * sinf(best_angle) + t.y * cosf(best_angle); + + vertex->tex.x = best_x + t.x + 0.5f; + vertex->tex.y = best_y + t.y + 0.5f; + + nvCheck(vertex->tex.x >= 0 && vertex->tex.y >= 0); + nvCheck(isFinite(vertex->tex.x) && isFinite(vertex->tex.y)); + } + +#if DEBUG_OUTPUT && 0 + StringBuilder fileName; + fileName.format("debug_packer_%d.tga", i); + //outputDebugBitmap(fileName.str(), m_bitmap, w, h); + outputDebugImage(fileName.str(), m_debug_bitmap, w, h); +#endif + } + +#else // 0 + + // Add sorted charts to bitmap. + for (uint i = 0; i < chartCount; i++) + { + uint c = ranks[chartCount - i - 1]; // largest chart first + + Chart * chart = m_atlas->chartAt(c); + + if (!chart->isDisk()) continue; + + Vector2 scale(1, 1); + +#if 0 // old method. + //m_padding_x = 2*padding; + //m_padding_y = 2*padding; +#else + //m_padding_x = 0; //padding; + //m_padding_y = 0; //padding; +#endif + + int bw = ftoi_ceil(chartExtents[c].x + 1); + int bh = ftoi_ceil(chartExtents[c].y + 1); + + if (chartExtents[c].x < 1.0f) { + scale.x = 0.01f; // @@ Ideally we would like to scale it to 0, but then our rasterizer would not touch any pixels. + bw = 1; + } + if (chartExtents[c].y < 1.0f) { + scale.y = 0.01f; + bh = 1; + } + + //BitMap chart_bitmap(iceil(chartExtents[c].x) + 1 + m_padding_x * 2, iceil(chartExtents[c].y) + 1 + m_padding_y * 2); + //BitMap chart_bitmap(ftoi_ceil(chartExtents[c].x/2)*2, ftoi_ceil(chartExtents[c].y/2)*2); + BitMap chart_bitmap(bw, bh); + chart_bitmap.clearAll(); + + Vector2 offset; + offset.x = 0; // (chart_bitmap.width() - chartExtents[c].x) * 0.5f; + offset.y = 0; // (chart_bitmap.height() - chartExtents[c].y) * 0.5f; + + drawChartBitmap(chart, &chart_bitmap, scale, offset); + + int best_x, best_y; + int best_cw, best_ch; + int best_r; + findChartLocation(quality, &chart_bitmap, chartExtents[c], w, h, &best_x, &best_y, &best_cw, &best_ch, &best_r); + + /*if (w < best_x + best_cw || h < best_y + best_ch) + { + nvDebug("Resize extents to (%d, %d).\n", best_x + best_cw, best_y + best_ch); + }*/ + + // Update parametric extents. + w = max(w, best_x + best_cw); + h = max(h, best_y + best_ch); + + // Resize bitmap if necessary. + if (uint(w) > m_bitmap.width() || uint(h) > m_bitmap.height()) + { + //nvDebug("Resize bitmap (%d, %d).\n", nextPowerOfTwo(w), nextPowerOfTwo(h)); + m_bitmap.resize(nextPowerOfTwo(w), nextPowerOfTwo(h), false); + m_debug_bitmap.resize(nextPowerOfTwo(w), nextPowerOfTwo(h)); + } + + //nvDebug("Add chart at (%d, %d).\n", best_x, best_y); + +#if 0 // old method. +#if _DEBUG + checkCanAddChart(chart, w, h, best_x, best_y, best_r); +#endif + + // Add chart. + addChart(chart, w, h, best_x, best_y, best_r); +#else + // Add chart reusing its bitmap. + addChart(&chart_bitmap, w, h, best_x, best_y, best_r); +#endif + + //float best_angle = 2 * PI * best_r; + + // Translate and rotate chart texture coordinates. + HalfEdge::Mesh * mesh = chart->chartMesh(); + const uint vertexCount = mesh->vertexCount(); + for (uint v = 0; v < vertexCount; v++) + { + HalfEdge::Vertex * vertex = mesh->vertexAt(v); + + Vector2 t = vertex->tex * scale + offset; + if (best_r) swap(t.x, t.y); + //vertex->tex.x = best_x + t.x * cosf(best_angle) - t.y * sinf(best_angle); + //vertex->tex.y = best_y + t.x * sinf(best_angle) + t.y * cosf(best_angle); + vertex->tex.x = best_x + t.x + 0.5f; + vertex->tex.y = best_y + t.y + 0.5f; + + nvCheck(vertex->tex.x >= 0 && vertex->tex.y >= 0); + } + +#if DEBUG_OUTPUT && 0 + StringBuilder fileName; + fileName.format("debug_packer_%d.tga", i); + //outputDebugBitmap(fileName.str(), m_bitmap, w, h); + outputDebugImage(fileName.str(), m_debug_bitmap, w, h); +#endif + } + +#endif // 0 + + //w -= padding - 1; // Leave one pixel border! + //h -= padding - 1; + + m_width = max(0, w); + m_height = max(0, h); + + nvCheck(isAligned(m_width, 4)); + nvCheck(isAligned(m_height, 4)); + + m_debug_bitmap.resize(m_width, m_height); + m_debug_bitmap.setFormat(Image::Format_ARGB); + +#if DEBUG_OUTPUT + //outputDebugBitmap("debug_packer_final.tga", m_bitmap, w, h); + //outputDebugImage("debug_packer_final.tga", m_debug_bitmap, w, h); + ImageIO::save("debug_packer_final.tga", &m_debug_bitmap); +#endif +} + + +// IC: Brute force is slow, and random may take too much time to converge. We start inserting large charts in a small atlas. Using brute force is lame, because most of the space +// is occupied at this point. At the end we have many small charts and a large atlas with sparse holes. Finding those holes randomly is slow. A better approach would be to +// start stacking large charts as if they were tetris pieces. Once charts get small try to place them randomly. It may be interesting to try a intermediate strategy, first try +// along one axis and then try exhaustively along that axis. +void AtlasPacker::findChartLocation(int quality, const BitMap * bitmap, Vector2::Arg extents, int w, int h, int * best_x, int * best_y, int * best_w, int * best_h, int * best_r) +{ + int attempts = 256; + if (quality == 1) attempts = 4096; + if (quality == 2) attempts = 2048; + if (quality == 3) attempts = 1024; + if (quality == 4) attempts = 512; + + if (quality == 0 || w*h < attempts) + { + findChartLocation_bruteForce(bitmap, extents, w, h, best_x, best_y, best_w, best_h, best_r); + } + else + { + findChartLocation_random(bitmap, extents, w, h, best_x, best_y, best_w, best_h, best_r, attempts); + } +} + +#define BLOCK_SIZE 4 + +void AtlasPacker::findChartLocation_bruteForce(const BitMap * bitmap, Vector2::Arg extents, int w, int h, int * best_x, int * best_y, int * best_w, int * best_h, int * best_r) +{ + int best_metric = INT_MAX; + + // Try two different orientations. + for (int r = 0; r < 2; r++) + { + int cw = bitmap->width(); + int ch = bitmap->height(); + if (r & 1) swap(cw, ch); + + for (int y = 0; y <= h + 1; y += BLOCK_SIZE) // + 1 to extend atlas in case atlas full. + { + for (int x = 0; x <= w + 1; x += BLOCK_SIZE) // + 1 not really necessary here. + { + // Early out. + int area = max(w, x+cw) * max(h, y+ch); + //int perimeter = max(w, x+cw) + max(h, y+ch); + int extents = max(max(w, x+cw), max(h, y+ch)); + + int metric = extents*extents + area; + + if (metric > best_metric) { + continue; + } + if (metric == best_metric && max(x, y) >= max(*best_x, *best_y)) { + // If metric is the same, pick the one closest to the origin. + continue; + } + + if (canAddChart(bitmap, w, h, x, y, r)) + { + best_metric = metric; + *best_x = x; + *best_y = y; + *best_w = cw; + *best_h = ch; + *best_r = r; + + if (area == w*h) + { + // Chart is completely inside, do not look at any other location. + goto done; + } + } + } + } + } + +done: + nvDebugCheck (best_metric != INT_MAX); +} + + +void AtlasPacker::findChartLocation_random(const BitMap * bitmap, Vector2::Arg extents, int w, int h, int * best_x, int * best_y, int * best_w, int * best_h, int * best_r, int minTrialCount) +{ + int best_metric = INT_MAX; + + for (int i = 0; i < minTrialCount || best_metric == INT_MAX; i++) + { + int r = m_rand.getRange(1); + int x = m_rand.getRange(w + 1); // + 1 to extend atlas in case atlas full. We may want to use a higher number to increase probability of extending atlas. + int y = m_rand.getRange(h + 1); // + 1 to extend atlas in case atlas full. + + x = align(x, BLOCK_SIZE); + y = align(y, BLOCK_SIZE); + + int cw = bitmap->width(); + int ch = bitmap->height(); + if (r & 1) swap(cw, ch); + + // Early out. + int area = max(w, x+cw) * max(h, y+ch); + //int perimeter = max(w, x+cw) + max(h, y+ch); + int extents = max(max(w, x+cw), max(h, y+ch)); + + int metric = extents*extents + area; + + if (metric > best_metric) { + continue; + } + if (metric == best_metric && min(x, y) > min(*best_x, *best_y)) { + // If metric is the same, pick the one closest to the origin. + continue; + } + + if (canAddChart(bitmap, w, h, x, y, r)) + { + best_metric = metric; + *best_x = x; + *best_y = y; + *best_w = cw; + *best_h = ch; + *best_r = r; + + if (area == w*h) + { + // Chart is completely inside, do not look at any other location. + break; + } + } + } +} + + +void AtlasPacker::drawChartBitmapDilate(const Chart * chart, BitMap * bitmap, int padding) +{ + const int w = bitmap->width(); + const int h = bitmap->height(); + const Vector2 extents = Vector2(float(w), float(h)); + + // Rasterize chart faces, check that all bits are not set. + const uint faceCount = chart->faceCount(); + for (uint f = 0; f < faceCount; f++) + { + const HalfEdge::Face * face = chart->chartMesh()->faceAt(f); + + Vector2 vertices[4]; + + uint edgeCount = 0; + for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + if (edgeCount < 4) + { + vertices[edgeCount] = it.vertex()->tex + Vector2(0.5) + Vector2(float(padding), float(padding)); + } + edgeCount++; + } + + if (edgeCount == 3) + { + Raster::drawTriangle(Raster::Mode_Antialiased, extents, true, vertices, AtlasPacker::setBitsCallback, bitmap); + } + else + { + Raster::drawQuad(Raster::Mode_Antialiased, extents, true, vertices, AtlasPacker::setBitsCallback, bitmap); + } + } + + // Expand chart by padding pixels. (dilation) + BitMap tmp(w, h); + for (int i = 0; i < padding; i++) { + tmp.clearAll(); + + for (int y = 0; y < h; y++) { + for (int x = 0; x < w; x++) { + bool b = bitmap->bitAt(x, y); + if (!b) { + if (x > 0) { + b |= bitmap->bitAt(x - 1, y); + if (y > 0) b |= bitmap->bitAt(x - 1, y - 1); + if (y < h-1) b |= bitmap->bitAt(x - 1, y + 1); + } + if (y > 0) b |= bitmap->bitAt(x, y - 1); + if (y < h-1) b |= bitmap->bitAt(x, y + 1); + if (x < w-1) { + b |= bitmap->bitAt(x + 1, y); + if (y > 0) b |= bitmap->bitAt(x + 1, y - 1); + if (y < h-1) b |= bitmap->bitAt(x + 1, y + 1); + } + } + if (b) tmp.setBitAt(x, y); + } + } + + swap(tmp, *bitmap); + } +} + + +void AtlasPacker::drawChartBitmap(const Chart * chart, BitMap * bitmap, const Vector2 & scale, const Vector2 & offset) +{ + const int w = bitmap->width(); + const int h = bitmap->height(); + const Vector2 extents = Vector2(float(w), float(h)); + + static const Vector2 pad[4] = { + Vector2(-0.5, -0.5), + Vector2(0.5, -0.5), + Vector2(-0.5, 0.5), + Vector2(0.5, 0.5) + }; + /*static const Vector2 pad[4] = { + Vector2(-1, -1), + Vector2(1, -1), + Vector2(-1, 1), + Vector2(1, 1) + };*/ + + // Rasterize 4 times to add proper padding. + for (int i = 0; i < 4; i++) { + + // Rasterize chart faces, check that all bits are not set. + const uint faceCount = chart->chartMesh()->faceCount(); + for (uint f = 0; f < faceCount; f++) + { + const HalfEdge::Face * face = chart->chartMesh()->faceAt(f); + + Vector2 vertices[4]; + + uint edgeCount = 0; + for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + if (edgeCount < 4) + { + vertices[edgeCount] = it.vertex()->tex * scale + offset + pad[i]; + nvCheck(ftoi_ceil(vertices[edgeCount].x) >= 0); + nvCheck(ftoi_ceil(vertices[edgeCount].y) >= 0); + nvCheck(ftoi_ceil(vertices[edgeCount].x) <= w); + nvCheck(ftoi_ceil(vertices[edgeCount].y) <= h); + } + edgeCount++; + } + + if (edgeCount == 3) + { + Raster::drawTriangle(Raster::Mode_Antialiased, extents, /*enableScissors=*/true, vertices, AtlasPacker::setBitsCallback, bitmap); + } + else + { + Raster::drawQuad(Raster::Mode_Antialiased, extents, /*enableScissors=*/true, vertices, AtlasPacker::setBitsCallback, bitmap); + } + } + } + + // @@ This only allows us to expand the size in texel intervals. + /*if (m_padding_x != 0 && m_padding_y != 0)*/ { + + // Expand chart by padding pixels. (dilation) + BitMap tmp(w, h); + //for (int i = 0; i < 1; i++) { + tmp.clearAll(); + + for (int y = 0; y < h; y++) { + for (int x = 0; x < w; x++) { + bool b = bitmap->bitAt(x, y); + if (!b) { + if (x > 0) { + b |= bitmap->bitAt(x - 1, y); + if (y > 0) b |= bitmap->bitAt(x - 1, y - 1); + if (y < h-1) b |= bitmap->bitAt(x - 1, y + 1); + } + if (y > 0) b |= bitmap->bitAt(x, y - 1); + if (y < h-1) b |= bitmap->bitAt(x, y + 1); + if (x < w-1) { + b |= bitmap->bitAt(x + 1, y); + if (y > 0) b |= bitmap->bitAt(x + 1, y - 1); + if (y < h-1) b |= bitmap->bitAt(x + 1, y + 1); + } + } + if (b) tmp.setBitAt(x, y); + } + } + + swap(tmp, *bitmap); + //} + } +} + +bool AtlasPacker::canAddChart(const BitMap * bitmap, int atlas_w, int atlas_h, int offset_x, int offset_y, int r) +{ + nvDebugCheck(r == 0 || r == 1); + + // Check whether the two bitmaps overlap. + + const int w = bitmap->width(); + const int h = bitmap->height(); + + if (r == 0) { + for (int y = 0; y < h; y++) { + int yy = y + offset_y; + if (yy >= 0) { + for (int x = 0; x < w; x++) { + int xx = x + offset_x; + if (xx >= 0) { + if (bitmap->bitAt(x, y)) { + if (xx < atlas_w && yy < atlas_h) { + if (m_bitmap.bitAt(xx, yy)) return false; + } + } + } + } + } + } + } + else if (r == 1) { + for (int y = 0; y < h; y++) { + int xx = y + offset_x; + if (xx >= 0) { + for (int x = 0; x < w; x++) { + int yy = x + offset_y; + if (yy >= 0) { + if (bitmap->bitAt(x, y)) { + if (xx < atlas_w && yy < atlas_h) { + if (m_bitmap.bitAt(xx, yy)) return false; + } + } + } + } + } + } + } + + return true; +} + +#if 0 +void AtlasPacker::checkCanAddChart(const Chart * chart, int w, int h, int x, int y, int r) +{ + nvDebugCheck(r == 0 || r == 1); + Vector2 extents = Vector2(float(w), float(h)); + Vector2 offset = Vector2(float(x), float(y)); + + // Rasterize chart faces, set bits. + const uint faceCount = chart->faceCount(); + for (uint f = 0; f < faceCount; f++) + { + const HalfEdge::Face * face = chart->chartMesh()->faceAt(f); + + Vector2 vertices[4]; + + uint edgeCount = 0; + for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + if (edgeCount < 4) + { + Vector2 t = it.vertex()->tex; + if (r == 1) swap(t.x, t.y); + vertices[edgeCount] = t + offset; + } + edgeCount++; + } + + if (edgeCount == 3) + { + Raster::drawTriangle(Raster::Mode_Antialiased, extents, /*enableScissors=*/true, vertices, AtlasPacker::checkBitsCallback, &m_bitmap); + } + else + { + Raster::drawQuad(Raster::Mode_Antialiased, extents, /*enableScissors=*/true, vertices, AtlasPacker::checkBitsCallback, &m_bitmap); + } + } +} +#endif // 0 + + +static Color32 chartColor = Color32(0); +static void selectRandomColor(MTRand & rand) { + // Pick random color for this chart. @@ Select random hue, but fixed saturation/luminance? + chartColor.r = 128 + rand.getRange(127); + chartColor.g = 128 + rand.getRange(127); + chartColor.b = 128 + rand.getRange(127); + chartColor.a = 255; +} +static bool debugDrawCallback(void * param, int x, int y, Vector3::Arg, Vector3::Arg, Vector3::Arg, float area) +{ + Image * image = (Image *)param; + + if (area > 0.0) { + Color32 c = image->pixel(x, y); + c.r = chartColor.r; + c.g = chartColor.g; + c.b = chartColor.b; + c.a += U8(ftoi_round(0.5f * area * 255)); + image->pixel(x, y) = c; + } + + return true; +} + +void AtlasPacker::addChart(const Chart * chart, int w, int h, int x, int y, int r, Image * debugOutput) +{ + nvDebugCheck(r == 0 || r == 1); + + nvDebugCheck(debugOutput != NULL); + selectRandomColor(m_rand); + + Vector2 extents = Vector2(float(w), float(h)); + Vector2 offset = Vector2(float(x), float(y)) + Vector2(0.5); + + // Rasterize chart faces, set bits. + const uint faceCount = chart->faceCount(); + for (uint f = 0; f < faceCount; f++) + { + const HalfEdge::Face * face = chart->chartMesh()->faceAt(f); + + Vector2 vertices[4]; + + uint edgeCount = 0; + for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + if (edgeCount < 4) + { + Vector2 t = it.vertex()->tex; + if (r == 1) swap(t.x, t.y); + vertices[edgeCount] = t + offset; + } + edgeCount++; + } + + if (edgeCount == 3) + { + Raster::drawTriangle(Raster::Mode_Antialiased, extents, /*enableScissors=*/true, vertices, debugDrawCallback, debugOutput); + } + else + { + Raster::drawQuad(Raster::Mode_Antialiased, extents, /*enableScissors=*/true, vertices, debugDrawCallback, debugOutput); + } + } +} + + +void AtlasPacker::addChart(const BitMap * bitmap, int atlas_w, int atlas_h, int offset_x, int offset_y, int r, Image * debugOutput) +{ + nvDebugCheck(r == 0 || r == 1); + + // Check whether the two bitmaps overlap. + + const int w = bitmap->width(); + const int h = bitmap->height(); + + if (debugOutput != NULL) { + selectRandomColor(m_rand); + } + + if (r == 0) { + for (int y = 0; y < h; y++) { + int yy = y + offset_y; + if (yy >= 0) { + for (int x = 0; x < w; x++) { + int xx = x + offset_x; + if (xx >= 0) { + if (bitmap->bitAt(x, y)) { + if (xx < atlas_w && yy < atlas_h) { + if (debugOutput) debugOutput->pixel(xx, yy) = chartColor; + else { + nvDebugCheck(m_bitmap.bitAt(xx, yy) == false); + m_bitmap.setBitAt(xx, yy); + } + } + } + } + } + } + } + } + else if (r == 1) { + for (int y = 0; y < h; y++) { + int xx = y + offset_x; + if (xx >= 0) { + for (int x = 0; x < w; x++) { + int yy = x + offset_y; + if (yy >= 0) { + if (bitmap->bitAt(x, y)) { + if (xx < atlas_w && yy < atlas_h) { + if (debugOutput) debugOutput->pixel(xx, yy) = chartColor; + else { + nvDebugCheck(m_bitmap.bitAt(xx, yy) == false); + m_bitmap.setBitAt(xx, yy); + } + } + } + } + } + } + } + } +} + + + +/*static*/ bool AtlasPacker::checkBitsCallback(void * param, int x, int y, Vector3::Arg, Vector3::Arg, Vector3::Arg, float) +{ + BitMap * bitmap = (BitMap * )param; + + nvDebugCheck(bitmap->bitAt(x, y) == false); + + return true; +} + +/*static*/ bool AtlasPacker::setBitsCallback(void * param, int x, int y, Vector3::Arg, Vector3::Arg, Vector3::Arg, float area) +{ + BitMap * bitmap = (BitMap * )param; + + if (area > 0.0) { + bitmap->setBitAt(x, y); + } + + return true; +} + + + +float AtlasPacker::computeAtlasUtilization() const { + const uint w = m_width; + const uint h = m_height; + nvDebugCheck(w <= m_bitmap.width()); + nvDebugCheck(h <= m_bitmap.height()); + + uint count = 0; + for (uint y = 0; y < h; y++) { + for (uint x = 0; x < w; x++) { + count += m_bitmap.bitAt(x, y); + } + } + + return float(count) / (w * h); +} diff --git a/thirdparty/thekla_atlas/nvmesh/param/AtlasPacker.h b/thirdparty/thekla_atlas/nvmesh/param/AtlasPacker.h new file mode 100644 index 0000000000..2d305f38cd --- /dev/null +++ b/thirdparty/thekla_atlas/nvmesh/param/AtlasPacker.h @@ -0,0 +1,63 @@ +// This code is in the public domain -- castano@gmail.com + +#pragma once +#ifndef NV_MESH_ATLASPACKER_H +#define NV_MESH_ATLASPACKER_H + +#include "nvcore/RadixSort.h" +#include "nvmath/Vector.h" +#include "nvmath/Random.h" +#include "nvimage/BitMap.h" +#include "nvimage/Image.h" + +#include "nvmesh/nvmesh.h" + + +namespace nv +{ + class Atlas; + class Chart; + + struct AtlasPacker + { + AtlasPacker(Atlas * atlas); + ~AtlasPacker(); + + void packCharts(int quality, float texelArea, bool blockAligned, bool conservative); + float computeAtlasUtilization() const; + + private: + + void findChartLocation(int quality, const BitMap * bitmap, Vector2::Arg extents, int w, int h, int * best_x, int * best_y, int * best_w, int * best_h, int * best_r); + void findChartLocation_bruteForce(const BitMap * bitmap, Vector2::Arg extents, int w, int h, int * best_x, int * best_y, int * best_w, int * best_h, int * best_r); + void findChartLocation_random(const BitMap * bitmap, Vector2::Arg extents, int w, int h, int * best_x, int * best_y, int * best_w, int * best_h, int * best_r, int minTrialCount); + + void drawChartBitmapDilate(const Chart * chart, BitMap * bitmap, int padding); + void drawChartBitmap(const Chart * chart, BitMap * bitmap, const Vector2 & scale, const Vector2 & offset); + + bool canAddChart(const BitMap * bitmap, int w, int h, int x, int y, int r); + void addChart(const BitMap * bitmap, int w, int h, int x, int y, int r, Image * debugOutput); + //void checkCanAddChart(const Chart * chart, int w, int h, int x, int y, int r); + void addChart(const Chart * chart, int w, int h, int x, int y, int r, Image * debugOutput); + + + static bool checkBitsCallback(void * param, int x, int y, Vector3::Arg bar, Vector3::Arg dx, Vector3::Arg dy, float coverage); + static bool setBitsCallback(void * param, int x, int y, Vector3::Arg bar, Vector3::Arg dx, Vector3::Arg dy, float coverage); + + private: + + Atlas * m_atlas; + BitMap m_bitmap; + Image m_debug_bitmap; + RadixSort m_radix; + + uint m_width; + uint m_height; + + MTRand m_rand; + + }; + +} // nv namespace + +#endif // NV_MESH_ATLASPACKER_H diff --git a/thirdparty/thekla_atlas/nvmesh/param/LeastSquaresConformalMap.cpp b/thirdparty/thekla_atlas/nvmesh/param/LeastSquaresConformalMap.cpp new file mode 100644 index 0000000000..cd1e8bbb7b --- /dev/null +++ b/thirdparty/thekla_atlas/nvmesh/param/LeastSquaresConformalMap.cpp @@ -0,0 +1,483 @@ +// Copyright NVIDIA Corporation 2008 -- Ignacio Castano <icastano@nvidia.com> + +#include "nvmesh.h" // pch + +#include "LeastSquaresConformalMap.h" +#include "ParameterizationQuality.h" +#include "Util.h" + +#include "nvmesh/halfedge/Mesh.h" +#include "nvmesh/halfedge/Vertex.h" +#include "nvmesh/halfedge/Face.h" + +#include "nvmath/Sparse.h" +#include "nvmath/Solver.h" +#include "nvmath/Vector.inl" + +#include "nvcore/Array.inl" + + +using namespace nv; +using namespace HalfEdge; + +namespace +{ + + // Test all pairs of vertices in the boundary and check distance. + static void findDiameterVertices(HalfEdge::Mesh * mesh, HalfEdge::Vertex ** a, HalfEdge::Vertex ** b) + { + nvDebugCheck(mesh != NULL); + nvDebugCheck(a != NULL); + nvDebugCheck(b != NULL); + + const uint vertexCount = mesh->vertexCount(); + + float maxLength = 0.0f; + + for (uint v0 = 1; v0 < vertexCount; v0++) + { + HalfEdge::Vertex * vertex0 = mesh->vertexAt(v0); + nvDebugCheck(vertex0 != NULL); + + if (!vertex0->isBoundary()) continue; + + for (uint v1 = 0; v1 < v0; v1++) + { + HalfEdge::Vertex * vertex1 = mesh->vertexAt(v1); + nvDebugCheck(vertex1 != NULL); + + if (!vertex1->isBoundary()) continue; + + float len = length(vertex0->pos - vertex1->pos); + + if (len > maxLength) + { + maxLength = len; + + *a = vertex0; + *b = vertex1; + } + } + } + + nvDebugCheck(*a != NULL && *b != NULL); + } + + // Fast sweep in 3 directions + static bool findApproximateDiameterVertices(HalfEdge::Mesh * mesh, HalfEdge::Vertex ** a, HalfEdge::Vertex ** b) + { + nvDebugCheck(mesh != NULL); + nvDebugCheck(a != NULL); + nvDebugCheck(b != NULL); + + const uint vertexCount = mesh->vertexCount(); + + HalfEdge::Vertex * minVertex[3]; + HalfEdge::Vertex * maxVertex[3]; + + minVertex[0] = minVertex[1] = minVertex[2] = NULL; + maxVertex[0] = maxVertex[1] = maxVertex[2] = NULL; + + for (uint v = 1; v < vertexCount; v++) + { + HalfEdge::Vertex * vertex = mesh->vertexAt(v); + nvDebugCheck(vertex != NULL); + + if (vertex->isBoundary()) + { + minVertex[0] = minVertex[1] = minVertex[2] = vertex; + maxVertex[0] = maxVertex[1] = maxVertex[2] = vertex; + break; + } + } + + if (minVertex[0] == NULL) + { + // Input mesh has not boundaries. + return false; + } + + for (uint v = 1; v < vertexCount; v++) + { + HalfEdge::Vertex * vertex = mesh->vertexAt(v); + nvDebugCheck(vertex != NULL); + + if (!vertex->isBoundary()) + { + // Skip interior vertices. + continue; + } + + if (vertex->pos.x < minVertex[0]->pos.x) minVertex[0] = vertex; + else if (vertex->pos.x > maxVertex[0]->pos.x) maxVertex[0] = vertex; + + if (vertex->pos.y < minVertex[1]->pos.y) minVertex[1] = vertex; + else if (vertex->pos.y > maxVertex[1]->pos.y) maxVertex[1] = vertex; + + if (vertex->pos.z < minVertex[2]->pos.z) minVertex[2] = vertex; + else if (vertex->pos.z > maxVertex[2]->pos.z) maxVertex[2] = vertex; + } + + float lengths[3]; + for (int i = 0; i < 3; i++) + { + lengths[i] = length(minVertex[i]->pos - maxVertex[i]->pos); + } + + if (lengths[0] > lengths[1] && lengths[0] > lengths[2]) + { + *a = minVertex[0]; + *b = maxVertex[0]; + } + else if (lengths[1] > lengths[2]) + { + *a = minVertex[1]; + *b = maxVertex[1]; + } + else + { + *a = minVertex[2]; + *b = maxVertex[2]; + } + + return true; + } + + // Conformal relations from Bruno Levy: + + // Computes the coordinates of the vertices of a triangle + // in a local 2D orthonormal basis of the triangle's plane. + static void project_triangle(Vector3::Arg p0, Vector3::Arg p1, Vector3::Arg p2, Vector2 * z0, Vector2 * z1, Vector2 * z2) + { + Vector3 X = normalize(p1 - p0, 0.0f); + Vector3 Z = normalize(cross(X, (p2 - p0)), 0.0f); + Vector3 Y = normalize(cross(Z, X), 0.0f); + + float x0 = 0.0f; + float y0 = 0.0f; + float x1 = length(p1 - p0); + float y1 = 0.0f; + float x2 = dot((p2 - p0), X); + float y2 = dot((p2 - p0), Y); + + *z0 = Vector2(x0, y0); + *z1 = Vector2(x1, y1); + *z2 = Vector2(x2, y2); + } + + // LSCM equation, geometric form : + // (Z1 - Z0)(U2 - U0) = (Z2 - Z0)(U1 - U0) + // Where Uk = uk + i.vk is the complex number + // corresponding to (u,v) coords + // Zk = xk + i.yk is the complex number + // corresponding to local (x,y) coords + // cool: no divide with this expression, + // makes it more numerically stable in + // the presence of degenerate triangles. + + static void setup_conformal_map_relations(SparseMatrix & A, int row, const HalfEdge::Vertex * v0, const HalfEdge::Vertex * v1, const HalfEdge::Vertex * v2) + { + int id0 = v0->id; + int id1 = v1->id; + int id2 = v2->id; + + Vector3 p0 = v0->pos; + Vector3 p1 = v1->pos; + Vector3 p2 = v2->pos; + + Vector2 z0, z1, z2; + project_triangle(p0, p1, p2, &z0, &z1, &z2); + + Vector2 z01 = z1 - z0; + Vector2 z02 = z2 - z0; + + float a = z01.x; + float b = z01.y; + float c = z02.x; + float d = z02.y; + nvCheck(b == 0.0f); + + // Note : 2*id + 0 --> u + // 2*id + 1 --> v + int u0_id = 2 * id0 + 0; + int v0_id = 2 * id0 + 1; + int u1_id = 2 * id1 + 0; + int v1_id = 2 * id1 + 1; + int u2_id = 2 * id2 + 0; + int v2_id = 2 * id2 + 1; + + // Note : b = 0 + + // Real part + A.setCoefficient(u0_id, 2 * row + 0, -a+c); + A.setCoefficient(v0_id, 2 * row + 0, b-d); + A.setCoefficient(u1_id, 2 * row + 0, -c); + A.setCoefficient(v1_id, 2 * row + 0, d); + A.setCoefficient(u2_id, 2 * row + 0, a); + + // Imaginary part + A.setCoefficient(u0_id, 2 * row + 1, -b+d); + A.setCoefficient(v0_id, 2 * row + 1, -a+c); + A.setCoefficient(u1_id, 2 * row + 1, -d); + A.setCoefficient(v1_id, 2 * row + 1, -c); + A.setCoefficient(v2_id, 2 * row + 1, a); + } + + + // Conformal relations from Brecht Van Lommel (based on ABF): + + static float vec_angle_cos(Vector3::Arg v1, Vector3::Arg v2, Vector3::Arg v3) + { + Vector3 d1 = v1 - v2; + Vector3 d2 = v3 - v2; + return clamp(dot(d1, d2) / (length(d1) * length(d2)), -1.0f, 1.0f); + } + + static float vec_angle(Vector3::Arg v1, Vector3::Arg v2, Vector3::Arg v3) + { + float dot = vec_angle_cos(v1, v2, v3); + return acosf(dot); + } + + static void triangle_angles(Vector3::Arg v1, Vector3::Arg v2, Vector3::Arg v3, float *a1, float *a2, float *a3) + { + *a1 = vec_angle(v3, v1, v2); + *a2 = vec_angle(v1, v2, v3); + *a3 = PI - *a2 - *a1; + } + + static void triangle_cosines(Vector3::Arg v1, Vector3::Arg v2, Vector3::Arg v3, float *a1, float *a2, float *a3) + { + *a1 = vec_angle_cos(v3, v1, v2); + *a2 = vec_angle_cos(v1, v2, v3); + *a3 = vec_angle_cos(v2, v3, v1); + } + + static void setup_abf_relations(SparseMatrix & A, int row, const HalfEdge::Vertex * v0, const HalfEdge::Vertex * v1, const HalfEdge::Vertex * v2) + { + int id0 = v0->id; + int id1 = v1->id; + int id2 = v2->id; + + Vector3 p0 = v0->pos; + Vector3 p1 = v1->pos; + Vector3 p2 = v2->pos; + +#if 1 + // @@ IC: Wouldn't it be more accurate to return cos and compute 1-cos^2? + // It does indeed seem to be a little bit more robust. + // @@ Need to revisit this more carefully! + + float a0, a1, a2; + triangle_angles(p0, p1, p2, &a0, &a1, &a2); + + float s0 = sinf(a0); + float s1 = sinf(a1); + float s2 = sinf(a2); + + /*// Hack for degenerate triangles. + if (equal(s0, 0) && equal(s1, 0) && equal(s2, 0)) { + if (equal(a0, 0)) a0 += 0.001f; + if (equal(a1, 0)) a1 += 0.001f; + if (equal(a2, 0)) a2 += 0.001f; + + if (equal(a0, PI)) a0 = PI - a1 - a2; + if (equal(a1, PI)) a1 = PI - a0 - a2; + if (equal(a2, PI)) a2 = PI - a0 - a1; + + s0 = sinf(a0); + s1 = sinf(a1); + s2 = sinf(a2); + }*/ + + if (s1 > s0 && s1 > s2) + { + swap(s1, s2); + swap(s0, s1); + + swap(a1, a2); + swap(a0, a1); + + swap(id1, id2); + swap(id0, id1); + } + else if (s0 > s1 && s0 > s2) + { + swap(s0, s2); + swap(s0, s1); + + swap(a0, a2); + swap(a0, a1); + + swap(id0, id2); + swap(id0, id1); + } + + float c0 = cosf(a0); +#else + float c0, c1, c2; + triangle_cosines(p0, p1, p2, &c0, &c1, &c2); + + float s0 = 1 - c0*c0; + float s1 = 1 - c1*c1; + float s2 = 1 - c2*c2; + + nvDebugCheck(s0 != 0 || s1 != 0 || s2 != 0); + + if (s1 > s0 && s1 > s2) + { + swap(s1, s2); + swap(s0, s1); + + swap(c1, c2); + swap(c0, c1); + + swap(id1, id2); + swap(id0, id1); + } + else if (s0 > s1 && s0 > s2) + { + swap(s0, s2); + swap(s0, s1); + + swap(c0, c2); + swap(c0, c1); + + swap(id0, id2); + swap(id0, id1); + } +#endif + + float ratio = (s2 == 0.0f) ? 1.0f: s1/s2; + float cosine = c0 * ratio; + float sine = s0 * ratio; + + // Note : 2*id + 0 --> u + // 2*id + 1 --> v + int u0_id = 2 * id0 + 0; + int v0_id = 2 * id0 + 1; + int u1_id = 2 * id1 + 0; + int v1_id = 2 * id1 + 1; + int u2_id = 2 * id2 + 0; + int v2_id = 2 * id2 + 1; + + // Real part + A.setCoefficient(u0_id, 2 * row + 0, cosine - 1.0f); + A.setCoefficient(v0_id, 2 * row + 0, -sine); + A.setCoefficient(u1_id, 2 * row + 0, -cosine); + A.setCoefficient(v1_id, 2 * row + 0, sine); + A.setCoefficient(u2_id, 2 * row + 0, 1); + + // Imaginary part + A.setCoefficient(u0_id, 2 * row + 1, sine); + A.setCoefficient(v0_id, 2 * row + 1, cosine - 1.0f); + A.setCoefficient(u1_id, 2 * row + 1, -sine); + A.setCoefficient(v1_id, 2 * row + 1, -cosine); + A.setCoefficient(v2_id, 2 * row + 1, 1); + } + +} // namespace + + +bool nv::computeLeastSquaresConformalMap(HalfEdge::Mesh * mesh) +{ + nvDebugCheck(mesh != NULL); + + // For this to work properly, mesh should not have colocals that have the same + // attributes, unless you want the vertices to actually have different texcoords. + + const uint vertexCount = mesh->vertexCount(); + const uint D = 2 * vertexCount; + const uint N = 2 * countMeshTriangles(mesh); + + // N is the number of equations (one per triangle) + // D is the number of variables (one per vertex; there are 2 pinned vertices). + if (N < D - 4) { + return false; + } + + SparseMatrix A(D, N); + FullVector b(N); + FullVector x(D); + + // Fill b: + b.fill(0.0f); + + // Fill x: + HalfEdge::Vertex * v0; + HalfEdge::Vertex * v1; + if (!findApproximateDiameterVertices(mesh, &v0, &v1)) + { + // Mesh has no boundaries. + return false; + } + if (v0->tex == v1->tex) + { + // LSCM expects an existing parameterization. + return false; + } + + for (uint v = 0; v < vertexCount; v++) + { + HalfEdge::Vertex * vertex = mesh->vertexAt(v); + nvDebugCheck(vertex != NULL); + + // Initial solution. + x[2 * v + 0] = vertex->tex.x; + x[2 * v + 1] = vertex->tex.y; + } + + // Fill A: + const uint faceCount = mesh->faceCount(); + for (uint f = 0, t = 0; f < faceCount; f++) + { + const HalfEdge::Face * face = mesh->faceAt(f); + nvDebugCheck(face != NULL); + nvDebugCheck(face->edgeCount() == 3); + + const HalfEdge::Vertex * vertex0 = NULL; + + for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + const HalfEdge::Edge * edge = it.current(); + nvCheck(edge != NULL); + + if (vertex0 == NULL) + { + vertex0 = edge->vertex; + } + else if (edge->next->vertex != vertex0) + { + const HalfEdge::Vertex * vertex1 = edge->from(); + const HalfEdge::Vertex * vertex2 = edge->to(); + + setup_abf_relations(A, t, vertex0, vertex1, vertex2); + //setup_conformal_map_relations(A, t, vertex0, vertex1, vertex2); + + t++; + } + } + } + + const uint lockedParameters[] = + { + 2 * v0->id + 0, + 2 * v0->id + 1, + 2 * v1->id + 0, + 2 * v1->id + 1 + }; + + // Solve + LeastSquaresSolver(A, b, x, lockedParameters, 4, 0.000001f); + + // Map x back to texcoords: + for (uint v = 0; v < vertexCount; v++) + { + HalfEdge::Vertex * vertex = mesh->vertexAt(v); + nvDebugCheck(vertex != NULL); + + vertex->tex = Vector2(x[2 * v + 0], x[2 * v + 1]); + } + + return true; +} diff --git a/thirdparty/thekla_atlas/nvmesh/param/LeastSquaresConformalMap.h b/thirdparty/thekla_atlas/nvmesh/param/LeastSquaresConformalMap.h new file mode 100644 index 0000000000..51fbf193c8 --- /dev/null +++ b/thirdparty/thekla_atlas/nvmesh/param/LeastSquaresConformalMap.h @@ -0,0 +1,15 @@ +// Copyright NVIDIA Corporation 2008 -- Ignacio Castano <icastano@nvidia.com> + +#pragma once +#ifndef NV_MESH_LEASTSQUARESCONFORMALMAP_H +#define NV_MESH_LEASTSQUARESCONFORMALMAP_H + +namespace nv +{ + namespace HalfEdge { class Mesh; } + + bool computeLeastSquaresConformalMap(HalfEdge::Mesh * mesh); + +} // nv namespace + +#endif // NV_MESH_LEASTSQUARESCONFORMALMAP_H diff --git a/thirdparty/thekla_atlas/nvmesh/param/OrthogonalProjectionMap.cpp b/thirdparty/thekla_atlas/nvmesh/param/OrthogonalProjectionMap.cpp new file mode 100644 index 0000000000..d6e5e30561 --- /dev/null +++ b/thirdparty/thekla_atlas/nvmesh/param/OrthogonalProjectionMap.cpp @@ -0,0 +1,99 @@ +// This code is in the public domain -- castano@gmail.com + +#include "nvmesh.h" // pch + +#include "OrthogonalProjectionMap.h" + +#include "nvcore/Array.inl" + +#include "nvmath/Fitting.h" +#include "nvmath/Vector.inl" +#include "nvmath/Box.inl" +#include "nvmath/Plane.inl" + +#include "nvmesh/halfedge/Mesh.h" +#include "nvmesh/halfedge/Vertex.h" +#include "nvmesh/halfedge/Face.h" +#include "nvmesh/geometry/Bounds.h" + + +using namespace nv; + +bool nv::computeOrthogonalProjectionMap(HalfEdge::Mesh * mesh) +{ + Vector3 axis[2]; + +#if 1 + + uint vertexCount = mesh->vertexCount(); + Array<Vector3> points(vertexCount); + points.resize(vertexCount); + + for (uint i = 0; i < vertexCount; i++) + { + points[i] = mesh->vertexAt(i)->pos; + } + +#if 0 + axis[0] = Fit::computePrincipalComponent_EigenSolver(vertexCount, points.buffer()); + axis[0] = normalize(axis[0]); + + Plane plane = Fit::bestPlane(vertexCount, points.buffer()); + + Vector3 n = plane.vector(); + + axis[1] = cross(axis[0], n); + axis[1] = normalize(axis[1]); +#else + // Avoid redundant computations. + float matrix[6]; + Fit::computeCovariance(vertexCount, points.buffer(), matrix); + + if (matrix[0] == 0 && matrix[3] == 0 && matrix[5] == 0) { + return false; + } + + float eigenValues[3]; + Vector3 eigenVectors[3]; + if (!nv::Fit::eigenSolveSymmetric3(matrix, eigenValues, eigenVectors)) { + return false; + } + + axis[0] = normalize(eigenVectors[0]); + axis[1] = normalize(eigenVectors[1]); +#endif + + +#else + + // IC: I thought this was generally more robust, but turns out it's not even guaranteed to return a valid projection. Imagine a narrow quad perpendicular to one plane, but rotated so that the shortest axis of + // the bounding box is in the direction of that plane. + + // Use the shortest box axis + Box box = MeshBounds::box(mesh); + Vector3 dir = box.extents(); + + if (fabs(dir.x) <= fabs(dir.y) && fabs(dir.x) <= fabs(dir.z)) { + axis[0] = Vector3(0, 1, 0); + axis[1] = Vector3(0, 0, 1); + } + else if (fabs(dir.y) <= fabs(dir.z)) { + axis[0] = Vector3(1, 0, 0); + axis[1] = Vector3(0, 0, 1); + } + else { + axis[0] = Vector3(1, 0, 0); + axis[1] = Vector3(0, 1, 0); + } +#endif + + // Project vertices to plane. + for (HalfEdge::Mesh::VertexIterator it(mesh->vertices()); !it.isDone(); it.advance()) + { + HalfEdge::Vertex * vertex = it.current(); + vertex->tex.x = dot(axis[0], vertex->pos); + vertex->tex.y = dot(axis[1], vertex->pos); + } + + return true; +} diff --git a/thirdparty/thekla_atlas/nvmesh/param/OrthogonalProjectionMap.h b/thirdparty/thekla_atlas/nvmesh/param/OrthogonalProjectionMap.h new file mode 100644 index 0000000000..54920413d5 --- /dev/null +++ b/thirdparty/thekla_atlas/nvmesh/param/OrthogonalProjectionMap.h @@ -0,0 +1,15 @@ +// This code is in the public domain -- castano@gmail.com + +#pragma once +#ifndef NV_MESH_ORTHOGONALPROJECTIONMAP_H +#define NV_MESH_ORTHOGONALPROJECTIONMAP_H + +namespace nv +{ + namespace HalfEdge { class Mesh; } + + bool computeOrthogonalProjectionMap(HalfEdge::Mesh * mesh); + +} // nv namespace + +#endif // NV_MESH_ORTHOGONALPROJECTIONMAP_H diff --git a/thirdparty/thekla_atlas/nvmesh/param/ParameterizationQuality.cpp b/thirdparty/thekla_atlas/nvmesh/param/ParameterizationQuality.cpp new file mode 100644 index 0000000000..683ee603cd --- /dev/null +++ b/thirdparty/thekla_atlas/nvmesh/param/ParameterizationQuality.cpp @@ -0,0 +1,323 @@ +// Copyright NVIDIA Corporation 2008 -- Ignacio Castano <icastano@nvidia.com> + +#include "nvmesh.h" // pch + +#include "ParameterizationQuality.h" + +#include "nvmesh/halfedge/Mesh.h" +#include "nvmesh/halfedge/Face.h" +#include "nvmesh/halfedge/Vertex.h" +#include "nvmesh/halfedge/Edge.h" + +#include "nvmath/Vector.inl" + +#include "nvcore/Debug.h" + +#include <float.h> + + +using namespace nv; + +#if 0 +/* +float triangleConformalEnergy(Vector3 q[3], Vector2 p[3]) +{ +const Vector3 v1 = q[0]; +const Vector3 v2 = q[1]; +const Vector3 v3 = q[2]; + +const Vector2 w1 = p[0]; +const Vector2 w2 = p[1]; +const Vector2 w3 = p[2]; + +float x1 = v2.x() - v1.x(); +float x2 = v3.x() - v1.x(); +float y1 = v2.y() - v1.y(); +float y2 = v3.y() - v1.y(); +float z1 = v2.z() - v1.z(); +float z2 = v3.z() - v1.z(); + +float s1 = w2.x() - w1.x(); +float s2 = w3.x() - w1.x(); +float t1 = w2.y() - w1.y(); +float t2 = w3.y() - w1.y(); + +float r = 1.0f / (s1 * t2 - s2 * t1); +Vector3 sdir((t2 * x1 - t1 * x2) * r, (t2 * y1 - t1 * y2) * r, (t2 * z1 - t1 * z2) * r); +Vector3 tdir((s1 * x2 - s2 * x1) * r, (s1 * y2 - s2 * y1) * r, (s1 * z2 - s2 * z1) * r); + +Vector3 N = cross(v3-v1, v2-v1); + +// Rotate 90 around N. +} +*/ + +static float triangleConformalEnergy(Vector3 q[3], Vector2 p[3]) +{ + // Using Denis formulas: + Vector3 c0 = q[1] - q[2]; + Vector3 c1 = q[2] - q[0]; + Vector3 c2 = q[0] - q[1]; + + Vector3 N = cross(-c0, c1); + float T = length(N); // 2T + N = normalize(N, 0); + + float cot_alpha0 = dot(-c1, c2) / length(cross(-c1, c2)); + float cot_alpha1 = dot(-c2, c0) / length(cross(-c2, c0)); + float cot_alpha2 = dot(-c0, c1) / length(cross(-c0, c1)); + + Vector3 t0 = -cot_alpha1 * c1 + cot_alpha2 * c2; + Vector3 t1 = -cot_alpha2 * c2 + cot_alpha0 * c0; + Vector3 t2 = -cot_alpha0 * c0 + cot_alpha1 * c1; + + nvCheck(equal(length(t0), length(c0))); + nvCheck(equal(length(t1), length(c1))); + nvCheck(equal(length(t2), length(c2))); + nvCheck(equal(dot(t0, c0), 0)); + nvCheck(equal(dot(t1, c1), 0)); + nvCheck(equal(dot(t2, c2), 0)); + + // Gradients + Vector3 grad_u = 1.0f / T * (p[0].x * t0 + p[1].x * t1 + p[2].x * t2); + Vector3 grad_v = 1.0f / T * (p[0].y * t0 + p[1].y * t1 + p[2].y * t2); + + // Rotated gradients + Vector3 Jgrad_u = 1.0f / T * (p[0].x * c0 + p[1].x * c1 + p[2].x * c2); + Vector3 Jgrad_v = 1.0f / T * (p[0].y * c0 + p[1].y * c1 + p[2].y * c2); + + // Using Lengyel's formulas: + { + const Vector3 v1 = q[0]; + const Vector3 v2 = q[1]; + const Vector3 v3 = q[2]; + + const Vector2 w1 = p[0]; + const Vector2 w2 = p[1]; + const Vector2 w3 = p[2]; + + float x1 = v2.x - v1.x; + float x2 = v3.x - v1.x; + float y1 = v2.y - v1.y; + float y2 = v3.y - v1.y; + float z1 = v2.z - v1.z; + float z2 = v3.z - v1.z; + + float s1 = w2.x - w1.x; + float s2 = w3.x - w1.x; + float t1 = w2.y - w1.y; + float t2 = w3.y - w1.y; + + float r = 1.0f / (s1 * t2 - s2 * t1); + Vector3 sdir((t2 * x1 - t1 * x2) * r, (t2 * y1 - t1 * y2) * r, (t2 * z1 - t1 * z2) * r); + Vector3 tdir((s1 * x2 - s2 * x1) * r, (s1 * y2 - s2 * y1) * r, (s1 * z2 - s2 * z1) * r); + + Vector3 Jsdir = cross(N, sdir); + Vector3 Jtdir = cross(N, tdir); + + float x = 3; + } + + // check: sdir == grad_u + // check: tdir == grad_v + + return length(grad_u - Jgrad_v); +} +#endif // 0 + + +ParameterizationQuality::ParameterizationQuality() +{ + m_totalTriangleCount = 0; + m_flippedTriangleCount = 0; + m_zeroAreaTriangleCount = 0; + + m_parametricArea = 0.0f; + m_geometricArea = 0.0f; + + m_stretchMetric = 0.0f; + m_maxStretchMetric = 0.0f; + + m_conformalMetric = 0.0f; + m_authalicMetric = 0.0f; +} + +ParameterizationQuality::ParameterizationQuality(const HalfEdge::Mesh * mesh) +{ + nvDebugCheck(mesh != NULL); + + m_totalTriangleCount = 0; + m_flippedTriangleCount = 0; + m_zeroAreaTriangleCount = 0; + + m_parametricArea = 0.0f; + m_geometricArea = 0.0f; + + m_stretchMetric = 0.0f; + m_maxStretchMetric = 0.0f; + + m_conformalMetric = 0.0f; + m_authalicMetric = 0.0f; + + const uint faceCount = mesh->faceCount(); + for (uint f = 0; f < faceCount; f++) + { + const HalfEdge::Face * face = mesh->faceAt(f); + const HalfEdge::Vertex * vertex0 = NULL; + + Vector3 p[3]; + Vector2 t[3]; + + for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + const HalfEdge::Edge * edge = it.current(); + + if (vertex0 == NULL) + { + vertex0 = edge->vertex; + + p[0] = vertex0->pos; + t[0] = vertex0->tex; + } + else if (edge->to() != vertex0) + { + p[1] = edge->from()->pos; + p[2] = edge->to()->pos; + t[1] = edge->from()->tex; + t[2] = edge->to()->tex; + + processTriangle(p, t); + } + } + } + + if (m_flippedTriangleCount + m_zeroAreaTriangleCount == faceCount) + { + // If all triangles are flipped, then none is. + m_flippedTriangleCount = 0; + } + + nvDebugCheck(isFinite(m_parametricArea) && m_parametricArea >= 0); + nvDebugCheck(isFinite(m_geometricArea) && m_geometricArea >= 0); + nvDebugCheck(isFinite(m_stretchMetric)); + nvDebugCheck(isFinite(m_maxStretchMetric)); + nvDebugCheck(isFinite(m_conformalMetric)); + nvDebugCheck(isFinite(m_authalicMetric)); +} + +bool ParameterizationQuality::isValid() const +{ + return m_flippedTriangleCount == 0; // @@ Does not test for self-overlaps. +} + +float ParameterizationQuality::rmsStretchMetric() const +{ + if (m_geometricArea == 0) return 0.0f; + float normFactor = sqrtf(m_parametricArea / m_geometricArea); + return sqrtf(m_stretchMetric / m_geometricArea) * normFactor; +} + +float ParameterizationQuality::maxStretchMetric() const +{ + if (m_geometricArea == 0) return 0.0f; + float normFactor = sqrtf(m_parametricArea / m_geometricArea); + return m_maxStretchMetric * normFactor; +} + +float ParameterizationQuality::rmsConformalMetric() const +{ + if (m_geometricArea == 0) return 0.0f; + return sqrtf(m_conformalMetric / m_geometricArea); +} + +float ParameterizationQuality::maxAuthalicMetric() const +{ + if (m_geometricArea == 0) return 0.0f; + return sqrtf(m_authalicMetric / m_geometricArea); +} + +void ParameterizationQuality::operator += (const ParameterizationQuality & pq) +{ + m_totalTriangleCount += pq.m_totalTriangleCount; + m_flippedTriangleCount += pq.m_flippedTriangleCount; + m_zeroAreaTriangleCount += pq.m_zeroAreaTriangleCount; + + m_parametricArea += pq.m_parametricArea; + m_geometricArea += pq.m_geometricArea; + + m_stretchMetric += pq.m_stretchMetric; + m_maxStretchMetric = max(m_maxStretchMetric, pq.m_maxStretchMetric); + + m_conformalMetric += pq.m_conformalMetric; + m_authalicMetric += pq.m_authalicMetric; +} + + +void ParameterizationQuality::processTriangle(Vector3 q[3], Vector2 p[3]) +{ + m_totalTriangleCount++; + + // Evaluate texture stretch metric. See: + // - "Texture Mapping Progressive Meshes", Sander, Snyder, Gortler & Hoppe + // - "Mesh Parameterization: Theory and Practice", Siggraph'07 Course Notes, Hormann, Levy & Sheffer. + + float t1 = p[0].x; + float s1 = p[0].y; + float t2 = p[1].x; + float s2 = p[1].y; + float t3 = p[2].x; + float s3 = p[2].y; + + float geometricArea = length(cross(q[1] - q[0], q[2] - q[0])) / 2; + float parametricArea = ((s2 - s1)*(t3 - t1) - (s3 - s1)*(t2 - t1)) / 2; + + if (isZero(parametricArea)) + { + m_zeroAreaTriangleCount++; + return; + } + + Vector3 Ss = (q[0] * (t2- t3) + q[1] * (t3 - t1) + q[2] * (t1 - t2)) / (2 * parametricArea); + Vector3 St = (q[0] * (s3- s2) + q[1] * (s1 - s3) + q[2] * (s2 - s1)) / (2 * parametricArea); + + float a = dot(Ss, Ss); // E + float b = dot(Ss, St); // F + float c = dot(St, St); // G + + // Compute eigen-values of the first fundamental form: + float sigma1 = sqrtf(0.5f * max(0.0f, a + c - sqrtf(square(a - c) + 4 * square(b)))); // gamma uppercase, min eigenvalue. + float sigma2 = sqrtf(0.5f * max(0.0f, a + c + sqrtf(square(a - c) + 4 * square(b)))); // gamma lowercase, max eigenvalue. + nvCheck(sigma2 >= sigma1); + + // isometric: sigma1 = sigma2 = 1 + // conformal: sigma1 / sigma2 = 1 + // authalic: sigma1 * sigma2 = 1 + + float rmsStretch = sqrtf((a + c) * 0.5f); + float rmsStretch2 = sqrtf((square(sigma1) + square(sigma2)) * 0.5f); + nvDebugCheck(equal(rmsStretch, rmsStretch2, 0.01f)); + + if (parametricArea < 0.0f) + { + // Count flipped triangles. + m_flippedTriangleCount++; + + parametricArea = fabsf(parametricArea); + } + + m_stretchMetric += square(rmsStretch) * geometricArea; + m_maxStretchMetric = max(m_maxStretchMetric, sigma2); + + if (!isZero(sigma1, 0.000001f)) { + // sigma1 is zero when geometricArea is zero. + m_conformalMetric += (sigma2 / sigma1) * geometricArea; + } + m_authalicMetric += (sigma1 * sigma2) * geometricArea; + + // Accumulate total areas. + m_geometricArea += geometricArea; + m_parametricArea += parametricArea; + + + //triangleConformalEnergy(q, p); +} diff --git a/thirdparty/thekla_atlas/nvmesh/param/ParameterizationQuality.h b/thirdparty/thekla_atlas/nvmesh/param/ParameterizationQuality.h new file mode 100644 index 0000000000..342e26b889 --- /dev/null +++ b/thirdparty/thekla_atlas/nvmesh/param/ParameterizationQuality.h @@ -0,0 +1,56 @@ +// Copyright NVIDIA Corporation 2008 -- Ignacio Castano <icastano@nvidia.com> + +#pragma once +#ifndef NV_MESH_PARAMETERIZATIONQUALITY_H +#define NV_MESH_PARAMETERIZATIONQUALITY_H + +#include <nvmesh/nvmesh.h> + +namespace nv +{ + class Vector2; + class Vector3; + + namespace HalfEdge { class Mesh; } + + // Estimate quality of existing parameterization. + NVMESH_CLASS class ParameterizationQuality + { + public: + ParameterizationQuality(); + ParameterizationQuality(const HalfEdge::Mesh * mesh); + + bool isValid() const; + + float rmsStretchMetric() const; + float maxStretchMetric() const; + + float rmsConformalMetric() const; + float maxAuthalicMetric() const; + + void operator += (const ParameterizationQuality & pq); + + private: + + void processTriangle(Vector3 p[3], Vector2 t[3]); + + private: + + uint m_totalTriangleCount; + uint m_flippedTriangleCount; + uint m_zeroAreaTriangleCount; + + float m_parametricArea; + float m_geometricArea; + + float m_stretchMetric; + float m_maxStretchMetric; + + float m_conformalMetric; + float m_authalicMetric; + + }; + +} // nv namespace + +#endif // NV_MESH_PARAMETERIZATIONQUALITY_H diff --git a/thirdparty/thekla_atlas/nvmesh/param/SingleFaceMap.cpp b/thirdparty/thekla_atlas/nvmesh/param/SingleFaceMap.cpp new file mode 100644 index 0000000000..4b205de8bf --- /dev/null +++ b/thirdparty/thekla_atlas/nvmesh/param/SingleFaceMap.cpp @@ -0,0 +1,53 @@ +// Copyright NVIDIA Corporation 2008 -- Ignacio Castano <icastano@nvidia.com> + +#include "nvmesh.h" // pch + +#include "SingleFaceMap.h" + +#include "nvmesh/halfedge/Mesh.h" +#include "nvmesh/halfedge/Vertex.h" +#include "nvmesh/halfedge/Face.h" + +#include "nvmath/Vector.inl" + +using namespace nv; + + + +void nv::computeSingleFaceMap(HalfEdge::Mesh * mesh) +{ + nvDebugCheck(mesh != NULL); + nvDebugCheck(mesh->faceCount() == 1); + + HalfEdge::Face * face = mesh->faceAt(0); + nvCheck(face != NULL); + + Vector3 p0 = face->edge->from()->pos; + Vector3 p1 = face->edge->to()->pos; + + Vector3 X = normalizeSafe(p1 - p0, Vector3(0.0f), 0.0f); + Vector3 Z = face->normal(); + Vector3 Y = normalizeSafe(cross(Z, X), Vector3(0.0f), 0.0f); + + uint i = 0; + for (HalfEdge::Face::EdgeIterator it(face->edges()); !it.isDone(); it.advance(), i++) + { + HalfEdge::Vertex * vertex = it.vertex(); + nvCheck(vertex != NULL); + + if (i == 0) + { + vertex->tex = Vector2(0); + } + else + { + Vector3 pn = vertex->pos; + + float xn = dot((pn - p0), X); + float yn = dot((pn - p0), Y); + + vertex->tex = Vector2(xn, yn); + } + } +} + diff --git a/thirdparty/thekla_atlas/nvmesh/param/SingleFaceMap.h b/thirdparty/thekla_atlas/nvmesh/param/SingleFaceMap.h new file mode 100644 index 0000000000..b70719f5d8 --- /dev/null +++ b/thirdparty/thekla_atlas/nvmesh/param/SingleFaceMap.h @@ -0,0 +1,18 @@ +// Copyright NVIDIA Corporation 2008 -- Ignacio Castano <icastano@nvidia.com> + +#pragma once +#ifndef NV_MESH_SINGLEFACEMAP_H +#define NV_MESH_SINGLEFACEMAP_H + +namespace nv +{ + namespace HalfEdge + { + class Mesh; + } + + void computeSingleFaceMap(HalfEdge::Mesh * mesh); + +} // nv namespace + +#endif // NV_MESH_SINGLEFACEMAP_H diff --git a/thirdparty/thekla_atlas/nvmesh/param/Util.cpp b/thirdparty/thekla_atlas/nvmesh/param/Util.cpp new file mode 100644 index 0000000000..fe7b58edf8 --- /dev/null +++ b/thirdparty/thekla_atlas/nvmesh/param/Util.cpp @@ -0,0 +1,326 @@ +// This code is in the public domain -- castano@gmail.com + +#include "nvmesh.h" // pch + +#include "Util.h" + +#include "nvmesh/halfedge/Mesh.h" +#include "nvmesh/halfedge/Face.h" +#include "nvmesh/halfedge/Vertex.h" + +#include "nvmath/Vector.inl" + +#include "nvcore/Array.inl" + + +using namespace nv; + +// Determine if the given mesh is a quad mesh. +bool nv::isQuadMesh(const HalfEdge::Mesh * mesh) +{ + nvDebugCheck(mesh != NULL); + + const uint faceCount = mesh->faceCount(); + for(uint i = 0; i < faceCount; i++) { + const HalfEdge::Face * face = mesh->faceAt(i); + if (face->edgeCount() != 4) { + return false; + } + } + + return true; +} + +bool nv::isTriangularMesh(const HalfEdge::Mesh * mesh) +{ + for (HalfEdge::Mesh::ConstFaceIterator it(mesh->faces()); !it.isDone(); it.advance()) + { + const HalfEdge::Face * face = it.current(); + if (face->edgeCount() != 3) return false; + } + return true; +} + + +uint nv::countMeshTriangles(const HalfEdge::Mesh * mesh) +{ + const uint faceCount = mesh->faceCount(); + + uint triangleCount = 0; + + for (uint f = 0; f < faceCount; f++) + { + const HalfEdge::Face * face = mesh->faceAt(f); + + uint edgeCount = face->edgeCount(); + nvDebugCheck(edgeCount > 2); + + triangleCount += edgeCount - 2; + } + + return triangleCount; +} + +const HalfEdge::Vertex * nv::findBoundaryVertex(const HalfEdge::Mesh * mesh) +{ + const uint vertexCount = mesh->vertexCount(); + + for (uint v = 0; v < vertexCount; v++) + { + const HalfEdge::Vertex * vertex = mesh->vertexAt(v); + if (vertex->isBoundary()) return vertex; + } + + return NULL; +} + + +HalfEdge::Mesh * nv::unifyVertices(const HalfEdge::Mesh * inputMesh) +{ + HalfEdge::Mesh * mesh = new HalfEdge::Mesh; + + // Only add the first colocal. + const uint vertexCount = inputMesh->vertexCount(); + for (uint v = 0; v < vertexCount; v++) { + const HalfEdge::Vertex * vertex = inputMesh->vertexAt(v); + + if (vertex->isFirstColocal()) { + mesh->addVertex(vertex->pos); + } + } + + nv::Array<uint> indexArray; + + // Add new faces pointing to first colocals. + uint faceCount = inputMesh->faceCount(); + for (uint f = 0; f < faceCount; f++) { + const HalfEdge::Face * face = inputMesh->faceAt(f); + + indexArray.clear(); + + for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) { + const HalfEdge::Edge * edge = it.current(); + const HalfEdge::Vertex * vertex = edge->vertex->firstColocal(); + + indexArray.append(vertex->id); + } + + mesh->addFace(indexArray); + } + + mesh->linkBoundary(); + + return mesh; +} + +#include "nvmath/Basis.h" + +static bool pointInTriangle(const Vector2 & p, const Vector2 & a, const Vector2 & b, const Vector2 & c) +{ + return triangleArea(a, b, p) >= 0.00001f && + triangleArea(b, c, p) >= 0.00001f && + triangleArea(c, a, p) >= 0.00001f; +} + + +// This is doing a simple ear-clipping algorithm that skips invalid triangles. Ideally, we should +// also sort the ears by angle, start with the ones that have the smallest angle and proceed in order. +HalfEdge::Mesh * nv::triangulate(const HalfEdge::Mesh * inputMesh) +{ + HalfEdge::Mesh * mesh = new HalfEdge::Mesh; + + // Add all vertices. + const uint vertexCount = inputMesh->vertexCount(); + for (uint v = 0; v < vertexCount; v++) { + const HalfEdge::Vertex * vertex = inputMesh->vertexAt(v); + mesh->addVertex(vertex->pos); + } + + Array<int> polygonVertices; + Array<float> polygonAngles; + Array<Vector2> polygonPoints; + + const uint faceCount = inputMesh->faceCount(); + for (uint f = 0; f < faceCount; f++) + { + const HalfEdge::Face * face = inputMesh->faceAt(f); + nvDebugCheck(face != NULL); + + const uint edgeCount = face->edgeCount(); + nvDebugCheck(edgeCount >= 3); + + polygonVertices.clear(); + polygonVertices.reserve(edgeCount); + + if (edgeCount == 3) { + // Simple case for triangles. + for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + const HalfEdge::Edge * edge = it.current(); + const HalfEdge::Vertex * vertex = edge->vertex; + polygonVertices.append(vertex->id); + } + + int v0 = polygonVertices[0]; + int v1 = polygonVertices[1]; + int v2 = polygonVertices[2]; + + mesh->addFace(v0, v1, v2); + } + else { + // Build 2D polygon projecting vertices onto normal plane. + // Faces are not necesarily planar, this is for example the case, when the face comes from filling a hole. In such cases + // it's much better to use the best fit plane. + const Vector3 fn = face->normal(); + + Basis basis; + basis.buildFrameForDirection(fn); + + polygonPoints.clear(); + polygonPoints.reserve(edgeCount); + polygonAngles.clear(); + polygonAngles.reserve(edgeCount); + + for (HalfEdge::Face::ConstEdgeIterator it(face->edges()); !it.isDone(); it.advance()) + { + const HalfEdge::Edge * edge = it.current(); + const HalfEdge::Vertex * vertex = edge->vertex; + polygonVertices.append(vertex->id); + + Vector2 p; + p.x = dot(basis.tangent, vertex->pos); + p.y = dot(basis.bitangent, vertex->pos); + + polygonPoints.append(p); + } + polygonAngles.resize(edgeCount); + + while (polygonVertices.size() > 2) { + uint size = polygonVertices.size(); + + // Update polygon angles. @@ Update only those that have changed. + float minAngle = 2 * PI; + uint bestEar = 0; // Use first one if none of them is valid. + bool bestIsValid = false; + for (uint i = 0; i < size; i++) { + uint i0 = i; + uint i1 = (i+1) % size; // Use Sean's polygon interation trick. + uint i2 = (i+2) % size; + + Vector2 p0 = polygonPoints[i0]; + Vector2 p1 = polygonPoints[i1]; + Vector2 p2 = polygonPoints[i2]; + + float d = clamp(dot(p0-p1, p2-p1) / (length(p0-p1) * length(p2-p1)), -1.0f, 1.0f); + float angle = acosf(d); + + float area = triangleArea(p0, p1, p2); + if (area < 0.0f) angle = 2.0f * PI - angle; + + polygonAngles[i1] = angle; + + if (angle < minAngle || !bestIsValid) { + + // Make sure this is a valid ear, if not, skip this point. + bool valid = true; + for (uint j = 0; j < size; j++) { + if (j == i0 || j == i1 || j == i2) continue; + Vector2 p = polygonPoints[j]; + + if (pointInTriangle(p, p0, p1, p2)) { + valid = false; + break; + } + } + + if (valid || !bestIsValid) { + minAngle = angle; + bestEar = i1; + bestIsValid = valid; + } + } + } + + nvDebugCheck(minAngle <= 2 * PI); + + // Clip best ear: + + uint i0 = (bestEar+size-1) % size; + uint i1 = (bestEar+0) % size; + uint i2 = (bestEar+1) % size; + + int v0 = polygonVertices[i0]; + int v1 = polygonVertices[i1]; + int v2 = polygonVertices[i2]; + + mesh->addFace(v0, v1, v2); + + polygonVertices.removeAt(i1); + polygonPoints.removeAt(i1); + polygonAngles.removeAt(i1); + } + } + +#if 0 + + uint i = 0; + while (polygonVertices.size() > 2 && i < polygonVertices.size()) { + uint size = polygonVertices.size(); + uint i0 = (i+0) % size; + uint i1 = (i+1) % size; + uint i2 = (i+2) % size; + + const HalfEdge::Vertex * v0 = polygonVertices[i0]; + const HalfEdge::Vertex * v1 = polygonVertices[i1]; + const HalfEdge::Vertex * v2 = polygonVertices[i2]; + + const Vector3 p0 = v0->pos; + const Vector3 p1 = v1->pos; + const Vector3 p2 = v2->pos; + + const Vector3 e0 = p2 - p1; + const Vector3 e1 = p0 - p1; + + // If this ear forms a valid triangle, setup relations, remove v1 and repeat. + Vector3 n = cross(e0, e1); + float len = dot(fn, n); // = sin(angle) + + float angle = asin(len); + + + if (len > 0.0f) { + mesh->addFace(v0->id(), v1->id(), v2->id()); + polygonVertices.removeAt(i1); + polygonAngles.removeAt(i1); + if (i2 > i1) i2--; + // @@ Update angles at i0 and i2 + } + else { + i++; + } + } + + // @@ Create a few degenerate triangles to avoid introducing holes. + i = 0; + const uint size = polygonVertices.size(); + while (i < size - 2) { + uint i0 = (i+0) % size; + uint i1 = (i+1) % size; + uint i2 = (i+2) % size; + + const HalfEdge::Vertex * v0 = polygonVertices[i0]; + const HalfEdge::Vertex * v1 = polygonVertices[i1]; + const HalfEdge::Vertex * v2 = polygonVertices[i2]; + + mesh->addFace(v0->id(), v1->id(), v2->id()); + i++; + } +#endif + } + + mesh->linkBoundary(); + + return mesh; +} + + diff --git a/thirdparty/thekla_atlas/nvmesh/param/Util.h b/thirdparty/thekla_atlas/nvmesh/param/Util.h new file mode 100644 index 0000000000..774563ac0b --- /dev/null +++ b/thirdparty/thekla_atlas/nvmesh/param/Util.h @@ -0,0 +1,18 @@ +// This code is in the public domain -- castano@gmail.com + +#include "nvmesh/nvmesh.h" + +namespace nv { + + namespace HalfEdge { class Mesh; class Vertex; } + + bool isQuadMesh(const HalfEdge::Mesh * mesh); + bool isTriangularMesh(const HalfEdge::Mesh * mesh); + + uint countMeshTriangles(const HalfEdge::Mesh * mesh); + const HalfEdge::Vertex * findBoundaryVertex(const HalfEdge::Mesh * mesh); + + HalfEdge::Mesh * unifyVertices(const HalfEdge::Mesh * inputMesh); + HalfEdge::Mesh * triangulate(const HalfEdge::Mesh * inputMesh); + +} // nv namespace |