diff options
Diffstat (limited to 'thirdparty/thekla_atlas/nvmesh/param')
16 files changed, 0 insertions, 6008 deletions
diff --git a/thirdparty/thekla_atlas/nvmesh/param/Atlas.cpp b/thirdparty/thekla_atlas/nvmesh/param/Atlas.cpp deleted file mode 100644 index 98f92cef96..0000000000 --- a/thirdparty/thekla_atlas/nvmesh/param/Atlas.cpp +++ /dev/null @@ -1,1519 +0,0 @@ -// 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 deleted file mode 100644 index 41cfaea9cb..0000000000 --- a/thirdparty/thekla_atlas/nvmesh/param/Atlas.h +++ /dev/null @@ -1,186 +0,0 @@ -// 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 deleted file mode 100644 index bd2140c2f3..0000000000 --- a/thirdparty/thekla_atlas/nvmesh/param/AtlasBuilder.cpp +++ /dev/null @@ -1,1320 +0,0 @@ -// 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 deleted file mode 100644 index f25c724f7e..0000000000 --- a/thirdparty/thekla_atlas/nvmesh/param/AtlasBuilder.h +++ /dev/null @@ -1,111 +0,0 @@ -// 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 deleted file mode 100644 index 11e635db17..0000000000 --- a/thirdparty/thekla_atlas/nvmesh/param/AtlasPacker.cpp +++ /dev/null @@ -1,1401 +0,0 @@ -// 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; - - // -- GODOT start -- - //m_debug_bitmap.allocate(256, 256); - //m_debug_bitmap.fill(Color32(0,0,0,0)); - // -- GODOT end -- -} - -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); - // -- GODOT start -- - //nvDebugBreak(); - m_atlas->setFailed(); - return; - // -- GODOT end -- - } - //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); - // -- GODOT start -- - //m_debug_bitmap.resize(approximateExtent, approximateExtent); - //m_debug_bitmap.fill(Color32(0,0,0,0)); - // -- GODOT end -- - } - - - 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); - // -- GODOT start -- - //m_debug_bitmap.resize(nextPowerOfTwo(U32(w)), nextPowerOfTwo(U32(h))); - // -- GODOT end -- - } - - //nvDebug("Add chart at (%d, %d).\n", best_x, best_y); - - addChart(&chart_bitmap, w, h, best_x, best_y, best_r, /*debugOutput=*/NULL); - - // -- GODOT start -- - // 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); - }*/ - // -- GODOT end -- - - //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)); - - // -- GODOT start -- - //m_debug_bitmap.resize(m_width, m_height); - //m_debug_bitmap.setFormat(Image::Format_ARGB); - // -- GODOT end -- - -#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 deleted file mode 100644 index 845dbfb6f3..0000000000 --- a/thirdparty/thekla_atlas/nvmesh/param/AtlasPacker.h +++ /dev/null @@ -1,65 +0,0 @@ -// 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; - // -- GODOT start -- - //Image m_debug_bitmap; - // -- GODOT end -- - 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 deleted file mode 100644 index cd1e8bbb7b..0000000000 --- a/thirdparty/thekla_atlas/nvmesh/param/LeastSquaresConformalMap.cpp +++ /dev/null @@ -1,483 +0,0 @@ -// 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 deleted file mode 100644 index 51fbf193c8..0000000000 --- a/thirdparty/thekla_atlas/nvmesh/param/LeastSquaresConformalMap.h +++ /dev/null @@ -1,15 +0,0 @@ -// 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 deleted file mode 100644 index d6e5e30561..0000000000 --- a/thirdparty/thekla_atlas/nvmesh/param/OrthogonalProjectionMap.cpp +++ /dev/null @@ -1,99 +0,0 @@ -// 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 deleted file mode 100644 index 54920413d5..0000000000 --- a/thirdparty/thekla_atlas/nvmesh/param/OrthogonalProjectionMap.h +++ /dev/null @@ -1,15 +0,0 @@ -// 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 deleted file mode 100644 index 683ee603cd..0000000000 --- a/thirdparty/thekla_atlas/nvmesh/param/ParameterizationQuality.cpp +++ /dev/null @@ -1,323 +0,0 @@ -// 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 deleted file mode 100644 index 342e26b889..0000000000 --- a/thirdparty/thekla_atlas/nvmesh/param/ParameterizationQuality.h +++ /dev/null @@ -1,56 +0,0 @@ -// 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 deleted file mode 100644 index 4b205de8bf..0000000000 --- a/thirdparty/thekla_atlas/nvmesh/param/SingleFaceMap.cpp +++ /dev/null @@ -1,53 +0,0 @@ -// 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 deleted file mode 100644 index b70719f5d8..0000000000 --- a/thirdparty/thekla_atlas/nvmesh/param/SingleFaceMap.h +++ /dev/null @@ -1,18 +0,0 @@ -// 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 deleted file mode 100644 index fe7b58edf8..0000000000 --- a/thirdparty/thekla_atlas/nvmesh/param/Util.cpp +++ /dev/null @@ -1,326 +0,0 @@ -// 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 deleted file mode 100644 index 774563ac0b..0000000000 --- a/thirdparty/thekla_atlas/nvmesh/param/Util.h +++ /dev/null @@ -1,18 +0,0 @@ -// 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 |