From 59780fd046f4bcecbf892a47e47e2f22068b5d7c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Verschelde?= Date: Wed, 20 May 2020 13:03:57 +0200 Subject: xatlas: Sync with upstream 470576d --- COPYRIGHT.txt | 2 +- modules/xatlas_unwrap/register_types.cpp | 2 +- thirdparty/README.md | 2 +- thirdparty/xatlas/LICENSE | 2 +- thirdparty/xatlas/xatlas.cpp | 5225 +++++++++++++++++------------- thirdparty/xatlas/xatlas.h | 63 +- 6 files changed, 2964 insertions(+), 2332 deletions(-) diff --git a/COPYRIGHT.txt b/COPYRIGHT.txt index 7741573039..2d32bf1fd9 100644 --- a/COPYRIGHT.txt +++ b/COPYRIGHT.txt @@ -368,7 +368,7 @@ License: Expat Files: ./thirdparty/xatlas/ Comment: xatlas -Copyright: 2018, Jonathan Young +Copyright: 2018-2020, Jonathan Young 2013, Thekla, Inc 2006, NVIDIA Corporation, Ignacio Castano License: Expat diff --git a/modules/xatlas_unwrap/register_types.cpp b/modules/xatlas_unwrap/register_types.cpp index b996be2a4c..6242009f67 100644 --- a/modules/xatlas_unwrap/register_types.cpp +++ b/modules/xatlas_unwrap/register_types.cpp @@ -145,7 +145,7 @@ bool xatlas_mesh_lightmap_unwrap_callback(float p_texel_size, const float *p_ver ERR_FAIL_COND_V_MSG(err != xatlas::AddMeshError::Enum::Success, false, xatlas::StringForEnum(err)); printf("Generate..\n"); - xatlas::Generate(atlas, chart_options, nullptr, pack_options); + xatlas::Generate(atlas, chart_options, xatlas::ParameterizeOptions(), pack_options); *r_size_hint_x = atlas->width; *r_size_hint_y = atlas->height; diff --git a/thirdparty/README.md b/thirdparty/README.md index 5b24d2b96d..64c2ce336d 100644 --- a/thirdparty/README.md +++ b/thirdparty/README.md @@ -616,7 +616,7 @@ File extracted from upstream release tarball: ## xatlas - Upstream: https://github.com/jpcy/xatlas -- Version: git (e12ea82, 2019) +- Version: git (470576d3516f7e6d8b4554e7c941194a935969fd, 2020) - License: MIT Files extracted from upstream source: diff --git a/thirdparty/xatlas/LICENSE b/thirdparty/xatlas/LICENSE index 9e61e15fb7..7af8b33238 100644 --- a/thirdparty/xatlas/LICENSE +++ b/thirdparty/xatlas/LICENSE @@ -1,6 +1,6 @@ MIT License -Copyright (c) 2018-2019 Jonathan Young +Copyright (c) 2018-2020 Jonathan Young Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/thirdparty/xatlas/xatlas.cpp b/thirdparty/xatlas/xatlas.cpp index 80cacf9746..43aec33a9f 100644 --- a/thirdparty/xatlas/xatlas.cpp +++ b/thirdparty/xatlas/xatlas.cpp @@ -1,7 +1,7 @@ /* MIT License -Copyright (c) 2018-2019 Jonathan Young +Copyright (c) 2018-2020 Jonathan Young Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal @@ -90,8 +90,8 @@ Copyright (c) 2012 Brandon Pelfrey #endif #define XA_ALLOC(tag, type) (type *)internal::Realloc(nullptr, sizeof(type), tag, __FILE__, __LINE__) -#define XA_ALLOC_ARRAY(tag, type, num) (type *)internal::Realloc(nullptr, sizeof(type) * num, tag, __FILE__, __LINE__) -#define XA_REALLOC(tag, ptr, type, num) (type *)internal::Realloc(ptr, sizeof(type) * num, tag, __FILE__, __LINE__) +#define XA_ALLOC_ARRAY(tag, type, num) (type *)internal::Realloc(nullptr, sizeof(type) * (num), tag, __FILE__, __LINE__) +#define XA_REALLOC(tag, ptr, type, num) (type *)internal::Realloc(ptr, sizeof(type) * (num), tag, __FILE__, __LINE__) #define XA_REALLOC_SIZE(tag, ptr, size) (uint8_t *)internal::Realloc(ptr, size, tag, __FILE__, __LINE__) #define XA_FREE(ptr) internal::Realloc(ptr, 0, internal::MemTag::Default, __FILE__, __LINE__) #define XA_NEW(tag, type) new (XA_ALLOC(tag, type)) type() @@ -122,11 +122,12 @@ Copyright (c) 2012 Brandon Pelfrey #define XA_DEBUG_HEAP 0 #define XA_DEBUG_SINGLE_CHART 0 +#define XA_DEBUG_ALL_CHARTS_INVALID 0 #define XA_DEBUG_EXPORT_ATLAS_IMAGES 0 #define XA_DEBUG_EXPORT_ATLAS_IMAGES_PER_CHART 0 // Export an atlas image after each chart is added. #define XA_DEBUG_EXPORT_BOUNDARY_GRID 0 #define XA_DEBUG_EXPORT_TGA (XA_DEBUG_EXPORT_ATLAS_IMAGES || XA_DEBUG_EXPORT_BOUNDARY_GRID) -#define XA_DEBUG_EXPORT_OBJ_SOURCE_MESHES 0 +#define XA_DEBUG_EXPORT_OBJ_FACE_GROUPS 0 #define XA_DEBUG_EXPORT_OBJ_CHART_GROUPS 0 #define XA_DEBUG_EXPORT_OBJ_PLANAR_REGIONS 0 #define XA_DEBUG_EXPORT_OBJ_CHARTS 0 @@ -137,7 +138,7 @@ Copyright (c) 2012 Brandon Pelfrey #define XA_DEBUG_EXPORT_OBJ_RECOMPUTED_CHARTS 0 #define XA_DEBUG_EXPORT_OBJ (0 \ - || XA_DEBUG_EXPORT_OBJ_SOURCE_MESHES \ + || XA_DEBUG_EXPORT_OBJ_FACE_GROUPS \ || XA_DEBUG_EXPORT_OBJ_CHART_GROUPS \ || XA_DEBUG_EXPORT_OBJ_PLANAR_REGIONS \ || XA_DEBUG_EXPORT_OBJ_CHARTS \ @@ -163,6 +164,83 @@ static FreeFunc s_free = free; static PrintFunc s_print = printf; static bool s_printVerbose = false; +#if XA_PROFILE +#define XA_PROFILE_START(var) const clock_t var##Start = clock(); +#define XA_PROFILE_END(var) internal::s_profile.var += clock() - var##Start; +#define XA_PROFILE_PRINT_AND_RESET(label, var) XA_PRINT("%s%.2f seconds (%g ms)\n", label, internal::clockToSeconds(internal::s_profile.var), internal::clockToMs(internal::s_profile.var)); internal::s_profile.var = 0; +#define XA_PROFILE_ALLOC 0 + +struct ProfileData +{ +#if XA_PROFILE_ALLOC + std::atomic alloc; +#endif + clock_t addMeshReal; + clock_t addMeshCopyData; + std::atomic addMeshThread; + std::atomic addMeshCreateColocals; + clock_t computeChartsReal; + std::atomic computeChartsThread; + std::atomic createFaceGroups; + std::atomic extractInvalidMeshGeometry; + std::atomic chartGroupComputeChartsReal; + std::atomic chartGroupComputeChartsThread; + std::atomic createChartGroupMesh; + std::atomic createChartGroupMeshColocals; + std::atomic createChartGroupMeshBoundaries; + std::atomic buildAtlas; + std::atomic buildAtlasInit; + std::atomic planarCharts; + std::atomic clusteredCharts; + std::atomic clusteredChartsPlaceSeeds; + std::atomic clusteredChartsPlaceSeedsBoundaryIntersection; + std::atomic clusteredChartsRelocateSeeds; + std::atomic clusteredChartsReset; + std::atomic clusteredChartsGrow; + std::atomic clusteredChartsGrowBoundaryIntersection; + std::atomic clusteredChartsMerge; + std::atomic clusteredChartsFillHoles; + std::atomic copyChartFaces; + clock_t parameterizeChartsReal; + std::atomic parameterizeChartsThread; + std::atomic createChartMesh; + std::atomic fixChartMeshTJunctions; + std::atomic closeChartMeshHoles; + std::atomic parameterizeChartsOrthogonal; + std::atomic parameterizeChartsLSCM; + std::atomic parameterizeChartsRecompute; + std::atomic parameterizeChartsPiecewise; + std::atomic parameterizeChartsPiecewiseBoundaryIntersection; + std::atomic parameterizeChartsEvaluateQuality; + clock_t packCharts; + clock_t packChartsAddCharts; + std::atomic packChartsAddChartsThread; + std::atomic packChartsAddChartsRestoreTexcoords; + clock_t packChartsRasterize; + clock_t packChartsDilate; + clock_t packChartsFindLocation; + clock_t packChartsBlit; + clock_t buildOutputMeshes; +}; + +static ProfileData s_profile; + +static double clockToMs(clock_t c) +{ + return c * 1000.0 / CLOCKS_PER_SEC; +} + +static double clockToSeconds(clock_t c) +{ + return c / (double)CLOCKS_PER_SEC; +} +#else +#define XA_PROFILE_START(var) +#define XA_PROFILE_END(var) +#define XA_PROFILE_PRINT_AND_RESET(label, var) +#define XA_PROFILE_ALLOC 0 +#endif + struct MemTag { enum @@ -170,7 +248,6 @@ struct MemTag Default, BitImage, BVH, - FullVector, Matrix, Mesh, MeshBoundaries, @@ -180,6 +257,7 @@ struct MemTag MeshNormals, MeshPositions, MeshTexcoords, + OpenNL, SegmentAtlasChartCandidates, SegmentAtlasChartFaces, SegmentAtlasMeshData, @@ -305,7 +383,6 @@ static void PrintMemoryUsage() "Default", "BitImage", "BVH", - "FullVector", "Matrix", "Mesh", "MeshBoundaries", @@ -315,6 +392,7 @@ static void PrintMemoryUsage() "MeshNormals", "MeshPositions", "MeshTexcoords", + "OpenNL", "SegmentAtlasChartCandidates", "SegmentAtlasChartFaces", "SegmentAtlasMeshData", @@ -335,79 +413,21 @@ static void *Realloc(void *ptr, size_t size, int /*tag*/, const char * /*file*/, s_free(ptr); return nullptr; } +#if XA_PROFILE_ALLOC + XA_PROFILE_START(alloc) +#endif void *mem = s_realloc(ptr, size); - if (size > 0) { - XA_DEBUG_ASSERT(mem); - } +#if XA_PROFILE_ALLOC + XA_PROFILE_END(alloc) +#endif + XA_DEBUG_ASSERT(size <= 0 || (size > 0 && mem)); return mem; } #define XA_PRINT_MEM_USAGE #endif -#if XA_PROFILE -#define XA_PROFILE_START(var) const clock_t var##Start = clock(); -#define XA_PROFILE_END(var) internal::s_profile.var += clock() - var##Start; -#define XA_PROFILE_PRINT_AND_RESET(label, var) XA_PRINT("%s%.2f seconds (%g ms)\n", label, internal::clockToSeconds(internal::s_profile.var), internal::clockToMs(internal::s_profile.var)); internal::s_profile.var = 0; - -struct ProfileData -{ - clock_t addMeshReal; - clock_t addMeshCopyData; - std::atomic addMeshThread; - std::atomic addMeshCreateColocals; - std::atomic addMeshCreateFaceGroups; - std::atomic addMeshCreateChartGroupsReal; - std::atomic addMeshCreateChartGroupsThread; - clock_t computeChartsReal; - std::atomic computeChartsThread; - std::atomic buildAtlas; - std::atomic buildAtlasInit; - std::atomic buildAtlasPlaceSeeds; - std::atomic buildAtlasRelocateSeeds; - std::atomic buildAtlasResetCharts; - std::atomic buildAtlasGrowCharts; - std::atomic buildAtlasMergeCharts; - std::atomic buildAtlasFillHoles; - std::atomic createChartMeshesReal; - std::atomic createChartMeshesThread; - std::atomic fixChartMeshTJunctions; - std::atomic closeChartMeshHoles; - clock_t parameterizeChartsReal; - std::atomic parameterizeChartsThread; - std::atomic parameterizeChartsOrthogonal; - std::atomic parameterizeChartsLSCM; - std::atomic parameterizeChartsEvaluateQuality; - clock_t packCharts; - clock_t packChartsAddCharts; - std::atomic packChartsAddChartsThread; - std::atomic packChartsAddChartsRestoreTexcoords; - clock_t packChartsRasterize; - clock_t packChartsDilate; - clock_t packChartsFindLocation; - clock_t packChartsBlit; - clock_t buildOutputMeshes; -}; - -static ProfileData s_profile; - -static double clockToMs(clock_t c) -{ - return c * 1000.0 / CLOCKS_PER_SEC; -} - -static double clockToSeconds(clock_t c) -{ - return c / (double)CLOCKS_PER_SEC; -} -#else -#define XA_PROFILE_START(var) -#define XA_PROFILE_END(var) -#define XA_PROFILE_PRINT_AND_RESET(label, var) -#endif - static constexpr float kPi = 3.14159265358979323846f; static constexpr float kPi2 = 6.28318530717958647692f; -static constexpr float kPi4 = 12.56637061435917295384f; static constexpr float kEpsilon = 0.0001f; static constexpr float kAreaEpsilon = FLT_EPSILON; static constexpr float kNormalEpsilon = 0.001f; @@ -517,33 +537,6 @@ static uint32_t nextPowerOfTwo(uint32_t x) return x + 1; } -static uint32_t sdbmHash(const void *data_in, uint32_t size, uint32_t h = 5381) -{ - const uint8_t *data = (const uint8_t *) data_in; - uint32_t i = 0; - while (i < size) { - h = (h << 16) + (h << 6) - h + (uint32_t ) data[i++]; - } - return h; -} - -template -static uint32_t hash(const T &t, uint32_t h = 5381) -{ - return sdbmHash(&t, sizeof(T), h); -} - -// Functors for hash table: -template struct Hash -{ - uint32_t operator()(const Key &k) const { return hash(k); } -}; - -template struct Equal -{ - bool operator()(const Key &k0, const Key &k1) const { return k0 == k1; } -}; - class Vector2 { public: @@ -860,6 +853,14 @@ struct Extents2 { Vector2 min, max; + Extents2() {} + + Extents2(Vector2 p1, Vector2 p2) + { + min = xatlas::internal::min(p1, p2); + max = xatlas::internal::max(p1, p2); + } + void reset() { min.x = min.y = FLT_MAX; @@ -877,7 +878,7 @@ struct Extents2 return Vector2(min.x + (max.x - min.x) * 0.5f, min.y + (max.y - min.y) * 0.5f); } - static bool intersect(Extents2 e1, Extents2 e2) + static bool intersect(const Extents2 &e1, const Extents2 &e2) { return e1.min.x <= e2.max.x && e1.max.x >= e2.min.x && e1.min.y <= e2.max.y && e1.max.y >= e2.min.y; } @@ -1123,6 +1124,15 @@ struct ArrayBase size--; } + // Element at index is swapped with the last element, then the array length is decremented. + void removeAtFast(uint32_t index) + { + XA_DEBUG_ASSERT(index >= 0 && index < size); + if (size != 1 && index != size - 1) + memcpy(buffer + elementSize * index, buffer + elementSize * (size - 1), elementSize); + size--; + } + void reserve(uint32_t desiredSize) { if (desiredSize > capacity) @@ -1164,9 +1174,9 @@ struct ArrayBase } #if XA_DEBUG_HEAP - void setMemTag(int memTag) + void setMemTag(int _memTag) { - this->memTag = memTag; + this->memTag = _memTag; } #endif @@ -1221,6 +1231,7 @@ public: void copyTo(Array &other) const { m_base.copyTo(other.m_base); } XA_INLINE const T *data() const { return (const T *)m_base.buffer; } XA_INLINE T *data() { return (T *)m_base.buffer; } + void destroy() { m_base.destroy(); } XA_INLINE T *end() { return (T *)m_base.buffer + m_base.size; } XA_INLINE bool isEmpty() const { return m_base.size == 0; } void insertAt(uint32_t index, const T &value) { m_base.insertAt(index, (const uint8_t *)&value); } @@ -1229,6 +1240,7 @@ public: void push_back(const Array &other) { m_base.push_back(other.m_base); } void pop_back() { m_base.pop_back(); } void removeAt(uint32_t index) { m_base.removeAt(index); } + void removeAtFast(uint32_t index) { m_base.removeAtFast(index); } void reserve(uint32_t desiredSize) { m_base.reserve(desiredSize); } void resize(uint32_t newSize) { m_base.resize(newSize, true); } @@ -1244,13 +1256,18 @@ public: ((T *)m_base.buffer)[i].~T(); } - void setAll(const T &value) + void fill(const T &value) { auto buffer = (T *)m_base.buffer; for (uint32_t i = 0; i < m_base.size; i++) buffer[i] = value; } + void fillBytes(uint8_t value) + { + memset(m_base.buffer, (int)value, m_base.size * m_base.elementSize); + } + #if XA_DEBUG_HEAP void setMemTag(int memTag) { m_base.setMemTag(memTag); } #endif @@ -1265,6 +1282,7 @@ private: template struct ArrayView { + ArrayView() : data(nullptr), length(0) {} ArrayView(Array &a) : data(a.data()), length(a.size()) {} ArrayView(T *data, uint32_t length) : data(data), length(length) {} ArrayView &operator=(Array &a) { data = a.data(); length = a.size(); return *this; } @@ -1276,6 +1294,7 @@ struct ArrayView template struct ConstArrayView { + ConstArrayView() : data(nullptr), length(0) {} ConstArrayView(const Array &a) : data(a.data()), length(a.size()) {} ConstArrayView(const T *data, uint32_t length) : data(data), length(length) {} ConstArrayView &operator=(const Array &a) { data = a.data(); length = a.size(); return *this; } @@ -1340,7 +1359,13 @@ public: void set(uint32_t index) { XA_DEBUG_ASSERT(index < m_size); - m_wordArray[index >> 5] |= (1 << (index & 31)); + m_wordArray[index >> 5] |= (1 << (index & 31)); + } + + void unset(uint32_t index) + { + XA_DEBUG_ASSERT(index < m_size); + m_wordArray[index >> 5] &= ~(1 << (index & 31)); } void zeroOutMemory() @@ -1927,26 +1952,38 @@ private: } }; -/// Fixed size vector class. -class FullVector +static uint32_t sdbmHash(const void *data_in, uint32_t size, uint32_t h = 5381) { -public: - FullVector(uint32_t dim) : m_array(MemTag::FullVector) { m_array.resize(dim); } - FullVector(const FullVector &v) : m_array(MemTag::FullVector) { v.m_array.copyTo(m_array); } - FullVector &operator=(const FullVector &v) = delete; - XA_INLINE uint32_t dimension() const { return m_array.size(); } - XA_INLINE const float &operator[](uint32_t index) const { return m_array[index]; } - XA_INLINE float &operator[](uint32_t index) { return m_array[index]; } - - void fill(float f) - { - const uint32_t dim = dimension(); - for (uint32_t i = 0; i < dim; i++) - m_array[i] = f; + const uint8_t *data = (const uint8_t *) data_in; + uint32_t i = 0; + while (i < size) { + h = (h << 16) + (h << 6) - h + (uint32_t ) data[i++]; } + return h; +} -private: - Array m_array; +template +static uint32_t hash(const T &t, uint32_t h = 5381) +{ + return sdbmHash(&t, sizeof(T), h); +} + +template +struct Hash +{ + uint32_t operator()(const Key &k) const { return hash(k); } +}; + +template +struct PassthroughHash +{ + uint32_t operator()(const Key &k) const { return (uint32_t)k; } +}; + +template +struct Equal +{ + bool operator()(const Key &k0, const Key &k1) const { return k0 == k1; } }; template, typename E = Equal > @@ -1963,7 +2000,17 @@ public: XA_FREE(m_slots); } - void add(const Key &key) + void destroy() + { + if (m_slots) { + XA_FREE(m_slots); + m_slots = nullptr; + } + m_keys.destroy(); + m_next.destroy(); + } + + uint32_t add(const Key &key) { if (!m_slots) alloc(); @@ -1971,6 +2018,7 @@ public: m_keys.push_back(key); m_next.push_back(m_slots[hash]); m_slots[hash] = m_next.size() - 1; + return m_keys.size() - 1; } uint32_t get(const Key &key) const @@ -2066,7 +2114,7 @@ public: y ^= (y << 5); uint64_t t = 698769069ULL * z + c; c = (t >> 32); - return (x + y + (z = (uint32_t)t)) % range; + return (x + y + (z = (uint32_t)t)) % (range + 1); } private: @@ -2079,45 +2127,38 @@ private: class RadixSort { public: - RadixSort() : m_size(0), m_ranks(nullptr), m_ranks2(nullptr), m_validRanks(false) {} - - ~RadixSort() - { - // Release everything - XA_FREE(m_ranks2); - XA_FREE(m_ranks); - } - - RadixSort &sort(const float *input, uint32_t count) + void sort(const float *input, uint32_t count) { - if (input == nullptr || count == 0) return *this; - // Resize lists if needed - if (count != m_size) { - if (count > m_size) { - m_ranks2 = XA_REALLOC(MemTag::Default, m_ranks2, uint32_t, count); - m_ranks = XA_REALLOC(MemTag::Default, m_ranks, uint32_t, count); - } - m_size = count; - m_validRanks = false; + if (input == nullptr || count == 0) { + m_buffer1.clear(); + m_buffer2.clear(); + m_ranks = m_buffer1.data(); + m_ranks2 = m_buffer2.data(); + return; } - if (count < 32) { + // Resize lists if needed + m_buffer1.resize(count); + m_buffer2.resize(count); + m_ranks = m_buffer1.data(); + m_ranks2 = m_buffer2.data(); + m_validRanks = false; + if (count < 32) insertionSort(input, count); - } else { + else { // @@ Avoid touching the input multiple times. for (uint32_t i = 0; i < count; i++) { - FloatFlip((uint32_t &)input[i]); + floatFlip((uint32_t &)input[i]); } radixSort((const uint32_t *)input, count); for (uint32_t i = 0; i < count; i++) { - IFloatFlip((uint32_t &)input[i]); + ifloatFlip((uint32_t &)input[i]); } } - return *this; } - RadixSort &sort(const Array &input) + void sort(const Array &input) { - return sort(input.data(), input.size()); + sort(input.data(), input.size()); } // Access to results. m_ranks is a list of indices in sorted order, i.e. in the order you may further process your data @@ -2127,25 +2168,18 @@ public: return m_ranks; } - uint32_t *ranks() - { - XA_DEBUG_ASSERT(m_validRanks); - return m_ranks; - } - private: - uint32_t m_size; - uint32_t *m_ranks; - uint32_t *m_ranks2; + uint32_t *m_ranks, *m_ranks2; + Array m_buffer1, m_buffer2; bool m_validRanks; - void FloatFlip(uint32_t &f) + void floatFlip(uint32_t &f) { int32_t mask = (int32_t(f) >> 31) | 0x80000000; // Warren Hunt, Manchor Ko. f ^= mask; } - void IFloatFlip(uint32_t &f) + void ifloatFlip(uint32_t &f) { uint32_t mask = ((f >> 31) - 1) | 0x80000000; // Michael Herf. f ^= mask; @@ -2179,7 +2213,8 @@ private: } } - template void insertionSort(const T *input, uint32_t count) + template + void insertionSort(const T *input, uint32_t count) { if (!m_validRanks) { m_ranks[0] = 0; @@ -2210,7 +2245,8 @@ private: } } - template void radixSort(const T *input, uint32_t count) + template + void radixSort(const T *input, uint32_t count) { const uint32_t P = sizeof(T); // pass count // Allocate histograms & offsets on the stack @@ -2247,9 +2283,8 @@ private: } // All values were equal, generate linear ranks. if (!m_validRanks) { - for (uint32_t i = 0; i < count; i++) { + for (uint32_t i = 0; i < count; i++) m_ranks[i] = i; - } m_validRanks = true; } } @@ -2327,9 +2362,8 @@ private: m_coords.resize(inputCount); for (uint32_t i = 0; i < inputCount; i++) m_coords[i] = input[i].x; - RadixSort radix; - radix.sort(m_coords); - const uint32_t *ranks = radix.ranks(); + m_radix.sort(m_coords); + const uint32_t *ranks = m_radix.ranks(); m_top.clear(); m_bottom.clear(); m_top.reserve(inputCount); @@ -2389,6 +2423,7 @@ private: Array m_boundaryVertices; Array m_coords; Array m_top, m_bottom, m_hull; + RadixSort m_radix; }; static uint32_t meshEdgeFace(uint32_t edge) { return edge / 3; } @@ -2404,9 +2439,8 @@ struct MeshFlags { enum { - HasFaceGroups = 1<<0, - HasIgnoredFaces = 1<<1, - HasNormals = 1<<2 + HasIgnoredFaces = 1<<0, + HasNormals = 1<<1 }; }; @@ -2416,20 +2450,17 @@ static void meshGetBoundaryLoops(const Mesh &mesh, Array &boundaryLoop class Mesh { public: - Mesh(float epsilon, uint32_t approxVertexCount, uint32_t approxFaceCount, uint32_t flags = 0, uint32_t id = UINT32_MAX) : m_epsilon(epsilon), m_flags(flags), m_id(id), m_faceIgnore(MemTag::Mesh), m_ignoredFaceCount(0), m_indices(MemTag::MeshIndices), m_positions(MemTag::MeshPositions), m_normals(MemTag::MeshNormals), m_texcoords(MemTag::MeshTexcoords), m_faceGroups(MemTag::Mesh), m_faceGroupFirstFace(MemTag::Mesh), m_faceGroupNextFace(MemTag::Mesh), m_faceGroupFaceCounts(MemTag::Mesh), m_colocalVertexCount(0), m_nextColocalVertex(MemTag::MeshColocals), m_boundaryEdges(MemTag::MeshBoundaries), m_oppositeEdges(MemTag::MeshBoundaries), m_nextBoundaryEdges(MemTag::MeshBoundaries), m_edgeMap(MemTag::MeshEdgeMap, approxFaceCount * 3) + Mesh(float epsilon, uint32_t approxVertexCount, uint32_t approxFaceCount, uint32_t flags = 0, uint32_t id = UINT32_MAX) : m_epsilon(epsilon), m_flags(flags), m_id(id), m_faceIgnore(MemTag::Mesh), m_indices(MemTag::MeshIndices), m_positions(MemTag::MeshPositions), m_normals(MemTag::MeshNormals), m_texcoords(MemTag::MeshTexcoords), m_nextColocalVertex(MemTag::MeshColocals), m_boundaryEdges(MemTag::MeshBoundaries), m_oppositeEdges(MemTag::MeshBoundaries), m_nextBoundaryEdges(MemTag::MeshBoundaries), m_edgeMap(MemTag::MeshEdgeMap, approxFaceCount * 3) { m_indices.reserve(approxFaceCount * 3); m_positions.reserve(approxVertexCount); m_texcoords.reserve(approxVertexCount); - if (m_flags & MeshFlags::HasFaceGroups) - m_faceGroups.reserve(approxFaceCount); if (m_flags & MeshFlags::HasIgnoredFaces) m_faceIgnore.reserve(approxFaceCount); if (m_flags & MeshFlags::HasNormals) m_normals.reserve(approxVertexCount); } - static constexpr uint16_t kInvalidFaceGroup = UINT16_MAX; uint32_t flags() const { return m_flags; } uint32_t id() const { return m_id; } @@ -2451,37 +2482,30 @@ public: }; }; - AddFaceResult::Enum addFace(uint32_t v0, uint32_t v1, uint32_t v2, bool ignore = false, bool hashEdge = true) + AddFaceResult::Enum addFace(uint32_t v0, uint32_t v1, uint32_t v2, bool ignore = false) { uint32_t indexArray[3]; indexArray[0] = v0; indexArray[1] = v1; indexArray[2] = v2; - return addFace(indexArray, ignore, hashEdge); + return addFace(indexArray, ignore); } - AddFaceResult::Enum addFace(const uint32_t *indices, bool ignore = false, bool hashEdge = true) + AddFaceResult::Enum addFace(const uint32_t *indices, bool ignore = false) { AddFaceResult::Enum result = AddFaceResult::OK; - if (m_flags & MeshFlags::HasFaceGroups) - m_faceGroups.push_back(kInvalidFaceGroup); - if (m_flags & MeshFlags::HasIgnoredFaces) { + if (m_flags & MeshFlags::HasIgnoredFaces) m_faceIgnore.push_back(ignore); - if (ignore) - m_ignoredFaceCount++; - } const uint32_t firstIndex = m_indices.size(); for (uint32_t i = 0; i < 3; i++) m_indices.push_back(indices[i]); - if (hashEdge) { - for (uint32_t i = 0; i < 3; i++) { - const uint32_t vertex0 = m_indices[firstIndex + i]; - const uint32_t vertex1 = m_indices[firstIndex + (i + 1) % 3]; - const EdgeKey key(vertex0, vertex1); - if (m_edgeMap.get(key) != UINT32_MAX) - result = AddFaceResult::DuplicateEdge; - m_edgeMap.add(key); - } + for (uint32_t i = 0; i < 3; i++) { + const uint32_t vertex0 = m_indices[firstIndex + i]; + const uint32_t vertex1 = m_indices[firstIndex + (i + 1) % 3]; + const EdgeKey key(vertex0, vertex1); + if (m_edgeMap.get(key) != UINT32_MAX) + result = AddFaceResult::DuplicateEdge; + m_edgeMap.add(key); } return result; } @@ -2496,10 +2520,8 @@ public: BVH bvh(aabbs); Array colocals(MemTag::MeshColocals); Array potential(MemTag::MeshColocals); - m_colocalVertexCount = 0; m_nextColocalVertex.resize(vertexCount); - for (uint32_t i = 0; i < vertexCount; i++) - m_nextColocalVertex[i] = UINT32_MAX; + m_nextColocalVertex.fillBytes(0xff); for (uint32_t i = 0; i < vertexCount; i++) { if (m_nextColocalVertex[i] != UINT32_MAX) continue; // Already linked. @@ -2517,7 +2539,6 @@ public: m_nextColocalVertex[i] = i; continue; } - m_colocalVertexCount += colocals.size(); // Link in ascending order. insertionSort(colocals.data(), colocals.size()); for (uint32_t j = 0; j < colocals.size(); j++) @@ -2526,150 +2547,57 @@ public: } } - // Check if the face duplicates any edges of any face already in the group. - bool faceDuplicatesGroupEdge(uint16_t group, uint32_t face) const + void createBoundaries() { - for (FaceEdgeIterator edgeIt(this, face); !edgeIt.isDone(); edgeIt.advance()) { - for (ColocalEdgeIterator colocalEdgeIt(this, edgeIt.vertex0(), edgeIt.vertex1()); !colocalEdgeIt.isDone(); colocalEdgeIt.advance()) { - if (m_faceGroups[meshEdgeFace(colocalEdgeIt.edge())] == group) - return true; + const uint32_t edgeCount = m_indices.size(); + const uint32_t vertexCount = m_positions.size(); + m_oppositeEdges.resize(edgeCount); + m_boundaryEdges.reserve(uint32_t(edgeCount * 0.1f)); + m_isBoundaryVertex.resize(vertexCount); + m_isBoundaryVertex.zeroOutMemory(); + for (uint32_t i = 0; i < edgeCount; i++) + m_oppositeEdges[i] = UINT32_MAX; + const uint32_t faceCount = m_indices.size() / 3; + for (uint32_t i = 0; i < faceCount; i++) { + if (isFaceIgnored(i)) + continue; + for (uint32_t j = 0; j < 3; j++) { + const uint32_t edge = i * 3 + j; + const uint32_t vertex0 = m_indices[edge]; + const uint32_t vertex1 = m_indices[i * 3 + (j + 1) % 3]; + // If there is an edge with opposite winding to this one, the edge isn't on a boundary. + const uint32_t oppositeEdge = findEdge(vertex1, vertex0); + if (oppositeEdge != UINT32_MAX) { + m_oppositeEdges[edge] = oppositeEdge; + } else { + m_boundaryEdges.push_back(edge); + m_isBoundaryVertex.set(vertex0); + m_isBoundaryVertex.set(vertex1); + } } } - return false; } - void createFaceGroups() + void linkBoundaries() { - uint32_t firstUnassignedFace = 0; - uint16_t group = 0; - Array growFaces; - const uint32_t n = faceCount(); - m_faceGroupNextFace.resize(n); + const uint32_t edgeCount = m_indices.size(); + HashMap vertexToEdgeMap(MemTag::Mesh, edgeCount); // Edge is index / 2 + for (uint32_t i = 0; i < edgeCount; i++) { + vertexToEdgeMap.add(m_indices[meshEdgeIndex0(i)]); + vertexToEdgeMap.add(m_indices[meshEdgeIndex1(i)]); + } + m_nextBoundaryEdges.resize(edgeCount); + for (uint32_t i = 0; i < edgeCount; i++) + m_nextBoundaryEdges[i] = UINT32_MAX; + uint32_t numBoundaryLoops = 0, numUnclosedBoundaries = 0; + BitArray linkedEdges(edgeCount); + linkedEdges.zeroOutMemory(); for (;;) { - // Find an unassigned face. - uint32_t face = UINT32_MAX; - for (uint32_t f = firstUnassignedFace; f < n; f++) { - if (m_faceGroups[f] == kInvalidFaceGroup && !isFaceIgnored(f)) { - face = f; - firstUnassignedFace = f + 1; - break; - } - } - if (face == UINT32_MAX) - break; // All faces assigned to a group (except ignored faces). - m_faceGroups[face] = group; - m_faceGroupNextFace[face] = UINT32_MAX; - m_faceGroupFirstFace.push_back(face); - growFaces.clear(); - growFaces.push_back(face); - uint32_t prevFace = face, groupFaceCount = 1; - // Find faces connected to the face and assign them to the same group as the face, unless they are already assigned to another group. - for (;;) { - if (growFaces.isEmpty()) - break; - const uint32_t f = growFaces.back(); - growFaces.pop_back(); - for (FaceEdgeIterator edgeIt(this, f); !edgeIt.isDone(); edgeIt.advance()) { - // Iterate opposite edges. There may be more than one - non-manifold geometry can have duplicate edges. - // Prioritize the one with exact vertex match, not just colocal. - // If *any* of the opposite edges are already assigned to this group, don't do anything. - bool alreadyAssignedToThisGroup = false; - uint32_t bestConnectedFace = UINT32_MAX; - for (ColocalEdgeIterator oppositeEdgeIt(this, edgeIt.vertex1(), edgeIt.vertex0()); !oppositeEdgeIt.isDone(); oppositeEdgeIt.advance()) { - const uint32_t oppositeEdge = oppositeEdgeIt.edge(); - const uint32_t oppositeFace = meshEdgeFace(oppositeEdge); - if (isFaceIgnored(oppositeFace)) - continue; // Don't add ignored faces to group. - if (m_faceGroups[oppositeFace] == group) { - alreadyAssignedToThisGroup = true; - break; - } - if (m_faceGroups[oppositeFace] != kInvalidFaceGroup) - continue; // Connected face is already assigned to another group. - if (faceDuplicatesGroupEdge(group, oppositeFace)) - continue; // Don't want duplicate edges in a group. - const uint32_t oppositeVertex0 = m_indices[meshEdgeIndex0(oppositeEdge)]; - const uint32_t oppositeVertex1 = m_indices[meshEdgeIndex1(oppositeEdge)]; - if (bestConnectedFace == UINT32_MAX || (oppositeVertex0 == edgeIt.vertex1() && oppositeVertex1 == edgeIt.vertex0())) - bestConnectedFace = oppositeFace; -#if 0 - else { - // Choose the opposite face with the smallest dihedral angle. - const float d1 = 1.0f - dot(computeFaceNormal(f), computeFaceNormal(bestConnectedFace)); - const float d2 = 1.0f - dot(computeFaceNormal(f), computeFaceNormal(oppositeFace)); - if (d2 < d1) - bestConnectedFace = oppositeFace; - } -#endif - } - if (!alreadyAssignedToThisGroup && bestConnectedFace != UINT32_MAX) { - m_faceGroups[bestConnectedFace] = group; - m_faceGroupNextFace[bestConnectedFace] = UINT32_MAX; - if (prevFace != UINT32_MAX) - m_faceGroupNextFace[prevFace] = bestConnectedFace; - prevFace = bestConnectedFace; - groupFaceCount++; - growFaces.push_back(bestConnectedFace); - } - } - } - m_faceGroupFaceCounts.push_back(groupFaceCount); - group++; - XA_ASSERT(group < kInvalidFaceGroup); - } - } - - void createBoundaries() - { - const uint32_t edgeCount = m_indices.size(); - const uint32_t vertexCount = m_positions.size(); - m_oppositeEdges.resize(edgeCount); - m_boundaryEdges.reserve(uint32_t(edgeCount * 0.1f)); - m_isBoundaryVertex.resize(vertexCount); - m_isBoundaryVertex.zeroOutMemory(); - for (uint32_t i = 0; i < edgeCount; i++) - m_oppositeEdges[i] = UINT32_MAX; - const uint32_t faceCount = m_indices.size() / 3; - for (uint32_t i = 0; i < faceCount; i++) { - if (isFaceIgnored(i)) - continue; - for (uint32_t j = 0; j < 3; j++) { - const uint32_t edge = i * 3 + j; - const uint32_t vertex0 = m_indices[edge]; - const uint32_t vertex1 = m_indices[i * 3 + (j + 1) % 3]; - // If there is an edge with opposite winding to this one, the edge isn't on a boundary. - const uint32_t oppositeEdge = findEdge(vertex1, vertex0); - if (oppositeEdge != UINT32_MAX) { - m_oppositeEdges[edge] = oppositeEdge; - } else { - m_boundaryEdges.push_back(edge); - m_isBoundaryVertex.set(vertex0); - m_isBoundaryVertex.set(vertex1); - } - } - } - } - - void linkBoundaries() - { - const uint32_t edgeCount = m_indices.size(); - HashMap vertexToEdgeMap(MemTag::Mesh, edgeCount); // Edge is index / 2 - for (uint32_t i = 0; i < edgeCount; i++) { - vertexToEdgeMap.add(m_indices[meshEdgeIndex0(i)]); - vertexToEdgeMap.add(m_indices[meshEdgeIndex1(i)]); - } - m_nextBoundaryEdges.resize(edgeCount); - for (uint32_t i = 0; i < edgeCount; i++) - m_nextBoundaryEdges[i] = UINT32_MAX; - uint32_t numBoundaryLoops = 0, numUnclosedBoundaries = 0; - BitArray linkedEdges(edgeCount); - linkedEdges.zeroOutMemory(); - for (;;) { - // Find the first boundary edge that hasn't been linked yet. - uint32_t firstEdge = UINT32_MAX; - for (uint32_t i = 0; i < edgeCount; i++) { - if (m_oppositeEdges[i] == UINT32_MAX && !linkedEdges.get(i)) { - firstEdge = i; + // Find the first boundary edge that hasn't been linked yet. + uint32_t firstEdge = UINT32_MAX; + for (uint32_t i = 0; i < edgeCount; i++) { + if (m_oppositeEdges[i] == UINT32_MAX && !linkedEdges.get(i)) { + firstEdge = i; break; } } @@ -2783,6 +2711,15 @@ public: return result; } + // Edge map can be destroyed when no longer used to reduce memory usage. It's used by: + // * Mesh::createBoundaries() + // * Mesh::ColocalEdgeIterator (used by MeshFaceGroups) + // * meshCloseHole() + void destroyEdgeMap() + { + m_edgeMap.destroy(); + } + #if XA_DEBUG_EXPORT_OBJ void writeObjVertices(FILE *file) const { @@ -2796,11 +2733,11 @@ public: fprintf(file, "vt %g %g\n", m_texcoords[i].x, m_texcoords[i].y); } - void writeObjFace(FILE *file, uint32_t face) const + void writeObjFace(FILE *file, uint32_t face, uint32_t offset = 0) const { fprintf(file, "f "); for (uint32_t j = 0; j < 3; j++) { - const uint32_t index = m_indices[face * 3 + j] + 1; // 1-indexed + const uint32_t index = m_indices[face * 3 + j] + 1 + offset; // 1-indexed fprintf(file, "%d/%d/%d%c", index, index, index, j == 2 ? '\n' : ' '); } } @@ -2866,12 +2803,13 @@ public: return area; } + // Returned value is always positive, even if some triangles are flipped. float computeParametricArea() const { float area = 0; for (uint32_t f = 0; f < faceCount(); f++) - area += computeFaceParametricArea(f); - return fabsf(area); // May be negative, depends on texcoord winding. + area += fabsf(computeFaceParametricArea(f)); // May be negative, depends on texcoord winding. + return area; } float computeFaceArea(uint32_t face) const @@ -2978,45 +2916,32 @@ public: XA_INLINE bool isBoundaryEdge(uint32_t edge) const { return m_oppositeEdges[edge] == UINT32_MAX; } XA_INLINE const Array &boundaryEdges() const { return m_boundaryEdges; } XA_INLINE bool isBoundaryVertex(uint32_t vertex) const { return m_isBoundaryVertex.get(vertex); } - XA_INLINE uint32_t colocalVertexCount() const { return m_colocalVertexCount; } XA_INLINE uint32_t vertexCount() const { return m_positions.size(); } XA_INLINE uint32_t vertexAt(uint32_t i) const { return m_indices[i]; } XA_INLINE const Vector3 &position(uint32_t vertex) const { return m_positions[vertex]; } + XA_INLINE const Vector3 *positions() const { return m_positions.data(); } XA_INLINE const Vector3 &normal(uint32_t vertex) const { XA_DEBUG_ASSERT(m_flags & MeshFlags::HasNormals); return m_normals[vertex]; } XA_INLINE const Vector2 &texcoord(uint32_t vertex) const { return m_texcoords[vertex]; } XA_INLINE Vector2 &texcoord(uint32_t vertex) { return m_texcoords[vertex]; } XA_INLINE const Vector2 *texcoords() const { return m_texcoords.data(); } XA_INLINE Vector2 *texcoords() { return m_texcoords.data(); } - XA_INLINE uint32_t ignoredFaceCount() const { return m_ignoredFaceCount; } XA_INLINE uint32_t faceCount() const { return m_indices.size() / 3; } - XA_INLINE uint16_t faceGroupAt(uint32_t face) const { XA_DEBUG_ASSERT(m_flags & MeshFlags::HasFaceGroups); return m_faceGroups[face]; } - XA_INLINE uint32_t faceGroupCount() const { XA_DEBUG_ASSERT(m_flags & MeshFlags::HasFaceGroups); return m_faceGroupFaceCounts.size(); } - XA_INLINE uint32_t faceGroupNextFace(uint32_t face) const { XA_DEBUG_ASSERT(m_flags & MeshFlags::HasFaceGroups); return m_faceGroupNextFace[face]; } - XA_INLINE uint32_t faceGroupFaceCount(uint32_t group) const { XA_DEBUG_ASSERT(m_flags & MeshFlags::HasFaceGroups); return m_faceGroupFaceCounts[group]; } XA_INLINE const uint32_t *indices() const { return m_indices.data(); } XA_INLINE uint32_t indexCount() const { return m_indices.size(); } + XA_INLINE bool isFaceIgnored(uint32_t face) const { return (m_flags & MeshFlags::HasIgnoredFaces) && m_faceIgnore[face]; } private: - bool isFaceIgnored(uint32_t face) const { return (m_flags & MeshFlags::HasIgnoredFaces) && m_faceIgnore[face]; } float m_epsilon; uint32_t m_flags; uint32_t m_id; Array m_faceIgnore; - uint32_t m_ignoredFaceCount; Array m_indices; Array m_positions; Array m_normals; Array m_texcoords; - // Populated by createFaceGroups - Array m_faceGroups; - Array m_faceGroupFirstFace; - Array m_faceGroupNextFace; // In: face. Out: the next face in the same group. - Array m_faceGroupFaceCounts; // In: face group. Out: number of faces in the group. - // Populated by createColocals - uint32_t m_colocalVertexCount; Array m_nextColocalVertex; // In: vertex index. Out: the vertex index of the next colocal position. // Populated by createBoundaries @@ -3029,7 +2954,6 @@ private: struct EdgeKey { - EdgeKey() {} EdgeKey(const EdgeKey &k) : v0(k.v0), v1(k.v1) {} EdgeKey(uint32_t v0, uint32_t v1) : v0(v0), v1(v1) {} bool operator==(const EdgeKey &k) const { return v0 == k.v0 && v1 == k.v1; } @@ -3252,19 +3176,123 @@ public: uint32_t m_edge; uint32_t m_relativeEdge; }; +}; + +struct MeshFaceGroups +{ + typedef uint32_t Handle; + static constexpr Handle kInvalid = UINT32_MAX; + + MeshFaceGroups(const Mesh *mesh) : m_mesh(mesh), m_groups(MemTag::Mesh), m_firstFace(MemTag::Mesh), m_nextFace(MemTag::Mesh), m_faceCount(MemTag::Mesh) {} + XA_INLINE Handle groupAt(uint32_t face) const { return m_groups[face]; } + XA_INLINE uint32_t groupCount() const { return m_faceCount.size(); } + XA_INLINE uint32_t nextFace(uint32_t face) const { return m_nextFace[face]; } + XA_INLINE uint32_t faceCount(uint32_t group) const { return m_faceCount[group]; } + + void compute() + { + m_groups.resize(m_mesh->faceCount()); + m_groups.fillBytes(0xff); // Set all faces to kInvalid + uint32_t firstUnassignedFace = 0; + Handle group = 0; + Array growFaces; + const uint32_t n = m_mesh->faceCount(); + m_nextFace.resize(n); + for (;;) { + // Find an unassigned face. + uint32_t face = UINT32_MAX; + for (uint32_t f = firstUnassignedFace; f < n; f++) { + if (m_groups[f] == kInvalid && !m_mesh->isFaceIgnored(f)) { + face = f; + firstUnassignedFace = f + 1; + break; + } + } + if (face == UINT32_MAX) + break; // All faces assigned to a group (except ignored faces). + m_groups[face] = group; + m_nextFace[face] = UINT32_MAX; + m_firstFace.push_back(face); + growFaces.clear(); + growFaces.push_back(face); + uint32_t prevFace = face, groupFaceCount = 1; + // Find faces connected to the face and assign them to the same group as the face, unless they are already assigned to another group. + for (;;) { + if (growFaces.isEmpty()) + break; + const uint32_t f = growFaces.back(); + growFaces.pop_back(); + for (Mesh::FaceEdgeIterator edgeIt(m_mesh, f); !edgeIt.isDone(); edgeIt.advance()) { + // Iterate opposite edges. There may be more than one - non-manifold geometry can have duplicate edges. + // Prioritize the one with exact vertex match, not just colocal. + // If *any* of the opposite edges are already assigned to this group, don't do anything. + bool alreadyAssignedToThisGroup = false; + uint32_t bestConnectedFace = UINT32_MAX; + for (Mesh::ColocalEdgeIterator oppositeEdgeIt(m_mesh, edgeIt.vertex1(), edgeIt.vertex0()); !oppositeEdgeIt.isDone(); oppositeEdgeIt.advance()) { + const uint32_t oppositeEdge = oppositeEdgeIt.edge(); + const uint32_t oppositeFace = meshEdgeFace(oppositeEdge); +#if 0 + // Reject opposite face if dihedral angle >= 90 degrees. + { + Vector3 a = m_mesh->computeFaceNormal(f); + Vector3 b = m_mesh->computeFaceNormal(oppositeFace); + if (dot(a, b) <= 0.0f) + continue; + } +#endif + if (m_mesh->isFaceIgnored(oppositeFace)) + continue; // Don't add ignored faces to group. + if (m_groups[oppositeFace] == group) { + alreadyAssignedToThisGroup = true; + break; + } + if (m_groups[oppositeFace] != kInvalid) + continue; // Connected face is already assigned to another group. + if (faceDuplicatesGroupEdge(group, oppositeFace)) + continue; // Don't want duplicate edges in a group. + const uint32_t oppositeVertex0 = m_mesh->vertexAt(meshEdgeIndex0(oppositeEdge)); + const uint32_t oppositeVertex1 = m_mesh->vertexAt(meshEdgeIndex1(oppositeEdge)); + if (bestConnectedFace == UINT32_MAX || (oppositeVertex0 == edgeIt.vertex1() && oppositeVertex1 == edgeIt.vertex0())) + bestConnectedFace = oppositeFace; +#if 0 + else { + // Choose the opposite face with the smallest dihedral angle. + const float d1 = 1.0f - dot(computeFaceNormal(f), computeFaceNormal(bestConnectedFace)); + const float d2 = 1.0f - dot(computeFaceNormal(f), computeFaceNormal(oppositeFace)); + if (d2 < d1) + bestConnectedFace = oppositeFace; + } +#endif + } + if (!alreadyAssignedToThisGroup && bestConnectedFace != UINT32_MAX) { + m_groups[bestConnectedFace] = group; + m_nextFace[bestConnectedFace] = UINT32_MAX; + if (prevFace != UINT32_MAX) + m_nextFace[prevFace] = bestConnectedFace; + prevFace = bestConnectedFace; + groupFaceCount++; + growFaces.push_back(bestConnectedFace); + } + } + } + m_faceCount.push_back(groupFaceCount); + group++; + XA_ASSERT(group < kInvalid); + } + } - class GroupFaceIterator + class Iterator { public: - GroupFaceIterator(const Mesh *mesh, uint32_t group) : m_mesh(mesh) + Iterator(const MeshFaceGroups *meshFaceGroups, Handle group) : m_meshFaceGroups(meshFaceGroups) { - XA_DEBUG_ASSERT(group != UINT32_MAX); - m_current = mesh->m_faceGroupFirstFace[group]; + XA_DEBUG_ASSERT(group != kInvalid); + m_current = m_meshFaceGroups->m_firstFace[group]; } void advance() { - m_current = m_mesh->m_faceGroupNextFace[m_current]; + m_current = m_meshFaceGroups->m_nextFace[m_current]; } bool isDone() const @@ -3278,12 +3306,31 @@ public: } private: - const Mesh *m_mesh; + const MeshFaceGroups *m_meshFaceGroups; uint32_t m_current; }; + +private: + // Check if the face duplicates any edges of any face already in the group. + bool faceDuplicatesGroupEdge(Handle group, uint32_t face) const + { + for (Mesh::FaceEdgeIterator edgeIt(m_mesh, face); !edgeIt.isDone(); edgeIt.advance()) { + for (Mesh::ColocalEdgeIterator colocalEdgeIt(m_mesh, edgeIt.vertex0(), edgeIt.vertex1()); !colocalEdgeIt.isDone(); colocalEdgeIt.advance()) { + if (m_groups[meshEdgeFace(colocalEdgeIt.edge())] == group) + return true; + } + } + return false; + } + + const Mesh *m_mesh; + Array m_groups; + Array m_firstFace; + Array m_nextFace; // In: face. Out: the next face in the same group. + Array m_faceCount; // In: face group. Out: number of faces in the group. }; -constexpr uint16_t Mesh::kInvalidFaceGroup; +constexpr MeshFaceGroups::Handle MeshFaceGroups::kInvalid; static bool meshCloseHole(Mesh *mesh, const Array &holeVertices, const Vector3 &normal) { @@ -3748,7 +3795,7 @@ public: return handle; } - void run(TaskGroupHandle handle, Task task) + void run(TaskGroupHandle handle, const Task &task) { XA_DEBUG_ASSERT(handle.value != UINT32_MAX); TaskGroup &group = m_groups[handle.value]; @@ -3810,9 +3857,9 @@ private: }; TaskGroup *m_groups; - uint32_t m_maxGroups; Array m_workers; std::atomic m_shutdown; + uint32_t m_maxGroups; static thread_local uint32_t m_threadIndex; static void workerThread(TaskScheduler *scheduler, Worker *worker, uint32_t threadIndex) @@ -4024,7 +4071,7 @@ public: bool intersect(Vector2 v1, Vector2 v2, float epsilon) { const uint32_t edgeCount = m_edges.size(); - bool bruteForce = edgeCount <= 64; + bool bruteForce = edgeCount <= 20; if (!bruteForce && m_cellDataOffsets.isEmpty()) bruteForce = !createGrid(); if (bruteForce) { @@ -4048,31 +4095,72 @@ public: return false; } - bool intersectSelf(float epsilon) + // If edges is empty, checks for intersection with all edges in the grid. + bool intersect(float epsilon, ConstArrayView edges = ConstArrayView(), ConstArrayView ignoreEdges = ConstArrayView()) { - const uint32_t edgeCount = m_edges.size(); - bool bruteForce = edgeCount <= 64; + bool bruteForce = m_edges.size() <= 20; if (!bruteForce && m_cellDataOffsets.isEmpty()) bruteForce = !createGrid(); - for (uint32_t i = 0; i < edgeCount; i++) { - const uint32_t edge1 = m_edges[i]; + const uint32_t *edges1, *edges2 = nullptr; + uint32_t edges1Count, edges2Count = 0; + if (edges.length == 0) { + edges1 = m_edges.data(); + edges1Count = m_edges.size(); + } else { + edges1 = edges.data; + edges1Count = edges.length; + } + if (bruteForce) { + edges2 = m_edges.data(); + edges2Count = m_edges.size(); + } + for (uint32_t i = 0; i < edges1Count; i++) { + const uint32_t edge1 = edges1[i]; + const uint32_t edge1Vertex[2] = { vertexAt(meshEdgeIndex0(edge1)), vertexAt(meshEdgeIndex1(edge1)) }; + const Vector2 &edge1Position1 = m_positions[edge1Vertex[0]]; + const Vector2 &edge1Position2 = m_positions[edge1Vertex[1]]; + const Extents2 edge1Extents(edge1Position1, edge1Position2); + uint32_t j = 0; if (bruteForce) { - for (uint32_t j = 0; j < edgeCount; j++) { - const uint32_t edge2 = m_edges[j]; - if (edgesIntersect(edge1, edge2, epsilon)) - return true; + // If checking against self, test each edge pair only once. + if (edges.length == 0) { + j = i + 1; + if (j == edges1Count) + break; } } else { computePotentialEdges(edgePosition0(edge1), edgePosition1(edge1)); - uint32_t prevEdge = UINT32_MAX; - for (uint32_t j = 0; j < m_potentialEdges.size(); j++) { - const uint32_t edge2 = m_potentialEdges[j]; - if (edge2 == prevEdge) - continue; - if (edgesIntersect(edge1, edge2, epsilon)) - return true; - prevEdge = edge2; + edges2 = m_potentialEdges.data(); + edges2Count = m_potentialEdges.size(); + } + uint32_t prevEdge = UINT32_MAX; // Handle potential edges duplicates. + for (; j < edges2Count; j++) { + const uint32_t edge2 = edges2[j]; + if (edge1 == edge2) + continue; + if (edge2 == prevEdge) + continue; + prevEdge = edge2; + // Check if edge2 is ignored. + bool ignore = false; + for (uint32_t k = 0; k < ignoreEdges.length; k++) { + if (edge2 == ignoreEdges[k]) { + ignore = true; + break; + } } + if (ignore) + continue; + const uint32_t edge2Vertex[2] = { vertexAt(meshEdgeIndex0(edge2)), vertexAt(meshEdgeIndex1(edge2)) }; + // Ignore connected edges, since they can't intersect (only overlap), and may be detected as false positives. + if (edge1Vertex[0] == edge2Vertex[0] || edge1Vertex[0] == edge2Vertex[1] || edge1Vertex[1] == edge2Vertex[0] || edge1Vertex[1] == edge2Vertex[1]) + continue; + const Vector2 &edge2Position1 = m_positions[edge2Vertex[0]]; + const Vector2 &edge2Position2 = m_positions[edge2Vertex[1]]; + if (!Extents2::intersect(edge1Extents, Extents2(edge2Position1, edge2Position2))) + continue; + if (linesIntersect(edge1Position1, edge1Position2, edge2Position1, edge2Position2, epsilon)) + return true; } } return false; @@ -4218,11 +4306,11 @@ private: } if (currentCell[0] >= m_gridWidth || currentCell[1] >= m_gridHeight) break; - if (stepX == 0 && currentCell[0] < lastCell[0]) + if (stepX == -1 && currentCell[0] < lastCell[0]) break; if (stepX == 1 && currentCell[0] > lastCell[0]) break; - if (stepY == 0 && currentCell[1] < lastCell[1]) + if (stepY == -1 && currentCell[1] < lastCell[1]) break; if (stepY == 1 && currentCell[1] > lastCell[1]) break; @@ -4230,18 +4318,6 @@ private: } } - bool edgesIntersect(uint32_t edge1, uint32_t edge2, float epsilon) const - { - if (edge1 == edge2) - return false; - const uint32_t ai[2] = { vertexAt(meshEdgeIndex0(edge1)), vertexAt(meshEdgeIndex1(edge1)) }; - const uint32_t bi[2] = { vertexAt(meshEdgeIndex0(edge2)), vertexAt(meshEdgeIndex1(edge2)) }; - // Ignore connected edges, since they can't intersect (only overlap), and may be detected as false positives. - if (ai[0] == bi[0] || ai[0] == bi[1] || ai[1] == bi[0] || ai[1] == bi[1]) - return false; - return linesIntersect(m_positions[ai[0]], m_positions[ai[1]], m_positions[bi[0]], m_positions[bi[1]], epsilon); - } - uint32_t cellX(float x) const { return min((uint32_t)max(0.0f, (x - m_gridOrigin.x) / m_cellSize), m_gridWidth - 1u); @@ -4301,50 +4377,852 @@ struct UvMeshInstance bool rotateCharts; }; -namespace raster { -class ClippedTriangle -{ -public: - ClippedTriangle(const Vector2 &a, const Vector2 &b, const Vector2 &c) - { - m_numVertices = 3; - m_activeVertexBuffer = 0; - m_verticesA[0] = a; - m_verticesA[1] = b; - m_verticesA[2] = c; - m_vertexBuffers[0] = m_verticesA; - m_vertexBuffers[1] = m_verticesB; - } +/* + * Copyright (c) 2004-2010, Bruno Levy + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * * Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * * Neither the name of the ALICE Project-Team nor the names of its + * contributors may be used to endorse or promote products derived from this + * software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * + * If you modify this software, you should include a notice giving the + * name of the person performing the modification, the date of modification, + * and the reason for such modification. + * + * Contact: Bruno Levy + * + * levy@loria.fr + * + * ALICE Project + * LORIA, INRIA Lorraine, + * Campus Scientifique, BP 239 + * 54506 VANDOEUVRE LES NANCY CEDEX + * FRANCE + */ +namespace opennl { +#define NL_NEW(T) XA_ALLOC(MemTag::OpenNL, T) +#define NL_NEW_ARRAY(T,NB) XA_ALLOC_ARRAY(MemTag::OpenNL, T, NB) +#define NL_RENEW_ARRAY(T,x,NB) XA_REALLOC(MemTag::OpenNL, x, T, NB) +#define NL_DELETE(x) XA_FREE(x); x = nullptr +#define NL_DELETE_ARRAY(x) XA_FREE(x); x = nullptr +#define NL_CLEAR(x, T) memset(x, 0, sizeof(T)); +#define NL_CLEAR_ARRAY(T,x,NB) memset(x, 0, (size_t)(NB)*sizeof(T)) +#define NL_NEW_VECTOR(dim) XA_ALLOC_ARRAY(MemTag::OpenNL, double, dim) +#define NL_DELETE_VECTOR(ptr) XA_FREE(ptr) + +struct NLMatrixStruct; +typedef NLMatrixStruct * NLMatrix; +typedef void (*NLDestroyMatrixFunc)(NLMatrix M); +typedef void (*NLMultMatrixVectorFunc)(NLMatrix M, const double* x, double* y); + +#define NL_MATRIX_SPARSE_DYNAMIC 0x1001 +#define NL_MATRIX_CRS 0x1002 +#define NL_MATRIX_OTHER 0x1006 + +struct NLMatrixStruct +{ + uint32_t m; + uint32_t n; + uint32_t type; + NLDestroyMatrixFunc destroy_func; + NLMultMatrixVectorFunc mult_func; +}; - void clipHorizontalPlane(float offset, float clipdirection) - { - Vector2 *v = m_vertexBuffers[m_activeVertexBuffer]; - m_activeVertexBuffer ^= 1; - Vector2 *v2 = m_vertexBuffers[m_activeVertexBuffer]; - v[m_numVertices] = v[0]; - float dy2, dy1 = offset - v[0].y; - int dy2in, dy1in = clipdirection * dy1 >= 0; - uint32_t p = 0; - for (uint32_t k = 0; k < m_numVertices; k++) { - dy2 = offset - v[k + 1].y; - dy2in = clipdirection * dy2 >= 0; - if (dy1in) v2[p++] = v[k]; - if ( dy1in + dy2in == 1 ) { // not both in/out - float dx = v[k + 1].x - v[k].x; - float dy = v[k + 1].y - v[k].y; - v2[p++] = Vector2(v[k].x + dy1 * (dx / dy), offset); - } - dy1 = dy2; - dy1in = dy2in; - } - m_numVertices = p; - } +/* Dynamic arrays for sparse row/columns */ - void clipVerticalPlane(float offset, float clipdirection) - { - Vector2 *v = m_vertexBuffers[m_activeVertexBuffer]; - m_activeVertexBuffer ^= 1; - Vector2 *v2 = m_vertexBuffers[m_activeVertexBuffer]; +struct NLCoeff +{ + uint32_t index; + double value; +}; + +struct NLRowColumn +{ + uint32_t size; + uint32_t capacity; + NLCoeff* coeff; +}; + +/* Compressed Row Storage */ + +struct NLCRSMatrix +{ + uint32_t m; + uint32_t n; + uint32_t type; + NLDestroyMatrixFunc destroy_func; + NLMultMatrixVectorFunc mult_func; + double* val; + uint32_t* rowptr; + uint32_t* colind; + uint32_t nslices; + uint32_t* sliceptr; +}; + +/* SparseMatrix data structure */ + +struct NLSparseMatrix +{ + uint32_t m; + uint32_t n; + uint32_t type; + NLDestroyMatrixFunc destroy_func; + NLMultMatrixVectorFunc mult_func; + uint32_t diag_size; + uint32_t diag_capacity; + NLRowColumn* row; + NLRowColumn* column; + double* diag; + uint32_t row_capacity; + uint32_t column_capacity; +}; + +/* NLContext data structure */ + +struct NLBufferBinding +{ + void* base_address; + uint32_t stride; +}; + +#define NL_BUFFER_ITEM(B,i) *(double*)((void*)((char*)((B).base_address)+((i)*(B).stride))) + +struct NLContext +{ + NLBufferBinding *variable_buffer; + double *variable_value; + bool *variable_is_locked; + uint32_t *variable_index; + uint32_t n; + NLMatrix M; + NLMatrix P; + NLMatrix B; + NLRowColumn af; + NLRowColumn al; + double *x; + double *b; + uint32_t nb_variables; + uint32_t nb_systems; + uint32_t current_row; + uint32_t max_iterations; + bool max_iterations_defined; + double threshold; + double omega; + uint32_t used_iterations; + double error; +}; + +static void nlDeleteMatrix(NLMatrix M) +{ + if (!M) + return; + M->destroy_func(M); + NL_DELETE(M); +} + +static void nlMultMatrixVector(NLMatrix M, const double* x, double* y) +{ + M->mult_func(M, x, y); +} + +static void nlRowColumnConstruct(NLRowColumn* c) +{ + c->size = 0; + c->capacity = 0; + c->coeff = nullptr; +} + +static void nlRowColumnDestroy(NLRowColumn* c) +{ + NL_DELETE_ARRAY(c->coeff); + c->size = 0; + c->capacity = 0; +} + +static void nlRowColumnGrow(NLRowColumn* c) +{ + if (c->capacity != 0) { + c->capacity = 2 * c->capacity; + c->coeff = NL_RENEW_ARRAY(NLCoeff, c->coeff, c->capacity); + } else { + c->capacity = 4; + c->coeff = NL_NEW_ARRAY(NLCoeff, c->capacity); + NL_CLEAR_ARRAY(NLCoeff, c->coeff, c->capacity); + } +} + +static void nlRowColumnAdd(NLRowColumn* c, uint32_t index, double value) +{ + for (uint32_t i = 0; i < c->size; i++) { + if (c->coeff[i].index == index) { + c->coeff[i].value += value; + return; + } + } + if (c->size == c->capacity) + nlRowColumnGrow(c); + c->coeff[c->size].index = index; + c->coeff[c->size].value = value; + c->size++; +} + +/* Does not check whether the index already exists */ +static void nlRowColumnAppend(NLRowColumn* c, uint32_t index, double value) +{ + if (c->size == c->capacity) + nlRowColumnGrow(c); + c->coeff[c->size].index = index; + c->coeff[c->size].value = value; + c->size++; +} + +static void nlRowColumnZero(NLRowColumn* c) +{ + c->size = 0; +} + +static void nlRowColumnClear(NLRowColumn* c) +{ + c->size = 0; + c->capacity = 0; + NL_DELETE_ARRAY(c->coeff); +} + +static int nlCoeffCompare(const void* p1, const void* p2) +{ + return (((NLCoeff*)(p2))->index < ((NLCoeff*)(p1))->index); +} + +static void nlRowColumnSort(NLRowColumn* c) +{ + qsort(c->coeff, c->size, sizeof(NLCoeff), nlCoeffCompare); +} + +/* CRSMatrix data structure */ + +static void nlCRSMatrixDestroy(NLCRSMatrix* M) +{ + NL_DELETE_ARRAY(M->val); + NL_DELETE_ARRAY(M->rowptr); + NL_DELETE_ARRAY(M->colind); + NL_DELETE_ARRAY(M->sliceptr); + M->m = 0; + M->n = 0; + M->nslices = 0; +} + +static void nlCRSMatrixMultSlice(NLCRSMatrix* M, const double* x, double* y, uint32_t Ibegin, uint32_t Iend) +{ + for (uint32_t i = Ibegin; i < Iend; ++i) { + double sum = 0.0; + for (uint32_t j = M->rowptr[i]; j < M->rowptr[i + 1]; ++j) + sum += M->val[j] * x[M->colind[j]]; + y[i] = sum; + } +} + +static void nlCRSMatrixMult(NLCRSMatrix* M, const double* x, double* y) +{ + int nslices = (int)(M->nslices); + for (int slice = 0; slice < nslices; ++slice) + nlCRSMatrixMultSlice(M, x, y, M->sliceptr[slice], M->sliceptr[slice + 1]); +} + +static void nlCRSMatrixConstruct(NLCRSMatrix* M, uint32_t m, uint32_t n, uint32_t nnz, uint32_t nslices) +{ + M->m = m; + M->n = n; + M->type = NL_MATRIX_CRS; + M->destroy_func = (NLDestroyMatrixFunc)nlCRSMatrixDestroy; + M->mult_func = (NLMultMatrixVectorFunc)nlCRSMatrixMult; + M->nslices = nslices; + M->val = NL_NEW_ARRAY(double, nnz); + NL_CLEAR_ARRAY(double, M->val, nnz); + M->rowptr = NL_NEW_ARRAY(uint32_t, m + 1); + NL_CLEAR_ARRAY(uint32_t, M->rowptr, m + 1); + M->colind = NL_NEW_ARRAY(uint32_t, nnz); + NL_CLEAR_ARRAY(uint32_t, M->colind, nnz); + M->sliceptr = NL_NEW_ARRAY(uint32_t, nslices + 1); + NL_CLEAR_ARRAY(uint32_t, M->sliceptr, nslices + 1); +} + +/* SparseMatrix data structure */ + +static void nlSparseMatrixDestroyRowColumns(NLSparseMatrix* M) +{ + for (uint32_t i = 0; i < M->m; i++) + nlRowColumnDestroy(&(M->row[i])); + NL_DELETE_ARRAY(M->row); +} + +static void nlSparseMatrixDestroy(NLSparseMatrix* M) +{ + XA_DEBUG_ASSERT(M->type == NL_MATRIX_SPARSE_DYNAMIC); + nlSparseMatrixDestroyRowColumns(M); + NL_DELETE_ARRAY(M->diag); +} + +static void nlSparseMatrixAdd(NLSparseMatrix* M, uint32_t i, uint32_t j, double value) +{ + XA_DEBUG_ASSERT(i >= 0 && i <= M->m - 1); + XA_DEBUG_ASSERT(j >= 0 && j <= M->n - 1); + if (i == j) + M->diag[i] += value; + nlRowColumnAdd(&(M->row[i]), j, value); +} + +/* Returns the number of non-zero coefficients */ +static uint32_t nlSparseMatrixNNZ(NLSparseMatrix* M) +{ + uint32_t nnz = 0; + for (uint32_t i = 0; i < M->m; i++) + nnz += M->row[i].size; + return nnz; +} + +static void nlSparseMatrixSort(NLSparseMatrix* M) +{ + for (uint32_t i = 0; i < M->m; i++) + nlRowColumnSort(&(M->row[i])); +} + +/* SparseMatrix x Vector routines, internal helper routines */ + +static void nlSparseMatrix_mult_rows(NLSparseMatrix* A, const double* x, double* y) +{ + /* + * Note: OpenMP does not like unsigned ints + * (causes some floating point exceptions), + * therefore I use here signed ints for all + * indices. + */ + int m = (int)(A->m); + NLCoeff* c = nullptr; + NLRowColumn* Ri = nullptr; + for (int i = 0; i < m; i++) { + Ri = &(A->row[i]); + y[i] = 0; + for (int ij = 0; ij < (int)(Ri->size); ij++) { + c = &(Ri->coeff[ij]); + y[i] += c->value * x[c->index]; + } + } +} + +static void nlSparseMatrixMult(NLSparseMatrix* A, const double* x, double* y) +{ + XA_DEBUG_ASSERT(A->type == NL_MATRIX_SPARSE_DYNAMIC); + nlSparseMatrix_mult_rows(A, x, y); +} + +static void nlSparseMatrixConstruct(NLSparseMatrix* M, uint32_t m, uint32_t n) +{ + M->m = m; + M->n = n; + M->type = NL_MATRIX_SPARSE_DYNAMIC; + M->destroy_func = (NLDestroyMatrixFunc)nlSparseMatrixDestroy; + M->mult_func = (NLMultMatrixVectorFunc)nlSparseMatrixMult; + M->row = NL_NEW_ARRAY(NLRowColumn, m); + NL_CLEAR_ARRAY(NLRowColumn, M->row, m); + M->row_capacity = m; + for (uint32_t i = 0; i < n; i++) + nlRowColumnConstruct(&(M->row[i])); + M->row_capacity = 0; + M->column = nullptr; + M->column_capacity = 0; + M->diag_size = min(m, n); + M->diag_capacity = M->diag_size; + M->diag = NL_NEW_ARRAY(double, M->diag_size); + NL_CLEAR_ARRAY(double, M->diag, M->diag_size); +} + +static NLMatrix nlCRSMatrixNewFromSparseMatrix(NLSparseMatrix* M) +{ + uint32_t nnz = nlSparseMatrixNNZ(M); + uint32_t nslices = 8; /* TODO: get number of cores */ + uint32_t slice, cur_bound, cur_NNZ, cur_row; + uint32_t k; + uint32_t slice_size = nnz / nslices; + NLCRSMatrix* CRS = NL_NEW(NLCRSMatrix); + NL_CLEAR(CRS, NLCRSMatrix); + nlCRSMatrixConstruct(CRS, M->m, M->n, nnz, nslices); + nlSparseMatrixSort(M); + /* Convert matrix to CRS format */ + k = 0; + for (uint32_t i = 0; i < M->m; ++i) { + NLRowColumn* Ri = &(M->row[i]); + CRS->rowptr[i] = k; + for (uint32_t ij = 0; ij < Ri->size; ij++) { + NLCoeff* c = &(Ri->coeff[ij]); + CRS->val[k] = c->value; + CRS->colind[k] = c->index; + ++k; + } + } + CRS->rowptr[M->m] = k; + /* Create "slices" to be used by parallel sparse matrix vector product */ + if (CRS->sliceptr) { + cur_bound = slice_size; + cur_NNZ = 0; + cur_row = 0; + CRS->sliceptr[0] = 0; + for (slice = 1; slice < nslices; ++slice) { + while (cur_NNZ < cur_bound && cur_row < M->m) { + ++cur_row; + cur_NNZ += CRS->rowptr[cur_row + 1] - CRS->rowptr[cur_row]; + } + CRS->sliceptr[slice] = cur_row; + cur_bound += slice_size; + } + CRS->sliceptr[nslices] = M->m; + } + return (NLMatrix)CRS; +} + +static void nlMatrixCompress(NLMatrix* M) +{ + NLMatrix CRS = nullptr; + if ((*M)->type != NL_MATRIX_SPARSE_DYNAMIC) + return; + CRS = nlCRSMatrixNewFromSparseMatrix((NLSparseMatrix*)*M); + nlDeleteMatrix(*M); + *M = CRS; +} + +static NLContext *nlNewContext() +{ + NLContext* result = NL_NEW(NLContext); + NL_CLEAR(result, NLContext); + result->max_iterations = 100; + result->threshold = 1e-6; + result->omega = 1.5; + result->nb_systems = 1; + return result; +} + +static void nlDeleteContext(NLContext *context) +{ + nlDeleteMatrix(context->M); + context->M = nullptr; + nlDeleteMatrix(context->P); + context->P = nullptr; + nlDeleteMatrix(context->B); + context->B = nullptr; + nlRowColumnDestroy(&context->af); + nlRowColumnDestroy(&context->al); + NL_DELETE_ARRAY(context->variable_value); + NL_DELETE_ARRAY(context->variable_buffer); + NL_DELETE_ARRAY(context->variable_is_locked); + NL_DELETE_ARRAY(context->variable_index); + NL_DELETE_ARRAY(context->x); + NL_DELETE_ARRAY(context->b); + NL_DELETE(context); +} + +static double ddot(int n, const double *x, const double *y) +{ + double sum = 0.0; + for (int i = 0; i < n; i++) + sum += x[i] * y[i]; + return sum; +} + +static void daxpy(int n, double a, const double *x, double *y) +{ + for (int i = 0; i < n; i++) + y[i] = a * x[i] + y[i]; +} + +static void dscal(int n, double a, double *x) +{ + for (int i = 0; i < n; i++) + x[i] *= a; +} + +/* + * The implementation of the solvers is inspired by + * the lsolver library, by Christian Badura, available from: + * http://www.mathematik.uni-freiburg.de + * /IAM/Research/projectskr/lin_solver/ + * + * About the Conjugate Gradient, details can be found in: + * Ashby, Manteuffel, Saylor + * A taxononmy for conjugate gradient methods + * SIAM J Numer Anal 27, 1542-1568 (1990) + * + * This version is completely abstract, the same code can be used for + * CPU/GPU, dense matrix / sparse matrix etc... + * Abstraction is realized through: + * - Abstract matrix interface (NLMatrix), that can implement different + * versions of matrix x vector product (CPU/GPU, sparse/dense ...) + */ + +static uint32_t nlSolveSystem_PRE_CG(NLMatrix M, NLMatrix P, double* b, double* x, double eps, uint32_t max_iter, double *sq_bnorm, double *sq_rnorm) +{ + int N = (int)M->n; + double* r = NL_NEW_VECTOR(N); + double* d = NL_NEW_VECTOR(N); + double* h = NL_NEW_VECTOR(N); + double *Ad = h; + uint32_t its = 0; + double rh, alpha, beta; + double b_square = ddot(N, b, b); + double err = eps * eps*b_square; + double curr_err; + nlMultMatrixVector(M, x, r); + daxpy(N, -1., b, r); + nlMultMatrixVector(P, r, d); + memcpy(h, d, N * sizeof(double)); + rh = ddot(N, r, h); + curr_err = ddot(N, r, r); + while (curr_err > err && its < max_iter) { + nlMultMatrixVector(M, d, Ad); + alpha = rh / ddot(N, d, Ad); + daxpy(N, -alpha, d, x); + daxpy(N, -alpha, Ad, r); + nlMultMatrixVector(P, r, h); + beta = 1. / rh; + rh = ddot(N, r, h); + beta *= rh; + dscal(N, beta, d); + daxpy(N, 1., h, d); + ++its; + curr_err = ddot(N, r, r); + } + NL_DELETE_VECTOR(r); + NL_DELETE_VECTOR(d); + NL_DELETE_VECTOR(h); + *sq_bnorm = b_square; + *sq_rnorm = curr_err; + return its; +} + +static uint32_t nlSolveSystemIterative(NLContext *context, NLMatrix M, NLMatrix P, double* b_in, double* x_in, double eps, uint32_t max_iter) +{ + uint32_t result = 0; + double rnorm = 0.0; + double bnorm = 0.0; + double* b = b_in; + double* x = x_in; + XA_DEBUG_ASSERT(M->m == M->n); + double sq_bnorm, sq_rnorm; + result = nlSolveSystem_PRE_CG(M, P, b, x, eps, max_iter, &sq_bnorm, &sq_rnorm); + /* Get residual norm and rhs norm */ + bnorm = sqrt(sq_bnorm); + rnorm = sqrt(sq_rnorm); + if (bnorm == 0.0) + context->error = rnorm; + else + context->error = rnorm / bnorm; + context->used_iterations = result; + return result; +} + +static bool nlSolveIterative(NLContext *context) +{ + double* b = context->b; + double* x = context->x; + uint32_t n = context->n; + NLMatrix M = context->M; + NLMatrix P = context->P; + for (uint32_t k = 0; k < context->nb_systems; ++k) { + nlSolveSystemIterative(context, M, P, b, x, context->threshold, context->max_iterations); + b += n; + x += n; + } + return true; +} + +struct NLJacobiPreconditioner +{ + uint32_t m; + uint32_t n; + uint32_t type; + NLDestroyMatrixFunc destroy_func; + NLMultMatrixVectorFunc mult_func; + double* diag_inv; +}; + +static void nlJacobiPreconditionerDestroy(NLJacobiPreconditioner* M) +{ + NL_DELETE_ARRAY(M->diag_inv); +} + +static void nlJacobiPreconditionerMult(NLJacobiPreconditioner* M, const double* x, double* y) +{ + for (uint32_t i = 0; i < M->n; ++i) + y[i] = x[i] * M->diag_inv[i]; +} + +static NLMatrix nlNewJacobiPreconditioner(NLMatrix M_in) +{ + NLSparseMatrix* M = nullptr; + NLJacobiPreconditioner* result = nullptr; + XA_DEBUG_ASSERT(M_in->type == NL_MATRIX_SPARSE_DYNAMIC); + XA_DEBUG_ASSERT(M_in->m == M_in->n); + M = (NLSparseMatrix*)M_in; + result = NL_NEW(NLJacobiPreconditioner); + NL_CLEAR(result, NLJacobiPreconditioner); + result->m = M->m; + result->n = M->n; + result->type = NL_MATRIX_OTHER; + result->destroy_func = (NLDestroyMatrixFunc)nlJacobiPreconditionerDestroy; + result->mult_func = (NLMultMatrixVectorFunc)nlJacobiPreconditionerMult; + result->diag_inv = NL_NEW_ARRAY(double, M->n); + NL_CLEAR_ARRAY(double, result->diag_inv, M->n); + for (uint32_t i = 0; i < M->n; ++i) + result->diag_inv[i] = (M->diag[i] == 0.0) ? 1.0 : 1.0 / M->diag[i]; + return (NLMatrix)result; +} + +#define NL_NB_VARIABLES 0x101 +#define NL_MAX_ITERATIONS 0x103 + +static void nlSolverParameteri(NLContext *context, uint32_t pname, int param) +{ + if (pname == NL_NB_VARIABLES) { + XA_DEBUG_ASSERT(param > 0); + context->nb_variables = (uint32_t)param; + } else if (pname == NL_MAX_ITERATIONS) { + XA_DEBUG_ASSERT(param > 0); + context->max_iterations = (uint32_t)param; + context->max_iterations_defined = true; + } +} + +static void nlSetVariable(NLContext *context, uint32_t index, double value) +{ + XA_DEBUG_ASSERT(index >= 0 && index <= context->nb_variables - 1); + NL_BUFFER_ITEM(context->variable_buffer[0], index) = value; +} + +static double nlGetVariable(NLContext *context, uint32_t index) +{ + XA_DEBUG_ASSERT(index >= 0 && index <= context->nb_variables - 1); + return NL_BUFFER_ITEM(context->variable_buffer[0], index); +} + +static void nlLockVariable(NLContext *context, uint32_t index) +{ + XA_DEBUG_ASSERT(index >= 0 && index <= context->nb_variables - 1); + context->variable_is_locked[index] = true; +} + +static void nlVariablesToVector(NLContext *context) +{ + uint32_t n = context->n; + XA_DEBUG_ASSERT(context->x); + for (uint32_t k = 0; k < context->nb_systems; ++k) { + for (uint32_t i = 0; i < context->nb_variables; ++i) { + if (!context->variable_is_locked[i]) { + uint32_t index = context->variable_index[i]; + XA_DEBUG_ASSERT(index < context->n); + double value = NL_BUFFER_ITEM(context->variable_buffer[k], i); + context->x[index + k * n] = value; + } + } + } +} + +static void nlVectorToVariables(NLContext *context) +{ + uint32_t n = context->n; + XA_DEBUG_ASSERT(context->x); + for (uint32_t k = 0; k < context->nb_systems; ++k) { + for (uint32_t i = 0; i < context->nb_variables; ++i) { + if (!context->variable_is_locked[i]) { + uint32_t index = context->variable_index[i]; + XA_DEBUG_ASSERT(index < context->n); + double value = context->x[index + k * n]; + NL_BUFFER_ITEM(context->variable_buffer[k], i) = value; + } + } + } +} + +static void nlCoefficient(NLContext *context, uint32_t index, double value) +{ + XA_DEBUG_ASSERT(index >= 0 && index <= context->nb_variables - 1); + if (context->variable_is_locked[index]) { + /* + * Note: in al, indices are NLvariable indices, + * within [0..nb_variables-1] + */ + nlRowColumnAppend(&(context->al), index, value); + } else { + /* + * Note: in af, indices are system indices, + * within [0..n-1] + */ + nlRowColumnAppend(&(context->af), context->variable_index[index], value); + } +} + +#define NL_SYSTEM 0x0 +#define NL_MATRIX 0x1 +#define NL_ROW 0x2 + +static void nlBegin(NLContext *context, uint32_t prim) +{ + if (prim == NL_SYSTEM) { + XA_DEBUG_ASSERT(context->nb_variables > 0); + context->variable_buffer = NL_NEW_ARRAY(NLBufferBinding, context->nb_systems); + NL_CLEAR_ARRAY(NLBufferBinding, context->variable_buffer, context->nb_systems); + context->variable_value = NL_NEW_ARRAY(double, context->nb_variables * context->nb_systems); + NL_CLEAR_ARRAY(double, context->variable_value, context->nb_variables * context->nb_systems); + for (uint32_t k = 0; k < context->nb_systems; ++k) { + context->variable_buffer[k].base_address = + context->variable_value + + k * context->nb_variables; + context->variable_buffer[k].stride = sizeof(double); + } + context->variable_is_locked = NL_NEW_ARRAY(bool, context->nb_variables); + NL_CLEAR_ARRAY(bool, context->variable_is_locked, context->nb_variables); + context->variable_index = NL_NEW_ARRAY(uint32_t, context->nb_variables); + NL_CLEAR_ARRAY(uint32_t, context->variable_index, context->nb_variables); + } else if (prim == NL_MATRIX) { + if (context->M) + return; + uint32_t n = 0; + for (uint32_t i = 0; i < context->nb_variables; i++) { + if (!context->variable_is_locked[i]) { + context->variable_index[i] = n; + n++; + } else + context->variable_index[i] = (uint32_t)~0; + } + context->n = n; + if (!context->max_iterations_defined) + context->max_iterations = n * 5; + context->M = (NLMatrix)(NL_NEW(NLSparseMatrix)); + NL_CLEAR(context->M, NLSparseMatrix); + nlSparseMatrixConstruct((NLSparseMatrix*)(context->M), n, n); + context->x = NL_NEW_ARRAY(double, n*context->nb_systems); + NL_CLEAR_ARRAY(double, context->x, n*context->nb_systems); + context->b = NL_NEW_ARRAY(double, n*context->nb_systems); + NL_CLEAR_ARRAY(double, context->b, n*context->nb_systems); + nlVariablesToVector(context); + nlRowColumnConstruct(&context->af); + nlRowColumnConstruct(&context->al); + context->current_row = 0; + } else if (prim == NL_ROW) { + nlRowColumnZero(&context->af); + nlRowColumnZero(&context->al); + } +} + +static void nlEnd(NLContext *context, uint32_t prim) +{ + if (prim == NL_MATRIX) { + nlRowColumnClear(&context->af); + nlRowColumnClear(&context->al); + } else if (prim == NL_ROW) { + NLRowColumn* af = &context->af; + NLRowColumn* al = &context->al; + NLSparseMatrix* M = (NLSparseMatrix*)context->M; + double* b = context->b; + uint32_t nf = af->size; + uint32_t nl = al->size; + uint32_t n = context->n; + double S; + /* + * least_squares : we want to solve + * A'A x = A'b + */ + for (uint32_t i = 0; i < nf; i++) { + for (uint32_t j = 0; j < nf; j++) { + nlSparseMatrixAdd(M, af->coeff[i].index, af->coeff[j].index, af->coeff[i].value * af->coeff[j].value); + } + } + for (uint32_t k = 0; k < context->nb_systems; ++k) { + S = 0.0; + for (uint32_t jj = 0; jj < nl; ++jj) { + uint32_t j = al->coeff[jj].index; + S += al->coeff[jj].value * NL_BUFFER_ITEM(context->variable_buffer[k], j); + } + for (uint32_t jj = 0; jj < nf; jj++) + b[k*n + af->coeff[jj].index] -= af->coeff[jj].value * S; + } + context->current_row++; + } +} + +static bool nlSolve(NLContext *context) +{ + nlDeleteMatrix(context->P); + context->P = nlNewJacobiPreconditioner(context->M); + nlMatrixCompress(&context->M); + bool result = nlSolveIterative(context); + nlVectorToVariables(context); + return result; +} +} // namespace opennl + +namespace raster { +class ClippedTriangle +{ +public: + ClippedTriangle(const Vector2 &a, const Vector2 &b, const Vector2 &c) + { + m_numVertices = 3; + m_activeVertexBuffer = 0; + m_verticesA[0] = a; + m_verticesA[1] = b; + m_verticesA[2] = c; + m_vertexBuffers[0] = m_verticesA; + m_vertexBuffers[1] = m_verticesB; + m_area = 0; + } + + void clipHorizontalPlane(float offset, float clipdirection) + { + Vector2 *v = m_vertexBuffers[m_activeVertexBuffer]; + m_activeVertexBuffer ^= 1; + Vector2 *v2 = m_vertexBuffers[m_activeVertexBuffer]; + v[m_numVertices] = v[0]; + float dy2, dy1 = offset - v[0].y; + int dy2in, dy1in = clipdirection * dy1 >= 0; + uint32_t p = 0; + for (uint32_t k = 0; k < m_numVertices; k++) { + dy2 = offset - v[k + 1].y; + dy2in = clipdirection * dy2 >= 0; + if (dy1in) v2[p++] = v[k]; + if ( dy1in + dy2in == 1 ) { // not both in/out + float dx = v[k + 1].x - v[k].x; + float dy = v[k + 1].y - v[k].y; + v2[p++] = Vector2(v[k].x + dy1 * (dx / dy), offset); + } + dy1 = dy2; + dy1in = dy2in; + } + m_numVertices = p; + } + + void clipVerticalPlane(float offset, float clipdirection) + { + Vector2 *v = m_vertexBuffers[m_activeVertexBuffer]; + m_activeVertexBuffer ^= 1; + Vector2 *v2 = m_vertexBuffers[m_activeVertexBuffer]; v[m_numVertices] = v[0]; float dx2, dx1 = offset - v[0].x; int dx2in, dx1in = clipdirection * dx1 >= 0; @@ -4409,12 +5287,8 @@ typedef bool (*SamplingCallback)(void *param, int x, int y); /// A triangle for rasterization. struct Triangle { - Triangle(const Vector2 &v0, const Vector2 &v1, const Vector2 &v2) + Triangle(const Vector2 &_v0, const Vector2 &_v1, const Vector2 &_v2) : v1(_v0), v2(_v2), v3(_v1) { - // Init vertices. - this->v1 = v0; - this->v2 = v2; - this->v3 = v1; // make sure every triangle is front facing. flipBackface(); // Compute deltas. @@ -4508,297 +5382,51 @@ struct Triangle } } } - return true; - } - -private: - void flipBackface() - { - // check if triangle is backfacing, if so, swap two vertices - if ( ((v3.x - v1.x) * (v2.y - v1.y) - (v3.y - v1.y) * (v2.x - v1.x)) < 0 ) { - Vector2 hv = v1; - v1 = v2; - v2 = hv; // swap pos - } - } - - // compute unit inward normals for each edge. - void computeUnitInwardNormals() - { - n1 = v1 - v2; - n1 = Vector2(-n1.y, n1.x); - n1 = n1 * (1.0f / sqrtf(dot(n1, n1))); - n2 = v2 - v3; - n2 = Vector2(-n2.y, n2.x); - n2 = n2 * (1.0f / sqrtf(dot(n2, n2))); - n3 = v3 - v1; - n3 = Vector2(-n3.y, n3.x); - n3 = n3 * (1.0f / sqrtf(dot(n3, n3))); - } - - // Vertices. - Vector2 v1, v2, v3; - Vector2 n1, n2, n3; // unit inward normals -}; - -// Process the given triangle. Returns false if rasterization was interrupted by the callback. -static bool drawTriangle(const Vector2 &extents, const Vector2 v[3], SamplingCallback cb, void *param) -{ - Triangle tri(v[0], v[1], v[2]); - // @@ It would be nice to have a conservative drawing mode that enlarges the triangle extents by one texel and is able to handle degenerate triangles. - // @@ Maybe the simplest thing to do would be raster triangle edges. - if (tri.isValid()) - return tri.drawAA(extents, cb, param); - return true; -} - -} // namespace raster - -// Full and sparse vector and matrix classes. BLAS subset. -// Pseudo-BLAS interface. -namespace sparse { - -/** -* Sparse matrix class. The matrix is assumed to be sparse and to have -* very few non-zero elements, for this reason it's stored in indexed -* format. To multiply column vectors efficiently, the matrix stores -* the elements in indexed-column order, there is a list of indexed -* elements for each row of the matrix. As with the FullVector the -* dimension of the matrix is constant. -**/ -class Matrix -{ -public: - // An element of the sparse array. - struct Coefficient - { - uint32_t x; // column - float v; // value - }; - - Matrix(uint32_t d) : m_width(d), m_array(MemTag::Matrix) - { - m_array.resize(d); - m_array.runCtors(); -#if XA_DEBUG_HEAP - for (uint32_t i = 0; i < d; i++) - m_array[i].setMemTag(MemTag::Matrix); -#endif - } - - Matrix(uint32_t w, uint32_t h) : m_width(w), m_array(MemTag::Matrix) - { - m_array.resize(h); - m_array.runCtors(); -#if XA_DEBUG_HEAP - for (uint32_t i = 0; i < h; i++) - m_array[i].setMemTag(MemTag::Matrix); -#endif - } - - ~Matrix() - { - m_array.runDtors(); - } - - Matrix(const Matrix &m) = delete; - Matrix &operator=(const Matrix &m) = delete; - uint32_t width() const { return m_width; } - uint32_t height() const { return m_array.size(); } - bool isSquare() const { return width() == height(); } - - // x is column, y is row - float getCoefficient(uint32_t x, uint32_t y) const - { - XA_DEBUG_ASSERT( x < width() ); - XA_DEBUG_ASSERT( y < height() ); - const uint32_t count = m_array[y].size(); - for (uint32_t i = 0; i < count; i++) { - if (m_array[y][i].x == x) return m_array[y][i].v; - } - return 0.0f; - } - - void setCoefficient(uint32_t x, uint32_t y, float f) - { - XA_DEBUG_ASSERT( x < width() ); - XA_DEBUG_ASSERT( y < height() ); - const uint32_t count = m_array[y].size(); - for (uint32_t i = 0; i < count; i++) { - if (m_array[y][i].x == x) { - m_array[y][i].v = f; - return; - } - } - if (f != 0.0f) { - Coefficient c = { x, f }; - m_array[y].push_back( c ); - } - } - - float dotRow(uint32_t y, const FullVector &v) const - { - XA_DEBUG_ASSERT( y < height() ); - const uint32_t count = m_array[y].size(); - float sum = 0; - for (uint32_t i = 0; i < count; i++) { - sum += m_array[y][i].v * v[m_array[y][i].x]; - } - return sum; - } - - void madRow(uint32_t y, float alpha, FullVector &v) const - { - XA_DEBUG_ASSERT(y < height()); - const uint32_t count = m_array[y].size(); - for (uint32_t i = 0; i < count; i++) { - v[m_array[y][i].x] += alpha * m_array[y][i].v; - } - } - - void clearRow(uint32_t y) - { - XA_DEBUG_ASSERT( y < height() ); - m_array[y].clear(); - } - - const Array &getRow(uint32_t y) const { return m_array[y]; } - -private: - /// Number of columns. - const uint32_t m_width; - - /// Array of matrix elements. - Array< Array > m_array; -}; - -// y = a * x + y -static void saxpy(float a, const FullVector &x, FullVector &y) -{ - XA_DEBUG_ASSERT(x.dimension() == y.dimension()); - const uint32_t dim = x.dimension(); - for (uint32_t i = 0; i < dim; i++) { - y[i] += a * x[i]; - } -} - -static void copy(const FullVector &x, FullVector &y) -{ - XA_DEBUG_ASSERT(x.dimension() == y.dimension()); - const uint32_t dim = x.dimension(); - for (uint32_t i = 0; i < dim; i++) { - y[i] = x[i]; - } -} - -static void scal(float a, FullVector &x) -{ - const uint32_t dim = x.dimension(); - for (uint32_t i = 0; i < dim; i++) { - x[i] *= a; - } -} - -static float dot(const FullVector &x, const FullVector &y) -{ - XA_DEBUG_ASSERT(x.dimension() == y.dimension()); - const uint32_t dim = x.dimension(); - float sum = 0; - for (uint32_t i = 0; i < dim; i++) { - sum += x[i] * y[i]; - } - return sum; -} - -// y = M * x -static void mult(const Matrix &M, const FullVector &x, FullVector &y) -{ - uint32_t w = M.width(); - uint32_t h = M.height(); - XA_DEBUG_ASSERT( w == x.dimension() ); - XA_UNUSED(w); - XA_DEBUG_ASSERT( h == y.dimension() ); - for (uint32_t i = 0; i < h; i++) - y[i] = M.dotRow(i, x); -} - -// y = alpha*A*x + beta*y -static void sgemv(float alpha, const Matrix &A, const FullVector &x, float beta, FullVector &y) -{ - const uint32_t w = A.width(); - const uint32_t h = A.height(); - XA_DEBUG_ASSERT( w == x.dimension() ); - XA_DEBUG_ASSERT( h == y.dimension() ); - XA_UNUSED(w); - XA_UNUSED(h); - for (uint32_t i = 0; i < h; i++) - y[i] = alpha * A.dotRow(i, x) + beta * y[i]; -} - -// dot y-row of A by x-column of B -static float dotRowColumn(int y, const Matrix &A, int x, const Matrix &B) -{ - const Array &row = A.getRow(y); - const uint32_t count = row.size(); - float sum = 0.0f; - for (uint32_t i = 0; i < count; i++) { - const Matrix::Coefficient &c = row[i]; - sum += c.v * B.getCoefficient(x, c.x); - } - return sum; -} - -static void transpose(const Matrix &A, Matrix &B) -{ - XA_DEBUG_ASSERT(A.width() == B.height()); - XA_DEBUG_ASSERT(B.width() == A.height()); - const uint32_t w = A.width(); - for (uint32_t x = 0; x < w; x++) { - B.clearRow(x); - } - const uint32_t h = A.height(); - for (uint32_t y = 0; y < h; y++) { - const Array &row = A.getRow(y); - const uint32_t count = row.size(); - for (uint32_t i = 0; i < count; i++) { - const Matrix::Coefficient &c = row[i]; - XA_DEBUG_ASSERT(c.x < w); - B.setCoefficient(y, c.x, c.v); - } + return true; } -} -static void sgemm(float alpha, const Matrix &A, const Matrix &B, float beta, Matrix &C) -{ - const uint32_t w = C.width(); - const uint32_t h = C.height(); -#if XA_DEBUG - const uint32_t aw = A.width(); - const uint32_t ah = A.height(); - const uint32_t bw = B.width(); - const uint32_t bh = B.height(); - XA_DEBUG_ASSERT(aw == bh); - XA_DEBUG_ASSERT(bw == ah); - XA_DEBUG_ASSERT(w == bw); - XA_DEBUG_ASSERT(h == ah); -#endif - for (uint32_t y = 0; y < h; y++) { - for (uint32_t x = 0; x < w; x++) { - float c = beta * C.getCoefficient(x, y); - // dot y-row of A by x-column of B. - c += alpha * dotRowColumn(y, A, x, B); - C.setCoefficient(x, y, c); +private: + void flipBackface() + { + // check if triangle is backfacing, if so, swap two vertices + if ( ((v3.x - v1.x) * (v2.y - v1.y) - (v3.y - v1.y) * (v2.x - v1.x)) < 0 ) { + Vector2 hv = v1; + v1 = v2; + v2 = hv; // swap pos } } -} -// C = A * B -static void mult(const Matrix &A, const Matrix &B, Matrix &C) + // compute unit inward normals for each edge. + void computeUnitInwardNormals() + { + n1 = v1 - v2; + n1 = Vector2(-n1.y, n1.x); + n1 = n1 * (1.0f / sqrtf(dot(n1, n1))); + n2 = v2 - v3; + n2 = Vector2(-n2.y, n2.x); + n2 = n2 * (1.0f / sqrtf(dot(n2, n2))); + n3 = v3 - v1; + n3 = Vector2(-n3.y, n3.x); + n3 = n3 * (1.0f / sqrtf(dot(n3, n3))); + } + + // Vertices. + Vector2 v1, v2, v3; + Vector2 n1, n2, n3; // unit inward normals +}; + +// Process the given triangle. Returns false if rasterization was interrupted by the callback. +static bool drawTriangle(const Vector2 &extents, const Vector2 v[3], SamplingCallback cb, void *param) { - sgemm(1.0f, A, B, 0.0f, C); + Triangle tri(v[0], v[1], v[2]); + // @@ It would be nice to have a conservative drawing mode that enlarges the triangle extents by one texel and is able to handle degenerate triangles. + // @@ Maybe the simplest thing to do would be raster triangle edges. + if (tri.isValid()) + return tri.drawAA(extents, cb, param); + return true; } -} // namespace sparse +} // namespace raster namespace segment { @@ -4866,88 +5494,92 @@ private: Array m_pairs; }; -struct Chart +struct AtlasData { - Chart() : faces(MemTag::SegmentAtlasChartFaces) {} - - int id = -1; - Basis basis; // Best fit normal. - float area = 0.0f; - float boundaryLength = 0.0f; - Vector3 centroidSum = Vector3(0.0f); // Sum of chart face centroids. - Vector3 centroid = Vector3(0.0f); // Average centroid of chart faces. - Array seeds; - Array faces; - Array failedPlanarRegions; - CostQueue candidates; -}; + ChartOptions options; + const Mesh *mesh = nullptr; + Array edgeDihedralAngles; + Array edgeLengths; + Array faceAreas; + Array faceNormals; + BitArray isFaceInChart; -struct Atlas -{ - Atlas() : m_edgeLengths(MemTag::SegmentAtlasMeshData), m_faceAreas(MemTag::SegmentAtlasMeshData), m_faceNormals(MemTag::SegmentAtlasMeshData), m_texcoords(MemTag::SegmentAtlasMeshData), m_bestTriangles(10), m_nextPlanarRegionFace(MemTag::SegmentAtlasPlanarRegions), m_facePlanarRegionId(MemTag::SegmentAtlasPlanarRegions) {} + AtlasData() : edgeDihedralAngles(MemTag::SegmentAtlasMeshData), edgeLengths(MemTag::SegmentAtlasMeshData), faceAreas(MemTag::SegmentAtlasMeshData), faceNormals(MemTag::SegmentAtlasMeshData) {} - ~Atlas() + void compute() { - const uint32_t chartCount = m_charts.size(); - for (uint32_t i = 0; i < chartCount; i++) { - m_charts[i]->~Chart(); - XA_FREE(m_charts[i]); + const uint32_t faceCount = mesh->faceCount(); + const uint32_t edgeCount = mesh->edgeCount(); + edgeDihedralAngles.resize(edgeCount); + edgeLengths.resize(edgeCount); + faceAreas.resize(faceCount); + faceNormals.resize(faceCount); + isFaceInChart.resize(faceCount); + isFaceInChart.zeroOutMemory(); + for (uint32_t f = 0; f < faceCount; f++) { + for (uint32_t i = 0; i < 3; i++) { + const uint32_t edge = f * 3 + i; + const Vector3 &p0 = mesh->position(mesh->vertexAt(meshEdgeIndex0(edge))); + const Vector3 &p1 = mesh->position(mesh->vertexAt(meshEdgeIndex1(edge))); + edgeLengths[edge] = length(p1 - p0); + XA_DEBUG_ASSERT(edgeLengths[edge] > 0.0f); + } + faceAreas[f] = mesh->computeFaceArea(f); + XA_DEBUG_ASSERT(faceAreas[f] > 0.0f); + faceNormals[f] = mesh->computeFaceNormal(f); + } + for (uint32_t face = 0; face < faceCount; face++) { + for (uint32_t i = 0; i < 3; i++) { + const uint32_t edge = face * 3 + i; + const uint32_t oedge = mesh->oppositeEdge(edge); + if (oedge == UINT32_MAX) + edgeDihedralAngles[edge] = FLT_MAX; + else { + const uint32_t oface = meshEdgeFace(oedge); + edgeDihedralAngles[edge] = edgeDihedralAngles[oedge] = dot(faceNormals[face], faceNormals[oface]); + } + } } } +}; - uint32_t facesLeft() const { return m_facesLeft; } +#if XA_DEBUG_EXPORT_OBJ_PLANAR_REGIONS +static uint32_t s_planarRegionsCurrentRegion; +static uint32_t s_planarRegionsCurrentVertex; +#endif + +struct PlanarCharts +{ + PlanarCharts(AtlasData &data) : m_data(data), m_nextRegionFace(MemTag::SegmentAtlasPlanarRegions), m_faceToRegionId(MemTag::SegmentAtlasPlanarRegions) {} + const Basis &chartBasis(uint32_t chartIndex) const { return m_chartBasis[chartIndex]; } uint32_t chartCount() const { return m_charts.size(); } - const Array &chartFaces(uint32_t i) const { return m_charts[i]->faces; } - const Basis &chartBasis(uint32_t chartIndex) const { return m_charts[chartIndex]->basis; } + + ConstArrayView chartFaces(uint32_t chartIndex) const + { + const Chart &chart = m_charts[chartIndex]; + return ConstArrayView(&m_chartFaces[chart.firstFace], chart.faceCount); + } - void reset(uint32_t meshId, uint32_t chartGroupId, const Mesh *mesh, const ChartOptions &options) + uint32_t regionIdFromFace(uint32_t face) const { return m_faceToRegionId[face]; } + uint32_t nextRegionFace(uint32_t face) const { return m_nextRegionFace[face]; } + float regionArea(uint32_t region) const { return m_regionAreas[region]; } + + void compute() { - XA_UNUSED(meshId); - XA_UNUSED(chartGroupId); - XA_PROFILE_START(buildAtlasInit) - m_mesh = mesh; - const uint32_t faceCount = m_mesh->faceCount(); - m_facesLeft = faceCount; - m_options = options; - m_rand.reset(); - const uint32_t chartCount = m_charts.size(); - for (uint32_t i = 0; i < chartCount; i++) { - m_charts[i]->~Chart(); - XA_FREE(m_charts[i]); - } - m_charts.clear(); - m_faceCharts.resize(faceCount); - m_faceCharts.setAll(-1); - m_texcoords.resize(faceCount * 3); - // Precompute edge lengths and face areas. - const uint32_t edgeCount = m_mesh->edgeCount(); - m_edgeLengths.resize(edgeCount); - m_faceAreas.resize(faceCount); - m_faceNormals.resize(faceCount); - for (uint32_t f = 0; f < faceCount; f++) { - for (uint32_t i = 0; i < 3; i++) { - const uint32_t edge = f * 3 + i; - const Vector3 &p0 = mesh->position(m_mesh->vertexAt(meshEdgeIndex0(edge))); - const Vector3 &p1 = mesh->position(m_mesh->vertexAt(meshEdgeIndex1(edge))); - m_edgeLengths[edge] = length(p1 - p0); - XA_DEBUG_ASSERT(m_edgeLengths[edge] > 0.0f); - } - m_faceAreas[f] = m_mesh->computeFaceArea(f); - XA_DEBUG_ASSERT(m_faceAreas[f] > 0.0f); - m_faceNormals[f] = m_mesh->computeFaceNormal(f); - } + const uint32_t faceCount = m_data.mesh->faceCount(); // Precompute regions of coplanar incident faces. - m_nextPlanarRegionFace.resize(faceCount); - m_facePlanarRegionId.resize(faceCount); + m_regionFirstFace.clear(); + m_nextRegionFace.resize(faceCount); + m_faceToRegionId.resize(faceCount); for (uint32_t f = 0; f < faceCount; f++) { - m_nextPlanarRegionFace[f] = f; - m_facePlanarRegionId[f] = UINT32_MAX; + m_nextRegionFace[f] = f; + m_faceToRegionId[f] = UINT32_MAX; } Array faceStack; faceStack.reserve(min(faceCount, 16u)); - uint32_t planarRegionCount = 0; + uint32_t regionCount = 0; for (uint32_t f = 0; f < faceCount; f++) { - if (m_nextPlanarRegionFace[f] != f) + if (m_nextRegionFace[f] != f) continue; // Already assigned. faceStack.clear(); faceStack.push_back(f); @@ -4955,49 +5587,207 @@ struct Atlas if (faceStack.isEmpty()) break; const uint32_t face = faceStack.back(); - m_facePlanarRegionId[face] = planarRegionCount; + m_faceToRegionId[face] = regionCount; faceStack.pop_back(); - for (Mesh::FaceEdgeIterator it(m_mesh, face); !it.isDone(); it.advance()) { + for (Mesh::FaceEdgeIterator it(m_data.mesh, face); !it.isDone(); it.advance()) { const uint32_t oface = it.oppositeFace(); if (it.isBoundary()) continue; - if (m_nextPlanarRegionFace[oface] != oface) + if (m_nextRegionFace[oface] != oface) continue; // Already assigned. - if (!equal(dot(m_faceNormals[face], m_faceNormals[oface]), 1.0f, kEpsilon)) + if (!equal(dot(m_data.faceNormals[face], m_data.faceNormals[oface]), 1.0f, kEpsilon)) continue; // Not coplanar. - const uint32_t next = m_nextPlanarRegionFace[face]; - m_nextPlanarRegionFace[face] = oface; - m_nextPlanarRegionFace[oface] = next; - m_facePlanarRegionId[oface] = planarRegionCount; + const uint32_t next = m_nextRegionFace[face]; + m_nextRegionFace[face] = oface; + m_nextRegionFace[oface] = next; + m_faceToRegionId[oface] = regionCount; faceStack.push_back(oface); } } - planarRegionCount++; + m_regionFirstFace.push_back(f); + regionCount++; } #if XA_DEBUG_EXPORT_OBJ_PLANAR_REGIONS - char filename[256]; - XA_SPRINTF(filename, sizeof(filename), "debug_mesh_%03u_chartgroup_%03u_planar_regions.obj", meshId, chartGroupId); - FILE *file; - XA_FOPEN(file, filename, "w"); - if (file) { - m_mesh->writeObjVertices(file); - fprintf(file, "s off\n"); - for (uint32_t i = 0; i < planarRegionCount; i++) { - fprintf(file, "o region%u\n", i); - for (uint32_t j = 0; j < faceCount; j++) { - if (m_facePlanarRegionId[j] == i) - m_mesh->writeObjFace(file, j); + static std::mutex s_mutex; + { + std::lock_guard lock(s_mutex); + FILE *file; + XA_FOPEN(file, "debug_mesh_planar_regions.obj", s_planarRegionsCurrentRegion == 0 ? "w" : "a"); + if (file) { + m_data.mesh->writeObjVertices(file); + fprintf(file, "s off\n"); + for (uint32_t i = 0; i < regionCount; i++) { + fprintf(file, "o region%u\n", s_planarRegionsCurrentRegion); + for (uint32_t j = 0; j < faceCount; j++) { + if (m_faceToRegionId[j] == i) + m_data.mesh->writeObjFace(file, j, s_planarRegionsCurrentVertex); + } + s_planarRegionsCurrentRegion++; } + s_planarRegionsCurrentVertex += m_data.mesh->vertexCount(); + fclose(file); } - fclose(file); } #endif - XA_PROFILE_END(buildAtlasInit) + // Precompute planar region areas. + m_regionAreas.resize(regionCount); + m_regionAreas.zeroOutMemory(); + for (uint32_t f = 0; f < faceCount; f++) + m_regionAreas[m_faceToRegionId[f]] += m_data.faceAreas[f]; + // Create charts from suitable planar regions. + // The dihedral angle of all boundary edges must be >= 90 degrees. + m_charts.clear(); + m_chartFaces.clear(); + for (uint32_t region = 0; region < regionCount; region++) { + const uint32_t firstRegionFace = m_regionFirstFace[region]; + uint32_t face = firstRegionFace; + bool createChart = true; + do { + for (Mesh::FaceEdgeIterator it(m_data.mesh, face); !it.isDone(); it.advance()) { + if (it.isBoundary()) + continue; // Ignore mesh boundary edges. + const uint32_t oface = it.oppositeFace(); + if (m_faceToRegionId[oface] == region) + continue; // Ignore internal edges. + const float angle = m_data.edgeDihedralAngles[it.edge()]; + if (angle > 0.0f && angle < FLT_MAX) { // FLT_MAX on boundaries. + createChart = false; + break; + } + } + if (!createChart) + break; + face = m_nextRegionFace[face]; + } + while (face != firstRegionFace); + // Create a chart. + if (createChart) { + Chart chart; + chart.firstFace = m_chartFaces.size(); + chart.faceCount = 0; + face = firstRegionFace; + do { + m_data.isFaceInChart.set(face); + m_chartFaces.push_back(face); + chart.faceCount++; + face = m_nextRegionFace[face]; + } + while (face != firstRegionFace); + m_charts.push_back(chart); + } + } + // Compute basis for each chart using the first face normal (all faces have the same normal). + m_chartBasis.resize(m_charts.size()); + for (uint32_t c = 0; c < m_charts.size(); c++) + { + const uint32_t face = m_chartFaces[m_charts[c].firstFace]; + Basis &basis = m_chartBasis[c]; + basis.normal = m_data.faceNormals[face]; + basis.tangent = Basis::computeTangent(basis.normal); + basis.bitangent = Basis::computeBitangent(basis.normal, basis.tangent); + } + } + +private: + struct Chart + { + uint32_t firstFace, faceCount; + }; + + AtlasData &m_data; + Array m_regionFirstFace; + Array m_nextRegionFace; + Array m_faceToRegionId; + Array m_regionAreas; + Array m_charts; + Array m_chartFaces; + Array m_chartBasis; +}; + +struct ClusteredCharts +{ + ClusteredCharts(AtlasData &data, const PlanarCharts &planarCharts) : m_data(data), m_planarCharts(planarCharts), m_texcoords(MemTag::SegmentAtlasMeshData), m_bestTriangles(10), m_placingSeeds(false) {} + + ~ClusteredCharts() + { + const uint32_t chartCount = m_charts.size(); + for (uint32_t i = 0; i < chartCount; i++) { + m_charts[i]->~Chart(); + XA_FREE(m_charts[i]); + } + } + + uint32_t chartCount() const { return m_charts.size(); } + ConstArrayView chartFaces(uint32_t chartIndex) const { return m_charts[chartIndex]->faces; } + const Basis &chartBasis(uint32_t chartIndex) const { return m_charts[chartIndex]->basis; } + + void compute() + { + const uint32_t faceCount = m_data.mesh->faceCount(); + m_facesLeft = 0; + for (uint32_t i = 0; i < faceCount; i++) { + if (!m_data.isFaceInChart.get(i)) + m_facesLeft++; + } + const uint32_t chartCount = m_charts.size(); + for (uint32_t i = 0; i < chartCount; i++) { + m_charts[i]->~Chart(); + XA_FREE(m_charts[i]); + } + m_charts.clear(); + m_faceCharts.resize(faceCount); + m_faceCharts.fill(-1); + m_texcoords.resize(faceCount * 3); + if (m_facesLeft == 0) + return; + // Create initial charts greedely. + placeSeeds(m_data.options.maxCost * 0.5f); + if (m_data.options.maxIterations == 0) { + XA_DEBUG_ASSERT(m_facesLeft == 0); + return; + } + relocateSeeds(); + resetCharts(); + // Restart process growing charts in parallel. + uint32_t iteration = 0; + for (;;) { + growCharts(m_data.options.maxCost); + // When charts cannot grow more: fill holes, merge charts, relocate seeds and start new iteration. + fillHoles(m_data.options.maxCost * 0.5f); +#if XA_MERGE_CHARTS + mergeCharts(); +#endif + if (++iteration == m_data.options.maxIterations) + break; + if (!relocateSeeds()) + break; + resetCharts(); + } + // Make sure no holes are left! + XA_DEBUG_ASSERT(m_facesLeft == 0); } +private: + struct Chart + { + Chart() : faces(MemTag::SegmentAtlasChartFaces) {} + + int id = -1; + Basis basis; // Best fit normal. + float area = 0.0f; + float boundaryLength = 0.0f; + Vector3 centroidSum = Vector3(0.0f); // Sum of chart face centroids. + Vector3 centroid = Vector3(0.0f); // Average centroid of chart faces. + Array faces; + Array failedPlanarRegions; + CostQueue candidates; + uint32_t seed; + }; + void placeSeeds(float threshold) { - XA_PROFILE_START(buildAtlasPlaceSeeds) + XA_PROFILE_START(clusteredChartsPlaceSeeds) + m_placingSeeds = true; // 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. @@ -5005,14 +5795,15 @@ struct Atlas // - those points can be found using a simple flood filling algorithm. // - how do we weight the probabilities? while (m_facesLeft > 0) - createRandomChart(threshold); - XA_PROFILE_END(buildAtlasPlaceSeeds) + createChart(threshold); + m_placingSeeds = false; + XA_PROFILE_END(clusteredChartsPlaceSeeds) } // Returns true if any of the charts can grow more. void growCharts(float threshold) { - XA_PROFILE_START(buildAtlasGrowCharts) + XA_PROFILE_START(clusteredChartsGrow) for (;;) { if (m_facesLeft == 0) break; @@ -5030,7 +5821,7 @@ struct Atlas break; cost = chart->candidates.peekCost(); face = chart->candidates.peekFace(); - if (m_faceCharts[face] == -1) + if (!m_data.isFaceInChart.get(face)) break; else { // Face belongs to another chart. Pop from queue so the next best candidate can be retrieved. @@ -5052,22 +5843,28 @@ struct Atlas Chart *chart = m_charts[bestChart]; chart->candidates.pop(); // Pop the selected candidate from the queue. if (!addFaceToChart(chart, bestFace)) - chart->failedPlanarRegions.push_back(m_facePlanarRegionId[bestFace]); + chart->failedPlanarRegions.push_back(m_planarCharts.regionIdFromFace(bestFace)); } - XA_PROFILE_END(buildAtlasGrowCharts) + XA_PROFILE_END(clusteredChartsGrow) } void resetCharts() { - XA_PROFILE_START(buildAtlasResetCharts) - const uint32_t faceCount = m_mesh->faceCount(); - for (uint32_t i = 0; i < faceCount; i++) + XA_PROFILE_START(clusteredChartsReset) + const uint32_t faceCount = m_data.mesh->faceCount(); + for (uint32_t i = 0; i < faceCount; i++) { + if (m_faceCharts[i] != -1) + m_data.isFaceInChart.unset(i); m_faceCharts[i] = -1; - m_facesLeft = faceCount; + } + m_facesLeft = 0; + for (uint32_t i = 0; i < faceCount; i++) { + if (!m_data.isFaceInChart.get(i)) + m_facesLeft++; + } const uint32_t chartCount = m_charts.size(); for (uint32_t i = 0; i < chartCount; i++) { Chart *chart = m_charts[i]; - const uint32_t seed = chart->seeds.back(); chart->area = 0.0f; chart->boundaryLength = 0.0f; chart->basis.normal = Vector3(0.0f); @@ -5078,14 +5875,14 @@ struct Atlas chart->faces.clear(); chart->candidates.clear(); chart->failedPlanarRegions.clear(); - addFaceToChart(chart, seed); + addFaceToChart(chart, chart->seed); } - XA_PROFILE_END(buildAtlasResetCharts) + XA_PROFILE_END(clusteredChartsReset) } bool relocateSeeds() { - XA_PROFILE_START(buildAtlasRelocateSeeds) + XA_PROFILE_START(clusteredChartsRelocateSeeds) bool anySeedChanged = false; const uint32_t chartCount = m_charts.size(); for (uint32_t i = 0; i < chartCount; i++) { @@ -5093,22 +5890,22 @@ struct Atlas anySeedChanged = true; } } - XA_PROFILE_END(buildAtlasRelocateSeeds) + XA_PROFILE_END(clusteredChartsRelocateSeeds) return anySeedChanged; } void fillHoles(float threshold) { - XA_PROFILE_START(buildAtlasFillHoles) + XA_PROFILE_START(clusteredChartsFillHoles) while (m_facesLeft > 0) - createRandomChart(threshold); - XA_PROFILE_END(buildAtlasFillHoles) + createChart(threshold); + XA_PROFILE_END(clusteredChartsFillHoles) } #if XA_MERGE_CHARTS void mergeCharts() { - XA_PROFILE_START(buildAtlasMergeCharts) + XA_PROFILE_START(clusteredChartsMerge) const uint32_t chartCount = m_charts.size(); // Merge charts progressively until there's none left to merge. for (;;) { @@ -5127,13 +5924,15 @@ struct Atlas const uint32_t faceCount = chart->faces.size(); for (uint32_t i = 0; i < faceCount; i++) { const uint32_t f = chart->faces[i]; - for (Mesh::FaceEdgeIterator it(m_mesh, f); !it.isDone(); it.advance()) { - const float l = m_edgeLengths[it.edge()]; + for (Mesh::FaceEdgeIterator it(m_data.mesh, f); !it.isDone(); it.advance()) { + const float l = m_data.edgeLengths[it.edge()]; if (it.isBoundary()) { externalBoundaryLength += l; } else { const int neighborChart = m_faceCharts[it.oppositeFace()]; - if (m_charts[neighborChart] != chart) { + if (neighborChart == -1) + externalBoundaryLength += l; + else if (m_charts[neighborChart] != chart) { if ((it.isSeam() && (isNormalSeam(it.edge()) || it.isTextureSeam()))) { externalBoundaryLength += l; } else { @@ -5158,9 +5957,9 @@ struct Atlas if (dot(chart2->basis.normal, chart->basis.normal) < XA_MERGE_CHARTS_MIN_NORMAL_DEVIATION) continue; // Obey max chart area and boundary length. - if (m_options.maxChartArea > 0.0f && chart->area + chart2->area > m_options.maxChartArea) + if (m_data.options.maxChartArea > 0.0f && chart->area + chart2->area > m_data.options.maxChartArea) continue; - if (m_options.maxBoundaryLength > 0.0f && chart->boundaryLength + chart2->boundaryLength - m_sharedBoundaryLengthsNoSeams[cc] > m_options.maxBoundaryLength) + if (m_data.options.maxBoundaryLength > 0.0f && chart->boundaryLength + chart2->boundaryLength - m_sharedBoundaryLengthsNoSeams[cc] > m_data.options.maxBoundaryLength) continue; // Merge if chart2 has a single face. // chart1 must have more than 1 face. @@ -5207,33 +6006,38 @@ struct Atlas c++; } } - XA_PROFILE_END(buildAtlasMergeCharts) + XA_PROFILE_END(clusteredChartsMerge) } #endif private: - void createRandomChart(float threshold) + void createChart(float threshold) { Chart *chart = XA_NEW(MemTag::Default, Chart); chart->id = (int)m_charts.size(); m_charts.push_back(chart); - // Pick random face that is not used by any chart yet. - uint32_t face = m_rand.getRange(m_mesh->faceCount() - 1); - while (m_faceCharts[face] != -1) { - if (++face >= m_mesh->faceCount()) - face = 0; - } - chart->seeds.push_back(face); - addFaceToChart(chart, face); + // Pick a face not used by any chart yet, belonging to the largest planar region. + chart->seed = 0; + float largestArea = 0.0f; + for (uint32_t f = 0; f < m_data.mesh->faceCount(); f++) { + if (m_data.isFaceInChart.get(f)) + continue; + const float area = m_planarCharts.regionArea(m_planarCharts.regionIdFromFace(f)); + if (area > largestArea) { + largestArea = area; + chart->seed = f; + } + } + addFaceToChart(chart, chart->seed); // Grow the chart as much as possible within the given threshold. for (;;) { if (chart->candidates.count() == 0 || chart->candidates.peekCost() > threshold) break; const uint32_t f = chart->candidates.pop(); - if (m_faceCharts[f] != -1) + if (m_data.isFaceInChart.get(f)) continue; if (!addFaceToChart(chart, f)) { - chart->failedPlanarRegions.push_back(m_facePlanarRegionId[f]); + chart->failedPlanarRegions.push_back(m_planarCharts.regionIdFromFace(f)); continue; } } @@ -5241,7 +6045,7 @@ private: bool isChartBoundaryEdge(const Chart *chart, uint32_t edge) const { - const uint32_t oppositeEdge = m_mesh->oppositeEdge(edge); + const uint32_t oppositeEdge = m_data.mesh->oppositeEdge(edge); const uint32_t oppositeFace = meshEdgeFace(oppositeEdge); return oppositeEdge == UINT32_MAX || m_faceCharts[oppositeFace] != chart->id; } @@ -5253,7 +6057,7 @@ private: for (uint32_t i = 0; i < faceCount; i++) { const uint32_t f = chart->faces[i]; for (uint32_t j = 0; j < 3; j++) - m_tempPoints[i * 3 + j] = m_mesh->position(m_mesh->vertexAt(f * 3 + j)); + m_tempPoints[i * 3 + j] = m_data.mesh->position(m_data.mesh->vertexAt(f * 3 + j)); } return Fit::computeBasis(m_tempPoints.data(), m_tempPoints.size(), basis); } @@ -5274,7 +6078,7 @@ private: const uint32_t face = chart->faces[i]; for (uint32_t j = 0; j < 3; j++) { const uint32_t offset = face * 3 + j; - const Vector3 &pos = m_mesh->position(m_mesh->vertexAt(offset)); + const Vector3 &pos = m_data.mesh->position(m_data.mesh->vertexAt(offset)); m_texcoords[offset] = Vector2(dot(chart->basis.tangent, pos), dot(chart->basis.bitangent, pos)); } } @@ -5293,6 +6097,8 @@ private: if (flippedFaceCount != 0 && flippedFaceCount != faceCount) return false; // Check for boundary intersection in the parameterization. + XA_PROFILE_START(clusteredChartsPlaceSeedsBoundaryIntersection) + XA_PROFILE_START(clusteredChartsGrowBoundaryIntersection) m_boundaryGrid.reset(m_texcoords.data()); for (uint32_t i = 0; i < faceCount; i++) { const uint32_t f = chart->faces[i]; @@ -5302,23 +6108,30 @@ private: m_boundaryGrid.append(edge); } } - if (m_boundaryGrid.intersectSelf(m_mesh->epsilon())) + const bool intersection = m_boundaryGrid.intersect(m_data.mesh->epsilon()); +#if XA_PROFILE + if (m_placingSeeds) + XA_PROFILE_END(clusteredChartsPlaceSeedsBoundaryIntersection) + else + XA_PROFILE_END(clusteredChartsGrowBoundaryIntersection) +#endif + if (intersection) return false; return true; } bool addFaceToChart(Chart *chart, uint32_t face) { - XA_DEBUG_ASSERT(m_faceCharts[face] == -1); + XA_DEBUG_ASSERT(!m_data.isFaceInChart.get(face)); const uint32_t oldFaceCount = chart->faces.size(); const bool firstFace = oldFaceCount == 0; // Append the face and any coplanar connected faces to the chart faces array. chart->faces.push_back(face); - uint32_t coplanarFace = m_nextPlanarRegionFace[face]; + uint32_t coplanarFace = m_planarCharts.nextRegionFace(face); while (coplanarFace != face) { - XA_DEBUG_ASSERT(m_faceCharts[coplanarFace] == -1); + XA_DEBUG_ASSERT(!m_data.isFaceInChart.get(coplanarFace)); chart->faces.push_back(coplanarFace); - coplanarFace = m_nextPlanarRegionFace[coplanarFace]; + coplanarFace = m_planarCharts.nextRegionFace(coplanarFace); } const uint32_t faceCount = chart->faces.size(); // Compute basis. @@ -5326,8 +6139,8 @@ private: if (firstFace) { // Use the first face normal. // Use any edge as the tangent vector. - basis.normal = m_faceNormals[face]; - basis.tangent = normalize(m_mesh->position(m_mesh->vertexAt(face * 3 + 0)) - m_mesh->position(m_mesh->vertexAt(face * 3 + 1)), kEpsilon); + basis.normal = m_data.faceNormals[face]; + basis.tangent = normalize(m_data.mesh->position(m_data.mesh->vertexAt(face * 3 + 0)) - m_data.mesh->position(m_data.mesh->vertexAt(face * 3 + 1)), kEpsilon); basis.bitangent = cross(basis.normal, basis.tangent); } else { // Use best fit normal. @@ -5335,7 +6148,7 @@ private: chart->faces.resize(oldFaceCount); return false; } - if (dot(basis.normal, m_faceNormals[face]) < 0.0f) // Flip normal if oriented in the wrong direction. + if (dot(basis.normal, m_data.faceNormals[face]) < 0.0f) // Flip normal if oriented in the wrong direction. basis.normal = -basis.normal; } if (!firstFace) { @@ -5358,7 +6171,8 @@ private: const uint32_t f = chart->faces[i]; m_faceCharts[f] = chart->id; m_facesLeft--; - chart->centroidSum += m_mesh->computeFaceCenter(f); + m_data.isFaceInChart.set(f); + chart->centroidSum += m_data.mesh->computeFaceCenter(f); } chart->centroid = chart->centroidSum / float(chart->faces.size()); // Refresh candidates. @@ -5368,15 +6182,15 @@ private: const uint32_t f = chart->faces[i]; for (uint32_t j = 0; j < 3; j++) { const uint32_t edge = f * 3 + j; - const uint32_t oedge = m_mesh->oppositeEdge(edge); + const uint32_t oedge = m_data.mesh->oppositeEdge(edge); if (oedge == UINT32_MAX) continue; // Boundary edge. const uint32_t oface = meshEdgeFace(oedge); - if (m_faceCharts[oface] != -1) + if (m_data.isFaceInChart.get(oface)) continue; // Face belongs to another chart. - if (chart->failedPlanarRegions.contains(m_facePlanarRegionId[oface])) + if (chart->failedPlanarRegions.contains(m_planarCharts.regionIdFromFace(oface))) continue; // Failed to add this faces planar region to the chart before. - const float cost = evaluateCost(chart, oface); + const float cost = computeCost(chart, oface); if (cost < FLT_MAX) chart->candidates.push(cost, oface); } @@ -5391,72 +6205,56 @@ private: const uint32_t faceCount = chart->faces.size(); m_bestTriangles.clear(); for (uint32_t i = 0; i < faceCount; i++) { - const float cost = evaluateProxyFitMetric(chart, chart->faces[i]); + const float cost = computeNormalDeviationMetric(chart, chart->faces[i]); m_bestTriangles.push(cost, chart->faces[i]); } - // Of those, choose the least central triangle. - uint32_t leastCentral = 0; - float maxDistance = -1; + // Of those, choose the most central triangle. + uint32_t mostCentral = 0; + float minDistance = FLT_MAX; for (;;) { if (m_bestTriangles.count() == 0) break; const uint32_t face = m_bestTriangles.pop(); - Vector3 faceCentroid = m_mesh->computeFaceCenter(face); + Vector3 faceCentroid = m_data.mesh->computeFaceCenter(face); const float distance = length(chart->centroid - faceCentroid); - if (distance > maxDistance) { - maxDistance = distance; - leastCentral = face; - } - } - XA_DEBUG_ASSERT(maxDistance >= 0); - // In order to prevent k-means cyles we record all the previously chosen seeds. - for (uint32_t i = 0; i < chart->seeds.size(); i++) { - // Treat seeds belong to the same planar region as equal. - if (chart->seeds[i] == leastCentral || m_facePlanarRegionId[chart->seeds[i]] == m_facePlanarRegionId[leastCentral]) { - // Move new seed to the end of the seed array. - uint32_t last = chart->seeds.size() - 1; - swap(chart->seeds[i], chart->seeds[last]); - return false; + if (distance < minDistance) { + minDistance = distance; + mostCentral = face; } } - // Append new seed. - chart->seeds.push_back(leastCentral); + XA_DEBUG_ASSERT(minDistance < FLT_MAX); + if (mostCentral == chart->seed) + return false; + chart->seed = mostCentral; return true; } - // Evaluate combined metric. - float evaluateCost(Chart *chart, uint32_t face) const + // Cost is combined metrics * weights. + float computeCost(Chart *chart, uint32_t face) const { - if (dot(m_faceNormals[face], chart->basis.normal) <= 0.26f) // ~75 degrees - return FLT_MAX; // Estimate boundary length and area: - float newChartArea = 0.0f, newBoundaryLength = 0.0f; - if (m_options.maxChartArea > 0.0f || m_options.roundnessMetricWeight > 0.0f) - newChartArea = computeArea(chart, face); - if (m_options.maxBoundaryLength > 0.0f || m_options.roundnessMetricWeight > 0.0f) - newBoundaryLength = computeBoundaryLength(chart, face); + const float newChartArea = computeArea(chart, face); + const float newBoundaryLength = computeBoundaryLength(chart, face); // Enforce limits strictly: - if (m_options.maxChartArea > 0.0f && newChartArea > m_options.maxChartArea) + if (m_data.options.maxChartArea > 0.0f && newChartArea > m_data.options.maxChartArea) return FLT_MAX; - if (m_options.maxBoundaryLength > 0.0f && newBoundaryLength > m_options.maxBoundaryLength) + if (m_data.options.maxBoundaryLength > 0.0f && newBoundaryLength > m_data.options.maxBoundaryLength) return FLT_MAX; + // Compute metrics. float cost = 0.0f; - if (m_options.normalSeamMetricWeight > 0.0f) { - // Penalize faces that cross seams, reward faces that close seams or reach boundaries. - // Make sure normal seams are fully respected: - const float N = evaluateNormalSeamMetric(chart, face); - if (m_options.normalSeamMetricWeight >= 1000.0f && N > 0.0f) - return FLT_MAX; - cost += m_options.normalSeamMetricWeight * N; - } - if (m_options.proxyFitMetricWeight > 0.0f) - cost += m_options.proxyFitMetricWeight * evaluateProxyFitMetric(chart, face); - if (m_options.roundnessMetricWeight > 0.0f) - cost += m_options.roundnessMetricWeight * evaluateRoundnessMetric(chart, newBoundaryLength, newChartArea); - if (m_options.straightnessMetricWeight > 0.0f) - cost += m_options.straightnessMetricWeight * evaluateStraightnessMetric(chart, face); - if (m_options.textureSeamMetricWeight > 0.0f) - cost += m_options.textureSeamMetricWeight * evaluateTextureSeamMetric(chart, face); + const float normalDeviation = computeNormalDeviationMetric(chart, face); + if (normalDeviation >= 0.707f) // ~75 degrees + return FLT_MAX; + cost += m_data.options.normalDeviationWeight * normalDeviation; + // Penalize faces that cross seams, reward faces that close seams or reach boundaries. + // Make sure normal seams are fully respected: + const float normalSeam = computeNormalSeamMetric(chart, face); + if (m_data.options.normalSeamWeight >= 1000.0f && normalSeam > 0.0f) + return FLT_MAX; + cost += m_data.options.normalSeamWeight * normalSeam; + cost += m_data.options.roundnessWeight * computeRoundnessMetric(chart, newBoundaryLength, newChartArea); + cost += m_data.options.straightnessWeight * computeStraightnessMetric(chart, face); + cost += m_data.options.textureSeamWeight * computeTextureSeamMetric(chart, face); //float R = evaluateCompletenessMetric(chart, face); //float D = evaluateDihedralAngleMetric(chart, face); // @@ Add a metric based on local dihedral angle. @@ -5467,105 +6265,108 @@ private: } // Returns a value in [0-1]. - float evaluateProxyFitMetric(Chart *chart, uint32_t face) const + // 0 if face normal is coplanar to the chart's best fit normal. + // 1 if face normal is perpendicular. + float computeNormalDeviationMetric(Chart *chart, uint32_t face) const { // All faces in coplanar regions have the same normal, can use any face. - const Vector3 faceNormal = m_faceNormals[face]; + const Vector3 faceNormal = m_data.faceNormals[face]; // Use plane fitting metric for now: - return 1 - dot(faceNormal, chart->basis.normal); // @@ normal deviations should be weighted by face area + return min(1.0f - dot(faceNormal, chart->basis.normal), 1.0f); // @@ normal deviations should be weighted by face area } - float evaluateRoundnessMetric(Chart *chart, float newBoundaryLength, float newChartArea) const + float computeRoundnessMetric(Chart *chart, float newBoundaryLength, float newChartArea) const { - const float roundness = square(chart->boundaryLength) / chart->area; - const float newBoundaryLengthSq = square(newBoundaryLength); - const float newRoundness = newBoundaryLengthSq / newChartArea; - if (newRoundness > roundness) - return newBoundaryLengthSq / (newChartArea * kPi4); - // Offer no impedance to faces that improve roundness. - return 0; + const float oldRoundness = square(chart->boundaryLength) / chart->area; + const float newRoundness = square(newBoundaryLength) / newChartArea; + return 1.0f - oldRoundness / newRoundness; } - float evaluateStraightnessMetric(Chart *chart, uint32_t firstFace) const + float computeStraightnessMetric(Chart *chart, uint32_t firstFace) const { - float l_out = 0.0f, l_in = 0.0f; - const uint32_t planarRegionId = m_facePlanarRegionId[firstFace]; + float l_out = 0.0f; // Length of firstFace planar region boundary that doesn't border the chart. + float l_in = 0.0f; // Length that does border the chart. + const uint32_t planarRegionId = m_planarCharts.regionIdFromFace(firstFace); uint32_t face = firstFace; for (;;) { - for (Mesh::FaceEdgeIterator it(m_mesh, face); !it.isDone(); it.advance()) { - const float l = m_edgeLengths[it.edge()]; + for (Mesh::FaceEdgeIterator it(m_data.mesh, face); !it.isDone(); it.advance()) { + const float l = m_data.edgeLengths[it.edge()]; if (it.isBoundary()) { l_out += l; - } else if (m_facePlanarRegionId[it.oppositeFace()] != planarRegionId) { + } else if (m_planarCharts.regionIdFromFace(it.oppositeFace()) != planarRegionId) { if (m_faceCharts[it.oppositeFace()] != chart->id) l_out += l; else l_in += l; } } - face = m_nextPlanarRegionFace[face]; + face = m_planarCharts.nextRegionFace(face); if (face == firstFace) break; } +#if 1 XA_DEBUG_ASSERT(l_in != 0.0f); // Candidate face must be adjacent to chart. @@ This is not true if the input mesh has zero-length edges. float ratio = (l_out - l_in) / (l_out + l_in); return min(ratio, 0.0f); // Only use the straightness metric to close gaps. +#else + return 1.0f - l_in / l_out; +#endif } bool isNormalSeam(uint32_t edge) const { - const uint32_t oppositeEdge = m_mesh->oppositeEdge(edge); + const uint32_t oppositeEdge = m_data.mesh->oppositeEdge(edge); if (oppositeEdge == UINT32_MAX) return false; // boundary edge - if (m_mesh->flags() & MeshFlags::HasNormals) { - const uint32_t v0 = m_mesh->vertexAt(meshEdgeIndex0(edge)); - const uint32_t v1 = m_mesh->vertexAt(meshEdgeIndex1(edge)); - const uint32_t ov0 = m_mesh->vertexAt(meshEdgeIndex0(oppositeEdge)); - const uint32_t ov1 = m_mesh->vertexAt(meshEdgeIndex1(oppositeEdge)); + if (m_data.mesh->flags() & MeshFlags::HasNormals) { + const uint32_t v0 = m_data.mesh->vertexAt(meshEdgeIndex0(edge)); + const uint32_t v1 = m_data.mesh->vertexAt(meshEdgeIndex1(edge)); + const uint32_t ov0 = m_data.mesh->vertexAt(meshEdgeIndex0(oppositeEdge)); + const uint32_t ov1 = m_data.mesh->vertexAt(meshEdgeIndex1(oppositeEdge)); if (v0 == ov1 && v1 == ov0) return false; - return !equal(m_mesh->normal(v0), m_mesh->normal(ov1), kNormalEpsilon) || !equal(m_mesh->normal(v1), m_mesh->normal(ov0), kNormalEpsilon); + return !equal(m_data.mesh->normal(v0), m_data.mesh->normal(ov1), kNormalEpsilon) || !equal(m_data.mesh->normal(v1), m_data.mesh->normal(ov0), kNormalEpsilon); } const uint32_t f0 = meshEdgeFace(edge); const uint32_t f1 = meshEdgeFace(oppositeEdge); - if (m_facePlanarRegionId[f0] == m_facePlanarRegionId[f1]) + if (m_planarCharts.regionIdFromFace(f0) == m_planarCharts.regionIdFromFace(f1)) return false; - return !equal(m_faceNormals[f0], m_faceNormals[f1], kNormalEpsilon); + return !equal(m_data.faceNormals[f0], m_data.faceNormals[f1], kNormalEpsilon); } - float evaluateNormalSeamMetric(Chart *chart, uint32_t firstFace) const + float computeNormalSeamMetric(Chart *chart, uint32_t firstFace) const { float seamFactor = 0.0f, totalLength = 0.0f; uint32_t face = firstFace; for (;;) { - for (Mesh::FaceEdgeIterator it(m_mesh, face); !it.isDone(); it.advance()) { + for (Mesh::FaceEdgeIterator it(m_data.mesh, face); !it.isDone(); it.advance()) { if (it.isBoundary()) continue; if (m_faceCharts[it.oppositeFace()] != chart->id) continue; - float l = m_edgeLengths[it.edge()]; + float l = m_data.edgeLengths[it.edge()]; totalLength += l; if (!it.isSeam()) continue; // Make sure it's a normal seam. if (isNormalSeam(it.edge())) { float d; - if (m_mesh->flags() & MeshFlags::HasNormals) { - const Vector3 &n0 = m_mesh->normal(it.vertex0()); - const Vector3 &n1 = m_mesh->normal(it.vertex1()); - const Vector3 &on0 = m_mesh->normal(m_mesh->vertexAt(meshEdgeIndex0(it.oppositeEdge()))); - const Vector3 &on1 = m_mesh->normal(m_mesh->vertexAt(meshEdgeIndex1(it.oppositeEdge()))); + if (m_data.mesh->flags() & MeshFlags::HasNormals) { + const Vector3 &n0 = m_data.mesh->normal(it.vertex0()); + const Vector3 &n1 = m_data.mesh->normal(it.vertex1()); + const Vector3 &on0 = m_data.mesh->normal(m_data.mesh->vertexAt(meshEdgeIndex0(it.oppositeEdge()))); + const Vector3 &on1 = m_data.mesh->normal(m_data.mesh->vertexAt(meshEdgeIndex1(it.oppositeEdge()))); const float d0 = clamp(dot(n0, on1), 0.0f, 1.0f); const float d1 = clamp(dot(n1, on0), 0.0f, 1.0f); d = (d0 + d1) * 0.5f; } else { - d = clamp(dot(m_faceNormals[face], m_faceNormals[meshEdgeFace(it.oppositeEdge())]), 0.0f, 1.0f); + d = clamp(dot(m_data.faceNormals[face], m_data.faceNormals[meshEdgeFace(it.oppositeEdge())]), 0.0f, 1.0f); } l *= 1 - d; seamFactor += l; } } - face = m_nextPlanarRegionFace[face]; + face = m_planarCharts.nextRegionFace(face); if (face == firstFace) break; } @@ -5574,17 +6375,17 @@ private: return seamFactor / totalLength; } - float evaluateTextureSeamMetric(Chart *chart, uint32_t firstFace) const + float computeTextureSeamMetric(Chart *chart, uint32_t firstFace) const { float seamLength = 0.0f, totalLength = 0.0f; uint32_t face = firstFace; for (;;) { - for (Mesh::FaceEdgeIterator it(m_mesh, face); !it.isDone(); it.advance()) { + for (Mesh::FaceEdgeIterator it(m_data.mesh, face); !it.isDone(); it.advance()) { if (it.isBoundary()) continue; if (m_faceCharts[it.oppositeFace()] != chart->id) continue; - float l = m_edgeLengths[it.edge()]; + float l = m_data.edgeLengths[it.edge()]; totalLength += l; if (!it.isSeam()) continue; @@ -5592,7 +6393,7 @@ private: if (it.isTextureSeam()) seamLength += l; } - face = m_nextPlanarRegionFace[face]; + face = m_planarCharts.nextRegionFace(face); if (face == firstFace) break; } @@ -5606,8 +6407,8 @@ private: float area = chart->area; uint32_t face = firstFace; for (;;) { - area += m_faceAreas[face]; - face = m_nextPlanarRegionFace[face]; + area += m_data.faceAreas[face]; + face = m_planarCharts.nextRegionFace(face); if (face == firstFace) break; } @@ -5618,21 +6419,21 @@ private: { float boundaryLength = chart->boundaryLength; // Add new edges, subtract edges shared with the chart. - const uint32_t planarRegionId = m_facePlanarRegionId[firstFace]; + const uint32_t planarRegionId = m_planarCharts.regionIdFromFace(firstFace); uint32_t face = firstFace; for (;;) { - for (Mesh::FaceEdgeIterator it(m_mesh, face); !it.isDone(); it.advance()) { - const float edgeLength = m_edgeLengths[it.edge()]; + for (Mesh::FaceEdgeIterator it(m_data.mesh, face); !it.isDone(); it.advance()) { + const float edgeLength = m_data.edgeLengths[it.edge()]; if (it.isBoundary()) { boundaryLength += edgeLength; - } else if (m_facePlanarRegionId[it.oppositeFace()] != planarRegionId) { + } else if (m_planarCharts.regionIdFromFace(it.oppositeFace()) != planarRegionId) { if (m_faceCharts[it.oppositeFace()] != chart->id) boundaryLength += edgeLength; else boundaryLength -= edgeLength; } } - face = m_nextPlanarRegionFace[face]; + face = m_planarCharts.nextRegionFace(face); if (face == firstFace) break; } @@ -5656,7 +6457,7 @@ private: m_faceCharts[chart->faces[i]] = chart->id; return false; } - if (dot(basis.normal, m_faceNormals[owner->faces[0]]) < 0.0f) // Flip normal if oriented in the wrong direction. + if (dot(basis.normal, m_data.faceNormals[owner->faces[0]]) < 0.0f) // Flip normal if oriented in the wrong direction. basis.normal = -basis.normal; // Compute orthogonal parameterization and check that it is valid. parameterizeChart(owner); @@ -5679,19 +6480,14 @@ private: return true; } - const Mesh *m_mesh; - Array m_edgeLengths; - Array m_faceAreas; - Array m_faceNormals; +private: + AtlasData &m_data; + const PlanarCharts &m_planarCharts; Array m_texcoords; uint32_t m_facesLeft; Array m_faceCharts; Array m_charts; CostQueue m_bestTriangles; - KISSRng m_rand; - ChartOptions m_options; - Array m_nextPlanarRegionFace; - Array m_facePlanarRegionId; Array m_tempPoints; UniformGrid2 m_boundaryGrid; #if XA_MERGE_CHARTS @@ -5700,231 +6496,63 @@ private: Array m_sharedBoundaryLengthsNoSeams; Array m_sharedBoundaryEdgeCountNoSeams; #endif + bool m_placingSeeds; }; -} // namespace segment - -namespace param { - -class JacobiPreconditioner +struct Atlas { -public: - JacobiPreconditioner(const sparse::Matrix &M, bool symmetric) : m_inverseDiagonal(M.width()) - { - XA_ASSERT(M.isSquare()); - for (uint32_t x = 0; x < M.width(); x++) { - float elem = M.getCoefficient(x, x); - //XA_DEBUG_ASSERT( elem != 0.0f ); // This can be zero in the presence of zero area triangles. - if (symmetric) { - m_inverseDiagonal[x] = (elem != 0) ? 1.0f / sqrtf(fabsf(elem)) : 1.0f; - } else { - m_inverseDiagonal[x] = (elem != 0) ? 1.0f / elem : 1.0f; - } - } - } + Atlas() : m_planarCharts(m_data), m_clusteredCharts(m_data, m_planarCharts) {} - void apply(const FullVector &x, FullVector &y) const + uint32_t chartCount() const { - XA_DEBUG_ASSERT(x.dimension() == m_inverseDiagonal.dimension()); - XA_DEBUG_ASSERT(y.dimension() == m_inverseDiagonal.dimension()); - // @@ Wrap vector component-wise product into a separate function. - const uint32_t D = x.dimension(); - for (uint32_t i = 0; i < D; i++) { - y[i] = m_inverseDiagonal[i] * x[i]; - } + return m_planarCharts.chartCount() + m_clusteredCharts.chartCount(); } -private: - FullVector m_inverseDiagonal; -}; + ConstArrayView chartFaces(uint32_t chartIndex) const + { + if (chartIndex < m_planarCharts.chartCount()) + return m_planarCharts.chartFaces(chartIndex); + chartIndex -= m_planarCharts.chartCount(); + return m_clusteredCharts.chartFaces(chartIndex); + } -// Linear solvers. -class Solver -{ -public: - // Solve the symmetric system: At·A·x = At·b - static bool LeastSquaresSolver(const sparse::Matrix &A, const FullVector &b, FullVector &x, float epsilon = 1e-5f) - { - XA_DEBUG_ASSERT(A.width() == x.dimension()); - XA_DEBUG_ASSERT(A.height() == b.dimension()); - XA_DEBUG_ASSERT(A.height() >= A.width()); // @@ If height == width we could solve it directly... - const uint32_t D = A.width(); - sparse::Matrix At(A.height(), A.width()); - sparse::transpose(A, At); - FullVector Atb(D); - sparse::mult(At, b, Atb); - sparse::Matrix AtA(D); - sparse::mult(At, A, AtA); - return SymmetricSolver(AtA, Atb, x, epsilon); - } - - // See section 10.4.3 in: Mesh Parameterization: Theory and Practice, Siggraph Course Notes, August 2007 - static bool LeastSquaresSolver(const sparse::Matrix &A, const FullVector &b, FullVector &x, const uint32_t *lockedParameters, uint32_t lockedCount, float epsilon = 1e-5f) - { - XA_DEBUG_ASSERT(A.width() == x.dimension()); - XA_DEBUG_ASSERT(A.height() == b.dimension()); - XA_DEBUG_ASSERT(A.height() >= A.width() - lockedCount); - // @@ This is not the most efficient way of building a system with reduced degrees of freedom. It would be faster to do it on the fly. - const uint32_t D = A.width() - lockedCount; - XA_DEBUG_ASSERT(D > 0); - // Compute: b - Al * xl - FullVector b_Alxl(b); - for (uint32_t y = 0; y < A.height(); y++) { - const uint32_t count = A.getRow(y).size(); - for (uint32_t e = 0; e < count; e++) { - uint32_t column = A.getRow(y)[e].x; - bool isFree = true; - for (uint32_t i = 0; i < lockedCount; i++) { - isFree &= (lockedParameters[i] != column); - } - if (!isFree) { - b_Alxl[y] -= x[column] * A.getRow(y)[e].v; - } - } - } - // Remove locked columns from A. - sparse::Matrix Af(D, A.height()); - for (uint32_t y = 0; y < A.height(); y++) { - const uint32_t count = A.getRow(y).size(); - for (uint32_t e = 0; e < count; e++) { - uint32_t column = A.getRow(y)[e].x; - uint32_t ix = column; - bool isFree = true; - for (uint32_t i = 0; i < lockedCount; i++) { - isFree &= (lockedParameters[i] != column); - if (column > lockedParameters[i]) ix--; // shift columns - } - if (isFree) { - Af.setCoefficient(ix, y, A.getRow(y)[e].v); - } - } - } - // Remove elements from x - FullVector xf(D); - for (uint32_t i = 0, j = 0; i < A.width(); i++) { - bool isFree = true; - for (uint32_t l = 0; l < lockedCount; l++) { - isFree &= (lockedParameters[l] != i); - } - if (isFree) { - xf[j++] = x[i]; - } - } - // Solve reduced system. - bool result = LeastSquaresSolver(Af, b_Alxl, xf, epsilon); - // Copy results back to x. - for (uint32_t i = 0, j = 0; i < A.width(); i++) { - bool isFree = true; - for (uint32_t l = 0; l < lockedCount; l++) { - isFree &= (lockedParameters[l] != i); - } - if (isFree) { - x[i] = xf[j++]; - } - } - return result; + const Basis &chartBasis(uint32_t chartIndex) const + { + if (chartIndex < m_planarCharts.chartCount()) + return m_planarCharts.chartBasis(chartIndex); + chartIndex -= m_planarCharts.chartCount(); + return m_clusteredCharts.chartBasis(chartIndex); } -private: - /** - * Compute the solution of the sparse linear system Ab=x using the Conjugate - * Gradient method. - * - * Solving sparse linear systems: - * (1) A·x = b - * - * The conjugate gradient algorithm solves (1) only in the case that A is - * symmetric and positive definite. It is based on the idea of minimizing the - * function - * - * (2) f(x) = 1/2·x·A·x - b·x - * - * This function is minimized when its gradient - * - * (3) df = A·x - b - * - * is zero, which is equivalent to (1). The minimization is carried out by - * generating a succession of search directions p.k and improved minimizers x.k. - * At each stage a quantity alfa.k is found that minimizes f(x.k + alfa.k·p.k), - * and x.k+1 is set equal to the new point x.k + alfa.k·p.k. The p.k and x.k are - * built up in such a way that x.k+1 is also the minimizer of f over the whole - * vector space of directions already taken, {p.1, p.2, . . . , p.k}. After N - * iterations you arrive at the minimizer over the entire vector space, i.e., the - * solution to (1). - * - * For a really good explanation of the method see: - * - * "An Introduction to the Conjugate Gradient Method Without the Agonizing Pain", - * Jonhathan Richard Shewchuk. - * - **/ - // Conjugate gradient with preconditioner. - static bool ConjugateGradientSolver(const JacobiPreconditioner &preconditioner, const sparse::Matrix &A, const FullVector &b, FullVector &x, float epsilon) - { - XA_DEBUG_ASSERT( A.isSquare() ); - XA_DEBUG_ASSERT( A.width() == b.dimension() ); - XA_DEBUG_ASSERT( A.width() == x.dimension() ); - int i = 0; - const int D = A.width(); - const int i_max = 4 * D; // Convergence should be linear, but in some cases, it's not. - FullVector r(D); // residual - FullVector p(D); // search direction - FullVector q(D); // - FullVector s(D); // preconditioned - float delta_0; - float delta_old; - float delta_new; - float alpha; - float beta; - // r = b - A·x - sparse::copy(b, r); - sparse::sgemv(-1, A, x, 1, r); - // p = M^-1 · r - preconditioner.apply(r, p); - delta_new = sparse::dot(r, p); - delta_0 = delta_new; - while (i < i_max && delta_new > epsilon * epsilon * delta_0) { - i++; - // q = A·p - sparse::mult(A, p, q); - // alpha = delta_new / p·q - const float pdotq = sparse::dot(p, q); - if (!isFinite(pdotq) || isNan(pdotq)) - alpha = 0.0f; - else - alpha = delta_new / pdotq; - // x = alfa·p + x - sparse::saxpy(alpha, p, x); - if ((i & 31) == 0) { // recompute r after 32 steps - // r = b - A·x - sparse::copy(b, r); - sparse::sgemv(-1, A, x, 1, r); - } else { - // r = r - alfa·q - sparse::saxpy(-alpha, q, r); - } - // s = M^-1 · r - preconditioner.apply(r, s); - delta_old = delta_new; - delta_new = sparse::dot( r, s ); - beta = delta_new / delta_old; - // p = s + beta·p - sparse::scal(beta, p); - sparse::saxpy(1, s, p); - } - return delta_new <= epsilon * epsilon * delta_0; + void reset(const Mesh *mesh, const ChartOptions &options) + { + XA_PROFILE_START(buildAtlasInit) + m_data.options = options; + m_data.mesh = mesh; + m_data.compute(); + XA_PROFILE_END(buildAtlasInit) } - static bool SymmetricSolver(const sparse::Matrix &A, const FullVector &b, FullVector &x, float epsilon = 1e-5f) + void compute() { - XA_DEBUG_ASSERT(A.height() == A.width()); - XA_DEBUG_ASSERT(A.height() == b.dimension()); - XA_DEBUG_ASSERT(b.dimension() == x.dimension()); - JacobiPreconditioner jacobi(A, true); - return ConjugateGradientSolver(jacobi, A, b, x, epsilon); + XA_PROFILE_START(planarCharts) + m_planarCharts.compute(); + XA_PROFILE_END(planarCharts) + XA_PROFILE_START(clusteredCharts) + m_clusteredCharts.compute(); + XA_PROFILE_END(clusteredCharts) } + +private: + AtlasData m_data; + PlanarCharts m_planarCharts; + ClusteredCharts m_clusteredCharts; }; +} // namespace segment + +namespace param { + // Fast sweep in 3 directions static bool findApproximateDiameterVertices(Mesh *mesh, uint32_t *a, uint32_t *b) { @@ -5982,134 +6610,95 @@ static bool findApproximateDiameterVertices(Mesh *mesh, uint32_t *a, uint32_t *b return true; } -// Conformal relations from Brecht Van Lommel (based on ABF): - -static float vec_angle_cos(const Vector3 &v1, const Vector3 &v2, const Vector3 &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(const Vector3 &v1, const Vector3 &v2, const Vector3 &v3) -{ - float dot = vec_angle_cos(v1, v2, v3); - return acosf(dot); -} - -static void triangle_angles(const Vector3 &v1, const Vector3 &v2, const Vector3 &v3, float *a1, float *a2, float *a3) -{ - *a1 = vec_angle(v3, v1, v2); - *a2 = vec_angle(v1, v2, v3); - *a3 = kPi - *a2 - *a1; -} - -static void setup_abf_relations(sparse::Matrix &A, int row, int id0, int id1, int id2, const Vector3 &p0, const Vector3 &p1, const Vector3 &p2) -{ - // @@ 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); - 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); - 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); +// From OpenNL LSCM example. +// Computes the coordinates of the vertices of a triangle in a local 2D orthonormal basis of the triangle's plane. +static void projectTriangle(Vector3 p0, Vector3 p1, Vector3 p2, Vector2 *z0, Vector2 *z1, Vector2 *z2, float epsilon) +{ + Vector3 X = normalize(p1 - p0, epsilon); + Vector3 Z = normalize(cross(X, p2 - p0), epsilon); + Vector3 Y = cross(Z, X); + Vector3 &O = p0; + *z0 = Vector2(0, 0); + *z1 = Vector2(length(p1 - O), 0); + *z2 = Vector2(dot(p2 - O, X), dot(p2 - O, Y)); } static bool computeLeastSquaresConformalMap(Mesh *mesh) { - // 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 uint32_t vertexCount = mesh->vertexCount(); - const uint32_t D = 2 * vertexCount; - const uint32_t N = 2 * mesh->faceCount(); - // 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; - } - sparse::Matrix A(D, N); - FullVector b(N); - FullVector x(D); - // Fill b: - b.fill(0.0f); - // Fill x: - uint32_t v0, v1; - if (!findApproximateDiameterVertices(mesh, &v0, &v1)) { + uint32_t lockedVertex0, lockedVertex1; + if (!findApproximateDiameterVertices(mesh, &lockedVertex0, &lockedVertex1)) { // Mesh has no boundaries. return false; } - if (mesh->texcoord(v0) == mesh->texcoord(v1)) { - // LSCM expects an existing parameterization. + const uint32_t vertexCount = mesh->vertexCount(); + opennl::NLContext *context = opennl::nlNewContext(); + opennl::nlSolverParameteri(context, NL_NB_VARIABLES, int(2 * vertexCount)); + opennl::nlSolverParameteri(context, NL_MAX_ITERATIONS, int(5 * vertexCount)); + opennl::nlBegin(context, NL_SYSTEM); + const Vector2 *texcoords = mesh->texcoords(); + for (uint32_t i = 0; i < vertexCount; i++) { + opennl::nlSetVariable(context, 2 * i, texcoords[i].x); + opennl::nlSetVariable(context, 2 * i + 1, texcoords[i].y); + if (i == lockedVertex0 || i == lockedVertex1) { + opennl::nlLockVariable(context, 2 * i); + opennl::nlLockVariable(context, 2 * i + 1); + } + } + opennl::nlBegin(context, NL_MATRIX); + const uint32_t faceCount = mesh->faceCount(); + const Vector3 *positions = mesh->positions(); + const uint32_t *indices = mesh->indices(); + for (uint32_t f = 0; f < faceCount; f++) { + const uint32_t v0 = indices[f * 3 + 0]; + const uint32_t v1 = indices[f * 3 + 1]; + const uint32_t v2 = indices[f * 3 + 2]; + Vector2 z0, z1, z2; + projectTriangle(positions[v0], positions[v1], positions[v2], &z0, &z1, &z2, mesh->epsilon()); + double a = z1.x - z0.x; + double b = z1.y - z0.y; + double c = z2.x - z0.x; + double d = z2.y - z0.y; + XA_DEBUG_ASSERT(b == 0.0); + // Note : 2*id + 0 --> u + // 2*id + 1 --> v + uint32_t u0_id = 2 * v0; + uint32_t v0_id = 2 * v0 + 1; + uint32_t u1_id = 2 * v1; + uint32_t v1_id = 2 * v1 + 1; + uint32_t u2_id = 2 * v2; + uint32_t v2_id = 2 * v2 + 1; + // Note : b = 0 + // Real part + opennl::nlBegin(context, NL_ROW); + opennl::nlCoefficient(context, u0_id, -a+c) ; + opennl::nlCoefficient(context, v0_id, b-d) ; + opennl::nlCoefficient(context, u1_id, -c) ; + opennl::nlCoefficient(context, v1_id, d) ; + opennl::nlCoefficient(context, u2_id, a); + opennl::nlEnd(context, NL_ROW); + // Imaginary part + opennl::nlBegin(context, NL_ROW); + opennl::nlCoefficient(context, u0_id, -b+d); + opennl::nlCoefficient(context, v0_id, -a+c); + opennl::nlCoefficient(context, u1_id, -d); + opennl::nlCoefficient(context, v1_id, -c); + opennl::nlCoefficient(context, v2_id, a); + opennl::nlEnd(context, NL_ROW); + } + opennl::nlEnd(context, NL_MATRIX); + opennl::nlEnd(context, NL_SYSTEM); + if (!opennl::nlSolve(context)) { + opennl::nlDeleteContext(context); return false; } - for (uint32_t v = 0; v < vertexCount; v++) { - // Initial solution. - x[2 * v + 0] = mesh->texcoord(v).x; - x[2 * v + 1] = mesh->texcoord(v).y; - } - // Fill A: - const uint32_t faceCount = mesh->faceCount(); - for (uint32_t f = 0, t = 0; f < faceCount; f++) { - const uint32_t vertex0 = mesh->vertexAt(f * 3 + 0); - const uint32_t vertex1 = mesh->vertexAt(f * 3 + 1); - const uint32_t vertex2 = mesh->vertexAt(f * 3 + 2); - setup_abf_relations(A, t, vertex0, vertex1, vertex2, mesh->position(vertex0), mesh->position(vertex1), mesh->position(vertex2)); - t++; - } - const uint32_t lockedParameters[] = { - 2 * v0 + 0, - 2 * v0 + 1, - 2 * v1 + 0, - 2 * v1 + 1 - }; - // Solve - Solver::LeastSquaresSolver(A, b, x, lockedParameters, 4, 0.000001f); - // Map x back to texcoords: - for (uint32_t v = 0; v < vertexCount; v++) { - mesh->texcoord(v) = Vector2(x[2 * v + 0], x[2 * v + 1]); - XA_DEBUG_ASSERT(!isNan(mesh->texcoord(v).x)); - XA_DEBUG_ASSERT(!isNan(mesh->texcoord(v).y)); + for (uint32_t i = 0; i < vertexCount; i++) { + const double u = opennl::nlGetVariable(context, 2 * i); + const double v = opennl::nlGetVariable(context, 2 * i + 1); + mesh->texcoord(i) = Vector2((float)u, (float)v); + XA_DEBUG_ASSERT(!isNan(mesh->texcoord(i).x)); + XA_DEBUG_ASSERT(!isNan(mesh->texcoord(i).y)); } + opennl::nlDeleteContext(context); return true; } @@ -6123,12 +6712,13 @@ struct PiecewiseParam const uint32_t vertexCount = m_mesh->vertexCount(); m_texcoords.resize(vertexCount); m_patch.reserve(m_faceCount); - m_faceAssigned.resize(m_faceCount); - m_faceAssigned.zeroOutMemory(); + m_candidates.reserve(m_faceCount); + m_faceInAnyPatch.resize(m_faceCount); + m_faceInAnyPatch.zeroOutMemory(); m_faceInvalid.resize(m_faceCount); m_faceInPatch.resize(m_faceCount); m_vertexInPatch.resize(vertexCount); - m_faceInCandidates.resize(m_faceCount); + m_faceToCandidate.resize(m_faceCount); } ConstArrayView chartFaces() const { return m_patch; } @@ -6136,19 +6726,20 @@ struct PiecewiseParam bool computeChart() { + // Clear per-patch state. m_patch.clear(); + m_candidates.clear(); + m_faceToCandidate.zeroOutMemory(); m_faceInvalid.zeroOutMemory(); m_faceInPatch.zeroOutMemory(); m_vertexInPatch.zeroOutMemory(); // Add the seed face (first unassigned face) to the patch. uint32_t seed = UINT32_MAX; for (uint32_t f = 0; f < m_faceCount; f++) { - if (m_faceAssigned.get(f)) + if (m_faceInAnyPatch.get(f)) continue; seed = f; - m_patch.push_back(seed); - m_faceInPatch.set(seed); - m_faceAssigned.set(seed); + // Add all 3 vertices. Vector2 texcoords[3]; orthoProjectFace(seed, texcoords); for (uint32_t i = 0; i < 3; i++) { @@ -6156,86 +6747,94 @@ struct PiecewiseParam m_vertexInPatch.set(vertex); m_texcoords[vertex] = texcoords[i]; } + addFaceToPatch(seed); + // Initialize the boundary grid. + m_boundaryGrid.reset(m_texcoords.data(), m_mesh->indices()); + for (Mesh::FaceEdgeIterator it(m_mesh, seed); !it.isDone(); it.advance()) + m_boundaryGrid.append(it.edge()); break; } if (seed == UINT32_MAX) return false; for (;;) { - findCandidates(); - if (m_candidates.isEmpty()) - break; - for (;;) { - // Find the candidate with the lowest cost. - float lowestCost = FLT_MAX; - uint32_t bestCandidate = UINT32_MAX; - for (uint32_t i = 0; i < m_candidates.size(); i++) { - const Candidate &candidate = m_candidates[i]; - if (m_faceInvalid.get(candidate.face)) // A candidate face may be invalidated after is was added. - continue; - if (candidate.maxCost < lowestCost) { - lowestCost = candidate.maxCost; - bestCandidate = i; - } + // Find the candidate with the lowest cost. + float lowestCost = FLT_MAX; + Candidate *bestCandidate = nullptr; + for (uint32_t i = 0; i < m_candidates.size(); i++) { + Candidate *candidate = m_candidates[i]; + if (candidate->maxCost < lowestCost) { + lowestCost = candidate->maxCost; + bestCandidate = candidate; } - if (bestCandidate == UINT32_MAX) + } + if (!bestCandidate) + break; + XA_DEBUG_ASSERT(!bestCandidate->prev); // Must be head of linked candidates. + // Compute the position by averaging linked candidates (candidates that share the same free vertex). + Vector2 position(0.0f); + uint32_t n = 0; + for (CandidateIterator it(bestCandidate); !it.isDone(); it.advance()) { + position += it.current()->position; + n++; + } + position *= 1.0f / (float)n; + const uint32_t freeVertex = bestCandidate->vertex; + XA_DEBUG_ASSERT(!isNan(position.x)); + XA_DEBUG_ASSERT(!isNan(position.y)); + m_texcoords[freeVertex] = position; + // Check for flipped faces. This is also done when candidates are first added, but the averaged position of the free vertex is different now, so check again. + bool invalid = false; + for (CandidateIterator it(bestCandidate); !it.isDone(); it.advance()) { + const uint32_t vertex0 = m_mesh->vertexAt(meshEdgeIndex0(it.current()->patchEdge)); + const uint32_t vertex1 = m_mesh->vertexAt(meshEdgeIndex1(it.current()->patchEdge)); + const float freeVertexOrient = orientToEdge(m_texcoords[vertex0], m_texcoords[vertex1], position); + if ((it.current()->patchVertexOrient < 0.0f && freeVertexOrient < 0.0f) || (it.current()->patchVertexOrient > 0.0f && freeVertexOrient > 0.0f)) { + invalid = true; break; - // Compute the position by averaging linked candidates (candidates that share the same free vertex). - Vector2 position(0.0f); - uint32_t n = 0; - for (CandidateIterator it(m_candidates, bestCandidate); !it.isDone(); it.advance()) { - position += it.current().position; - n++; } - position *= 1.0f / (float)n; - const uint32_t freeVertex = m_candidates[bestCandidate].vertex; - XA_DEBUG_ASSERT(!isNan(position.x)); - XA_DEBUG_ASSERT(!isNan(position.y)); - m_texcoords[freeVertex] = position; - // Check for flipped faces. This is also done when candidates are first added, but the averaged position of the free vertex is different now, so check again. - bool invalid = false; - for (CandidateIterator it(m_candidates, bestCandidate); !it.isDone(); it.advance()) { - const uint32_t vertex0 = m_mesh->vertexAt(meshEdgeIndex0(it.current().patchEdge)); - const uint32_t vertex1 = m_mesh->vertexAt(meshEdgeIndex1(it.current().patchEdge)); - const float freeVertexOrient = orientToEdge(m_texcoords[vertex0], m_texcoords[vertex1], position); - if ((it.current().patchVertexOrient < 0.0f && freeVertexOrient < 0.0f) || (it.current().patchVertexOrient > 0.0f && freeVertexOrient > 0.0f)) { - invalid = true; - break; - } - } - // Check for boundary intersection. - if (!invalid) { - m_boundaryGrid.reset(m_texcoords.data(), m_mesh->indices()); - // Add edges on the patch boundary to the grid. - // Temporarily adding candidate faces to the patch makes it simpler to detect which edges are on the boundary. - const uint32_t oldPatchSize = m_patch.size(); - for (CandidateIterator it(m_candidates, bestCandidate); !it.isDone(); it.advance()) - m_patch.push_back(it.current().face); - for (uint32_t i = 0; i < m_patch.size(); i++) { - for (Mesh::FaceEdgeIterator it(m_mesh, m_patch[i]); !it.isDone(); it.advance()) { - const uint32_t oface = it.oppositeFace(); - if (oface == UINT32_MAX || oface >= m_faceCount || !m_faceInPatch.get(oface)) - m_boundaryGrid.append(it.edge()); - } + } + // Check for boundary intersection. + if (!invalid) { + XA_PROFILE_START(parameterizeChartsPiecewiseBoundaryIntersection) + // Test candidate edges that would form part of the new patch boundary. + // Ignore boundary edges that would become internal if the candidate faces were added to the patch. + Array newBoundaryEdges, ignoreEdges; + for (CandidateIterator candidateIt(bestCandidate); !candidateIt.isDone(); candidateIt.advance()) { + for (Mesh::FaceEdgeIterator it(m_mesh, candidateIt.current()->face); !it.isDone(); it.advance()) { + const uint32_t oface = it.oppositeFace(); + if (oface == UINT32_MAX || oface >= m_faceCount || !m_faceInPatch.get(oface)) + newBoundaryEdges.push_back(it.edge()); + if (oface != UINT32_MAX && oface < m_faceCount && m_faceInPatch.get(oface)) + ignoreEdges.push_back(it.oppositeEdge()); } - invalid = m_boundaryGrid.intersectSelf(m_mesh->epsilon()); - m_patch.resize(oldPatchSize); - } - if (invalid) { - // Mark all faces of linked candidates as invalid. - for (CandidateIterator it(m_candidates, bestCandidate); !it.isDone(); it.advance()) - m_faceInvalid.set(it.current().face); - continue; - } - // Add faces to the patch. - for (CandidateIterator it(m_candidates, bestCandidate); !it.isDone(); it.advance()) { - m_patch.push_back(it.current().face); - m_faceInPatch.set(it.current().face); - m_faceAssigned.set(it.current().face); } + invalid = m_boundaryGrid.intersect(m_mesh->epsilon(), newBoundaryEdges, ignoreEdges); + XA_PROFILE_END(parameterizeChartsPiecewiseBoundaryIntersection) + } + if (invalid) { + // Mark all faces of linked candidates as invalid. + for (CandidateIterator it(bestCandidate); !it.isDone(); it.advance()) + m_faceInvalid.set(it.current()->face); + removeLinkedCandidates(bestCandidate); + } else { // Add vertex to the patch. m_vertexInPatch.set(freeVertex); + // Add faces to the patch. + for (CandidateIterator it(bestCandidate); !it.isDone(); it.advance()) + addFaceToPatch(it.current()->face); // Successfully added candidate face(s) to patch. - break; + removeLinkedCandidates(bestCandidate); + // Reset the grid with all edges on the patch boundary. + XA_PROFILE_START(parameterizeChartsPiecewiseBoundaryIntersection) + m_boundaryGrid.reset(m_texcoords.data(), m_mesh->indices()); + for (uint32_t i = 0; i < m_patch.size(); i++) { + for (Mesh::FaceEdgeIterator it(m_mesh, m_patch[i]); !it.isDone(); it.advance()) { + const uint32_t oface = it.oppositeFace(); + if (oface == UINT32_MAX || oface >= m_faceCount || !m_faceInPatch.get(oface)) + m_boundaryGrid.append(it.edge()); + } + } + XA_PROFILE_END(parameterizeChartsPiecewiseBoundaryIntersection) } } return true; @@ -6245,7 +6844,7 @@ private: struct Candidate { uint32_t face, vertex; - uint32_t next; // The next candidate with the same vertex. + Candidate *prev, *next; // The previous/next candidate with the same vertex. Vector2 position; float cost; float maxCost; // Of all linked candidates. @@ -6255,88 +6854,68 @@ private: struct CandidateIterator { - CandidateIterator(Array &candidates, uint32_t first) : m_candidates(candidates), m_current(first) {} - void advance() { if (m_current != UINT32_MAX) m_current = m_candidates[m_current].next; } - bool isDone() const { return m_current == UINT32_MAX; } - Candidate ¤t() { return m_candidates[m_current]; } + CandidateIterator(Candidate *head) : m_current(head) { XA_DEBUG_ASSERT(!head->prev); } + void advance() { if (m_current != nullptr) { m_current = m_current->next; } } + bool isDone() const { return !m_current; } + Candidate *current() { return m_current; } private: - Array &m_candidates; - uint32_t m_current; + Candidate *m_current; }; const Mesh *m_mesh; uint32_t m_faceCount; Array m_texcoords; - Array m_candidates; - BitArray m_faceInCandidates; - Array m_patch; - BitArray m_faceAssigned; // Face is assigned to a previous chart or the current patch. - BitArray m_faceInPatch, m_vertexInPatch; + BitArray m_faceInAnyPatch; // Face is in a previous chart patch or the current patch. + Array m_candidates; // Incident faces to the patch. + Array m_faceToCandidate; + Array m_patch; // The current chart patch. + BitArray m_faceInPatch, m_vertexInPatch; // Face/vertex is in the current patch. BitArray m_faceInvalid; // Face cannot be added to the patch - flipped, cost too high or causes boundary intersection. UniformGrid2 m_boundaryGrid; - // Find candidate faces on the patch front. - void findCandidates() - { - m_candidates.clear(); - m_faceInCandidates.zeroOutMemory(); - for (uint32_t i = 0; i < m_patch.size(); i++) { - for (Mesh::FaceEdgeIterator it(m_mesh, m_patch[i]); !it.isDone(); it.advance()) { - const uint32_t oface = it.oppositeFace(); - if (oface == UINT32_MAX || oface >= m_faceCount || m_faceAssigned.get(oface) || m_faceInCandidates.get(oface)) - continue; - // Found an active edge on the patch front. - // Find the free vertex (the vertex that isn't on the active edge). - // Compute the orientation of the other patch face vertex to the active edge. - uint32_t freeVertex = UINT32_MAX; - float orient = 0.0f; - for (uint32_t j = 0; j < 3; j++) { - const uint32_t vertex = m_mesh->vertexAt(oface * 3 + j); - if (vertex != it.vertex0() && vertex != it.vertex1()) { - freeVertex = vertex; - orient = orientToEdge(m_texcoords[it.vertex0()], m_texcoords[it.vertex1()], m_texcoords[m_mesh->vertexAt(m_patch[i] * 3 + j)]); - break; - } - } - XA_DEBUG_ASSERT(freeVertex != UINT32_MAX); - // If the free vertex is already in the patch, the face is enclosed by the patch. Add the face to the patch - don't need to assign texcoords. - if (m_vertexInPatch.get(freeVertex)) { - freeVertex = UINT32_MAX; - m_patch.push_back(oface); - m_faceAssigned.set(oface); - continue; - } - // Check this here rather than above so faces enclosed by the patch are always added. - if (m_faceInvalid.get(oface)) - continue; - addCandidateFace(it.edge(), orient, oface, it.oppositeEdge(), freeVertex); - } - } - // Link candidates that share the same vertex. - for (uint32_t i = 0; i < m_candidates.size(); i++) { - if (m_candidates[i].next != UINT32_MAX) + void addFaceToPatch(uint32_t face) + { + XA_DEBUG_ASSERT(!m_faceInPatch.get(face)); + XA_DEBUG_ASSERT(!m_faceInAnyPatch.get(face)); + m_patch.push_back(face); + m_faceInPatch.set(face); + m_faceInAnyPatch.set(face); + // Find new candidate faces on the patch incident to the newly added face. + for (Mesh::FaceEdgeIterator it(m_mesh, face); !it.isDone(); it.advance()) { + const uint32_t oface = it.oppositeFace(); + if (oface == UINT32_MAX || oface >= m_faceCount || m_faceInAnyPatch.get(oface) || m_faceToCandidate[oface]) continue; - uint32_t current = i; - for (uint32_t j = i + 1; j < m_candidates.size(); j++) { - if (m_candidates[j].vertex == m_candidates[current].vertex) { - m_candidates[current].next = j; - current = j; + // Found an active edge on the patch front. + // Find the free vertex (the vertex that isn't on the active edge). + // Compute the orientation of the other patch face vertex to the active edge. + uint32_t freeVertex = UINT32_MAX; + float orient = 0.0f; + for (uint32_t j = 0; j < 3; j++) { + const uint32_t vertex = m_mesh->vertexAt(oface * 3 + j); + if (vertex != it.vertex0() && vertex != it.vertex1()) { + freeVertex = vertex; + orient = orientToEdge(m_texcoords[it.vertex0()], m_texcoords[it.vertex1()], m_texcoords[m_mesh->vertexAt(face * 3 + j)]); + break; } } - } - // Set max cost for linked candidates. - for (uint32_t i = 0; i < m_candidates.size(); i++) { - float maxCost = 0.0f; - for (CandidateIterator it(m_candidates, i); !it.isDone(); it.advance()) - maxCost = max(maxCost, it.current().cost); - for (CandidateIterator it(m_candidates, i); !it.isDone(); it.advance()) - it.current().maxCost = maxCost; + XA_DEBUG_ASSERT(freeVertex != UINT32_MAX); + // If the free vertex is already in the patch, the face is enclosed by the patch. Add the face to the patch - don't need to assign texcoords. + /*if (m_vertexInPatch.get(freeVertex)) { + freeVertex = UINT32_MAX; + addFaceToPatch(oface, false); + continue; + }*/ + // Check this here rather than above so faces enclosed by the patch are always added. + if (m_faceInvalid.get(oface)) + continue; + addCandidateFace(it.edge(), orient, oface, it.oppositeEdge(), freeVertex); } } void addCandidateFace(uint32_t patchEdge, float patchVertexOrient, uint32_t face, uint32_t edge, uint32_t freeVertex) { + XA_DEBUG_ASSERT(!m_faceToCandidate[face]); Vector2 texcoords[3]; orthoProjectFace(face, texcoords); // Find corresponding vertices between the patch edge and candidate edge. @@ -6357,8 +6936,9 @@ private: const Vector2 localEdgeVec = texcoords[localVertex1] - texcoords[localVertex0]; const float len1 = length(patchEdgeVec); const float len2 = length(localEdgeVec); + if (len1 <= 0.0f || len2 <= 0.0f) + return; // Zero length edge. const float scale = len1 / len2; - XA_ASSERT(scale > 0.0f); for (uint32_t i = 0; i < 3; i++) texcoords[i] *= scale; // Translate to the first vertex on the patch edge. @@ -6380,6 +6960,8 @@ private: uv.x = x + texcoords[localVertex0].x; uv.y = y + texcoords[localVertex0].y; } + if (isNan(texcoords[localFreeVertex].x) || isNan(texcoords[localFreeVertex].y)) + return; // Check for local overlap (flipped triangle). // The patch face vertex that isn't on the active edge and the free vertex should be oriented on opposite sides to the active edge. const float freeVertexOrient = orientToEdge(m_texcoords[vertex0], m_texcoords[vertex1], texcoords[localFreeVertex]); @@ -6400,16 +6982,68 @@ private: } #endif // Add the candidate. - Candidate candidate; - candidate.face = face; - candidate.vertex = freeVertex; - candidate.position = texcoords[localFreeVertex]; - candidate.next = UINT32_MAX; - candidate.cost = cost; - candidate.patchEdge = patchEdge; - candidate.patchVertexOrient = patchVertexOrient; + Candidate *candidate = XA_ALLOC(MemTag::Default, Candidate); + candidate->face = face; + candidate->vertex = freeVertex; + candidate->position = texcoords[localFreeVertex]; + candidate->prev = candidate->next = nullptr; + candidate->cost = candidate->maxCost = cost; + candidate->patchEdge = patchEdge; + candidate->patchVertexOrient = patchVertexOrient; m_candidates.push_back(candidate); - m_faceInCandidates.set(face); + m_faceToCandidate[face] = candidate; + // Link with candidates that share the same vertex. Append to tail. + for (uint32_t i = 0; i < m_candidates.size() - 1; i++) { + if (m_candidates[i]->vertex == candidate->vertex) { + Candidate *tail = m_candidates[i]; + for (;;) { + if (tail->next) + tail = tail->next; + else + break; + } + candidate->prev = tail; + candidate->next = nullptr; + tail->next = candidate; + break; + } + } + // Set max cost for linked candidates. + Candidate *head = linkedCandidateHead(candidate); + float maxCost = 0.0f; + for (CandidateIterator it(head); !it.isDone(); it.advance()) + maxCost = max(maxCost, it.current()->cost); + for (CandidateIterator it(head); !it.isDone(); it.advance()) + it.current()->maxCost = maxCost; + } + + Candidate *linkedCandidateHead(Candidate *candidate) + { + Candidate *current = candidate; + for (;;) { + if (!current->prev) + break; + current = current->prev; + } + return current; + } + + void removeLinkedCandidates(Candidate *head) + { + XA_DEBUG_ASSERT(!head->prev); + Candidate *current = head; + while (current) { + Candidate *next = current->next; + m_faceToCandidate[current->face] = nullptr; + for (uint32_t i = 0; i < m_candidates.size(); i++) { + if (m_candidates[i] == current) { + m_candidates.removeAt(i); + break; + } + } + XA_FREE(current); + current = next; + } } void orthoProjectFace(uint32_t face, Vector2 *texcoords) const @@ -6479,7 +7113,7 @@ struct Quality boundaryGrid.reset(mesh->texcoords(), mesh->indices(), boundaryEdgeCount); for (uint32_t i = 0; i < boundaryEdgeCount; i++) boundaryGrid.append(boundaryEdges[i]); - boundaryIntersection = boundaryGrid.intersectSelf(mesh->epsilon()); + boundaryIntersection = boundaryGrid.intersect(mesh->epsilon()); #if XA_DEBUG_EXPORT_BOUNDARY_GRID static int exportIndex = 0; char filename[256]; @@ -6636,41 +7270,34 @@ struct ChartCtorBuffers Array boundaryLoops; }; -/// A chart is a connected set of faces with a certain topology (usually a disk). class Chart { public: - Chart(ChartCtorBuffers &buffers, const Basis &basis, ConstArrayView faces, const Mesh *originalMesh, uint32_t meshId, uint32_t chartGroupId, uint32_t chartId) : m_basis(basis), m_mesh(nullptr), m_unifiedMesh(nullptr), m_unmodifiedUnifiedMesh(nullptr), m_type(ChartType::LSCM), m_warningFlags(0), m_closedHolesCount(0), m_fixedTJunctionsCount(0) + Chart(ChartCtorBuffers &buffers, const ParameterizeOptions &options, const Basis &basis, ConstArrayView faces, const Mesh *sourceMesh, uint32_t chartGroupId, uint32_t chartId) : m_basis(basis), m_mesh(nullptr), m_unifiedMesh(nullptr), m_unmodifiedUnifiedMesh(nullptr), m_type(ChartType::LSCM), m_warningFlags(0), m_closedHolesCount(0), m_fixedTJunctionsCount(0), m_isInvalid(false) { - XA_UNUSED(meshId); XA_UNUSED(chartGroupId); XA_UNUSED(chartId); - m_faceArray.copyFrom(faces.data, faces.length); - // Copy face indices. - m_mesh = XA_NEW_ARGS(MemTag::Mesh, Mesh, originalMesh->epsilon(), m_faceArray.size() * 3, m_faceArray.size()); - m_unifiedMesh = XA_NEW_ARGS(MemTag::Mesh, Mesh, originalMesh->epsilon(), m_faceArray.size() * 3, m_faceArray.size()); - Array &chartMeshIndices = buffers.chartMeshIndices; - chartMeshIndices.resize(originalMesh->vertexCount()); - chartMeshIndices.setAll(UINT32_MAX); - Array &unifiedMeshIndices = buffers.unifiedMeshIndices; - unifiedMeshIndices.resize(originalMesh->vertexCount()); - unifiedMeshIndices.setAll(UINT32_MAX); + m_faceToSourceFaceMap.copyFrom(faces.data, faces.length); + const uint32_t approxVertexCount = min(faces.length * 3, sourceMesh->vertexCount()); + m_mesh = XA_NEW_ARGS(MemTag::Mesh, Mesh, sourceMesh->epsilon(), approxVertexCount, faces.length); + m_unifiedMesh = XA_NEW_ARGS(MemTag::Mesh, Mesh, sourceMesh->epsilon(), approxVertexCount, faces.length); + HashMap> sourceVertexToUnifiedVertexMap(MemTag::Mesh, approxVertexCount), sourceVertexToChartVertexMap(MemTag::Mesh, approxVertexCount); // Add vertices. - const uint32_t faceCount = m_initialFaceCount = m_faceArray.size(); + const uint32_t faceCount = m_initialFaceCount = faces.length; for (uint32_t f = 0; f < faceCount; f++) { for (uint32_t i = 0; i < 3; i++) { - const uint32_t vertex = originalMesh->vertexAt(m_faceArray[f] * 3 + i); - const uint32_t unifiedVertex = originalMesh->firstColocal(vertex); - if (unifiedMeshIndices[unifiedVertex] == (uint32_t)~0) { - unifiedMeshIndices[unifiedVertex] = m_unifiedMesh->vertexCount(); - XA_DEBUG_ASSERT(equal(originalMesh->position(vertex), originalMesh->position(unifiedVertex), originalMesh->epsilon())); - m_unifiedMesh->addVertex(originalMesh->position(vertex)); + const uint32_t sourceVertex = sourceMesh->vertexAt(m_faceToSourceFaceMap[f] * 3 + i); + const uint32_t sourceUnifiedVertex = sourceMesh->firstColocal(sourceVertex); + uint32_t unifiedVertex = sourceVertexToUnifiedVertexMap.get(sourceUnifiedVertex); + if (unifiedVertex == UINT32_MAX) { + unifiedVertex = sourceVertexToUnifiedVertexMap.add(sourceUnifiedVertex); + m_unifiedMesh->addVertex(sourceMesh->position(sourceVertex)); } - if (chartMeshIndices[vertex] == (uint32_t)~0) { - chartMeshIndices[vertex] = m_mesh->vertexCount(); - m_chartToOriginalMap.push_back(vertex); - m_chartToUnifiedMap.push_back(unifiedMeshIndices[unifiedVertex]); - m_mesh->addVertex(originalMesh->position(vertex), Vector3(0.0f), originalMesh->texcoord(vertex)); + if (sourceVertexToChartVertexMap.get(sourceVertex) == UINT32_MAX) { + sourceVertexToChartVertexMap.add(sourceVertex); + m_vertexToSourceVertexMap.push_back(sourceVertex); + m_chartVertexToUnifiedVertexMap.push_back(unifiedVertex); + m_mesh->addVertex(sourceMesh->position(sourceVertex), Vector3(0.0f), sourceMesh->texcoord(sourceVertex)); } } } @@ -6678,9 +7305,12 @@ public: for (uint32_t f = 0; f < faceCount; f++) { uint32_t indices[3], unifiedIndices[3]; for (uint32_t i = 0; i < 3; i++) { - const uint32_t vertex = originalMesh->vertexAt(m_faceArray[f] * 3 + i); - indices[i] = chartMeshIndices[vertex]; - unifiedIndices[i] = unifiedMeshIndices[originalMesh->firstColocal(vertex)]; + const uint32_t sourceVertex = sourceMesh->vertexAt(m_faceToSourceFaceMap[f] * 3 + i); + const uint32_t sourceUnifiedVertex = sourceMesh->firstColocal(sourceVertex); + indices[i] = sourceVertexToChartVertexMap.get(sourceVertex); + XA_DEBUG_ASSERT(indices[i] != UINT32_MAX); + unifiedIndices[i] = sourceVertexToUnifiedVertexMap.get(sourceUnifiedVertex); + XA_DEBUG_ASSERT(unifiedIndices[i] != UINT32_MAX); } Mesh::AddFaceResult::Enum result = m_mesh->addFace(indices); XA_UNUSED(result); @@ -6698,15 +7328,18 @@ public: XA_DEBUG_ASSERT(result == Mesh::AddFaceResult::OK); } m_mesh->createBoundaries(); // For AtlasPacker::computeBoundingBox + m_mesh->destroyEdgeMap(); // Only needed it for createBoundaries. m_unifiedMesh->createBoundaries(); - if (meshIsPlanar(*m_unifiedMesh)) + if (meshIsPlanar(*m_unifiedMesh)) { m_type = ChartType::Planar; - else { - m_unifiedMesh->linkBoundaries(); + return; + } + m_unifiedMesh->linkBoundaries(); #if XA_DEBUG_EXPORT_OBJ_BEFORE_FIX_TJUNCTION - m_unifiedMesh->writeObjFile("debug_before_fix_tjunction.obj"); + m_unifiedMesh->writeObjFile("debug_before_fix_tjunction.obj"); #endif - bool duplicatedEdge = false, failed = false; + bool duplicatedEdge = false, failed = false; + if (options.fixTJunctions) { XA_PROFILE_START(fixChartMeshTJunctions) Mesh *fixedUnifiedMesh = meshFixTJunctions(*m_unifiedMesh, &duplicatedEdge, &failed, &m_fixedTJunctionsCount); XA_PROFILE_END(fixChartMeshTJunctions) @@ -6721,6 +7354,8 @@ public: m_unifiedMesh->linkBoundaries(); m_initialFaceCount = m_unifiedMesh->faceCount(); // Fixing t-junctions rewrites faces. } + } + if (options.closeHoles) { // See if there are any holes that need closing. Array &boundaryLoops = buffers.boundaryLoops; meshGetBoundaryLoops(*m_unifiedMesh, boundaryLoops); @@ -6751,7 +7386,7 @@ public: #if XA_DEBUG_EXPORT_OBJ_CLOSE_HOLES_ERROR if (m_warningFlags & ChartWarningFlags::CloseHolesFailed) { char filename[256]; - XA_SPRINTF(filename, sizeof(filename), "debug_mesh_%03u_chartgroup_%03u_chart_%03u_close_holes_error.obj", meshId, chartGroupId, chartId); + XA_SPRINTF(filename, sizeof(filename), "debug_mesh_%03u_chartgroup_%03u_chart_%03u_close_holes_error.obj", sourceMesh->id(), chartGroupId, chartId); FILE *file; XA_FOPEN(file, filename, "w"); if (file) { @@ -6780,69 +7415,44 @@ public: } #if XA_RECOMPUTE_CHARTS - Chart(ChartCtorBuffers &buffers, const Chart *parent, const Mesh *parentMesh, ConstArrayView faces, const Vector2 *texcoords, const Mesh *originalMesh, uint32_t meshId, uint32_t chartGroupId, uint32_t chartId) : m_mesh(nullptr), m_unifiedMesh(nullptr), m_unmodifiedUnifiedMesh(nullptr), m_type(ChartType::Piecewise), m_warningFlags(0), m_closedHolesCount(0), m_fixedTJunctionsCount(0) + Chart(ChartCtorBuffers &buffers, const Chart *parent, const Mesh *parentMesh, ConstArrayView faces, const Vector2 *texcoords, const Mesh *sourceMesh) : m_mesh(nullptr), m_unifiedMesh(nullptr), m_unmodifiedUnifiedMesh(nullptr), m_type(ChartType::Piecewise), m_warningFlags(0), m_closedHolesCount(0), m_fixedTJunctionsCount(0), m_isInvalid(false) { - XA_UNUSED(meshId); - XA_UNUSED(chartGroupId); - XA_UNUSED(chartId); const uint32_t faceCount = m_initialFaceCount = faces.length; - m_faceArray.resize(faceCount); + m_faceToSourceFaceMap.resize(faceCount); for (uint32_t i = 0; i < faceCount; i++) - m_faceArray[i] = parent->m_faceArray[faces[i]]; // Map faces to parent chart original mesh. + m_faceToSourceFaceMap[i] = parent->m_faceToSourceFaceMap[faces[i]]; // Map faces to parent chart source mesh. // Copy face indices. - m_mesh = XA_NEW_ARGS(MemTag::Mesh, Mesh, originalMesh->epsilon(), m_faceArray.size() * 3, m_faceArray.size()); - m_unifiedMesh = XA_NEW_ARGS(MemTag::Mesh, Mesh, originalMesh->epsilon(), m_faceArray.size() * 3, m_faceArray.size()); + m_mesh = XA_NEW_ARGS(MemTag::Mesh, Mesh, sourceMesh->epsilon(), m_faceToSourceFaceMap.size() * 3, m_faceToSourceFaceMap.size()); Array &chartMeshIndices = buffers.chartMeshIndices; - chartMeshIndices.resize(originalMesh->vertexCount()); - chartMeshIndices.setAll(UINT32_MAX); - Array &unifiedMeshIndices = buffers.unifiedMeshIndices; - unifiedMeshIndices.resize(originalMesh->vertexCount()); - unifiedMeshIndices.setAll(UINT32_MAX); + chartMeshIndices.resize(sourceMesh->vertexCount()); + chartMeshIndices.fillBytes(0xff); // Add vertices. for (uint32_t f = 0; f < faceCount; f++) { for (uint32_t i = 0; i < 3; i++) { - const uint32_t vertex = originalMesh->vertexAt(m_faceArray[f] * 3 + i); - const uint32_t unifiedVertex = originalMesh->firstColocal(vertex); + const uint32_t vertex = sourceMesh->vertexAt(m_faceToSourceFaceMap[f] * 3 + i); const uint32_t parentVertex = parentMesh->vertexAt(faces[f] * 3 + i); - if (unifiedMeshIndices[unifiedVertex] == (uint32_t)~0) { - unifiedMeshIndices[unifiedVertex] = m_unifiedMesh->vertexCount(); - XA_DEBUG_ASSERT(equal(originalMesh->position(vertex), originalMesh->position(unifiedVertex), originalMesh->epsilon())); - m_unifiedMesh->addVertex(originalMesh->position(vertex), Vector3(0.0f), texcoords[parentVertex]); - } if (chartMeshIndices[vertex] == (uint32_t)~0) { chartMeshIndices[vertex] = m_mesh->vertexCount(); - m_chartToOriginalMap.push_back(vertex); - m_chartToUnifiedMap.push_back(unifiedMeshIndices[unifiedVertex]); - m_mesh->addVertex(originalMesh->position(vertex), Vector3(0.0f), texcoords[parentVertex]); + m_vertexToSourceVertexMap.push_back(vertex); + m_mesh->addVertex(sourceMesh->position(vertex), Vector3(0.0f), texcoords[parentVertex]); } } } // Add faces. for (uint32_t f = 0; f < faceCount; f++) { - uint32_t indices[3], unifiedIndices[3]; + uint32_t indices[3]; for (uint32_t i = 0; i < 3; i++) { - const uint32_t vertex = originalMesh->vertexAt(m_faceArray[f] * 3 + i); + const uint32_t vertex = sourceMesh->vertexAt(m_faceToSourceFaceMap[f] * 3 + i); indices[i] = chartMeshIndices[vertex]; - unifiedIndices[i] = unifiedMeshIndices[originalMesh->firstColocal(vertex)]; } Mesh::AddFaceResult::Enum result = m_mesh->addFace(indices); XA_UNUSED(result); XA_DEBUG_ASSERT(result == Mesh::AddFaceResult::OK); -#if XA_DEBUG - // Unifying colocals may create degenerate edges. e.g. if two triangle vertices are colocal. - for (int i = 0; i < 3; i++) { - const uint32_t index1 = unifiedIndices[i]; - const uint32_t index2 = unifiedIndices[(i + 1) % 3]; - XA_DEBUG_ASSERT(index1 != index2); - } -#endif - result = m_unifiedMesh->addFace(unifiedIndices); - XA_UNUSED(result); - XA_DEBUG_ASSERT(result == Mesh::AddFaceResult::OK); } m_mesh->createBoundaries(); // For AtlasPacker::computeBoundingBox - m_unifiedMesh->createBoundaries(); - m_unifiedMesh->linkBoundaries(); + m_mesh->destroyEdgeMap(); // Only needed it for createBoundaries. + // Need to store texcoords for backup/restore so packing can be run multiple times. + backupTexcoords(); } #endif @@ -6852,17 +7462,10 @@ public: m_mesh->~Mesh(); XA_FREE(m_mesh); } - if (m_unifiedMesh) { - m_unifiedMesh->~Mesh(); - XA_FREE(m_unifiedMesh); - } - if (m_unmodifiedUnifiedMesh) { - m_unmodifiedUnifiedMesh->~Mesh(); - XA_FREE(m_unmodifiedUnifiedMesh); - } + destroyUnifiedMesh(); } - const Basis &basis() const { return m_basis; } + bool isInvalid() const { return m_isInvalid; } ChartType::Enum type() const { return m_type; } uint32_t warningFlags() const { return m_warningFlags; } uint32_t closedHolesCount() const { return m_closedHolesCount; } @@ -6872,45 +7475,67 @@ public: #if XA_DEBUG_EXPORT_OBJ_INVALID_PARAMETERIZATION const Array ¶mFlippedFaces() const { return m_paramFlippedFaces; } #endif - uint32_t mapFaceToSourceFace(uint32_t i) const { return m_faceArray[i]; } + uint32_t mapFaceToSourceFace(uint32_t i) const { return m_faceToSourceFaceMap[i]; } + uint32_t mapChartVertexToSourceVertex(uint32_t i) const { return m_vertexToSourceVertexMap[i]; } const Mesh *mesh() const { return m_mesh; } Mesh *mesh() { return m_mesh; } const Mesh *unifiedMesh() const { return m_unifiedMesh; } - Mesh *unifiedMesh() { return m_unifiedMesh; } const Mesh *unmodifiedUnifiedMesh() const { return m_unmodifiedUnifiedMesh; } - uint32_t mapChartVertexToOriginalVertex(uint32_t i) const { return m_chartToOriginalMap[i]; } - - void evaluateOrthoQuality(UniformGrid2 &boundaryGrid) - { - XA_PROFILE_START(parameterizeChartsEvaluateQuality) - m_quality.computeBoundaryIntersection(m_unifiedMesh, boundaryGrid); - m_quality.computeFlippedFaces(m_unifiedMesh, m_initialFaceCount, nullptr); - m_quality.computeMetrics(m_unifiedMesh, m_initialFaceCount); - XA_PROFILE_END(parameterizeChartsEvaluateQuality) - // Use orthogonal parameterization if quality is acceptable. - if (!m_quality.boundaryIntersection && m_quality.totalGeometricArea > 0.0f && m_quality.stretchMetric <= 1.1f && m_quality.maxStretchMetric <= 1.25f) - m_type = ChartType::Ortho; - } - void evaluateQuality(UniformGrid2 &boundaryGrid) + void parameterize(const ParameterizeOptions &options, UniformGrid2 &boundaryGrid) { - XA_PROFILE_START(parameterizeChartsEvaluateQuality) - m_quality.computeBoundaryIntersection(m_unifiedMesh, boundaryGrid); + XA_PROFILE_START(parameterizeChartsOrthogonal) + { + // Project vertices to plane. + const uint32_t vertexCount = m_unifiedMesh->vertexCount(); + for (uint32_t i = 0; i < vertexCount; i++) + m_unifiedMesh->texcoord(i) = Vector2(dot(m_basis.tangent, m_unifiedMesh->position(i)), dot(m_basis.bitangent, m_unifiedMesh->position(i))); + } + XA_PROFILE_END(parameterizeChartsOrthogonal) + // Computing charts checks for flipped triangles and boundary intersection. Don't need to do that again here if chart is planar. + if (m_type != ChartType::Planar) { + XA_PROFILE_START(parameterizeChartsEvaluateQuality) + m_quality.computeBoundaryIntersection(m_unifiedMesh, boundaryGrid); + m_quality.computeFlippedFaces(m_unifiedMesh, m_initialFaceCount, nullptr); + m_quality.computeMetrics(m_unifiedMesh, m_initialFaceCount); + XA_PROFILE_END(parameterizeChartsEvaluateQuality) + // Use orthogonal parameterization if quality is acceptable. + if (!m_quality.boundaryIntersection && m_quality.flippedTriangleCount == 0 && m_quality.totalGeometricArea > 0.0f && m_quality.stretchMetric <= 1.1f && m_quality.maxStretchMetric <= 1.25f) + m_type = ChartType::Ortho; + } + if (m_type == ChartType::LSCM) { + XA_PROFILE_START(parameterizeChartsLSCM) + if (options.func) { + options.func(&m_unifiedMesh->position(0).x, &m_unifiedMesh->texcoord(0).x, m_unifiedMesh->vertexCount(), m_unifiedMesh->indices(), m_unifiedMesh->indexCount()); + } + else + computeLeastSquaresConformalMap(m_unifiedMesh); + XA_PROFILE_END(parameterizeChartsLSCM) + XA_PROFILE_START(parameterizeChartsEvaluateQuality) + m_quality.computeBoundaryIntersection(m_unifiedMesh, boundaryGrid); #if XA_DEBUG_EXPORT_OBJ_INVALID_PARAMETERIZATION - m_quality.computeFlippedFaces(m_unifiedMesh, m_initialFaceCount, &m_paramFlippedFaces); + m_quality.computeFlippedFaces(m_unifiedMesh, m_initialFaceCount, &m_paramFlippedFaces); #else - m_quality.computeFlippedFaces(m_unifiedMesh, m_initialFaceCount, nullptr); + m_quality.computeFlippedFaces(m_unifiedMesh, m_initialFaceCount, nullptr); #endif - // Don't need to call computeMetrics here, that's only used in evaluateOrthoQuality to determine if quality is acceptable enough to use ortho projection. - XA_PROFILE_END(parameterizeChartsEvaluateQuality) - } - - // Transfer parameterization from unified mesh to chart mesh. - void transferParameterization() - { + // Don't need to call computeMetrics here, that's only used in evaluateOrthoQuality to determine if quality is acceptable enough to use ortho projection. + if (m_quality.boundaryIntersection || m_quality.flippedTriangleCount > 0) + m_isInvalid = true; + XA_PROFILE_END(parameterizeChartsEvaluateQuality) + } +#if XA_DEBUG_ALL_CHARTS_INVALID + m_isInvalid = true; +#endif + // Transfer parameterization from unified mesh to chart mesh. const uint32_t vertexCount = m_mesh->vertexCount(); for (uint32_t v = 0; v < vertexCount; v++) - m_mesh->texcoord(v) = m_unifiedMesh->texcoord(m_chartToUnifiedMap[v]); + m_mesh->texcoord(v) = m_unifiedMesh->texcoord(m_chartVertexToUnifiedVertexMap[v]); + // Can destroy unified mesh now. + // But not if the parameterization is invalid, the unified mesh will be needed for PiecewiseParameterization. + if (!m_isInvalid) + destroyUnifiedMesh(); + // Need to store texcoords for backup/restore so packing can be run multiple times. + backupTexcoords(); } Vector2 computeParametricBounds() const @@ -6925,7 +7550,34 @@ public: return (maxCorner - minCorner) * 0.5f; } + void restoreTexcoords() + { + memcpy(m_mesh->texcoords(), m_backupTexcoords.data(), m_mesh->vertexCount() * sizeof(Vector2)); + } + private: + void backupTexcoords() + { + m_backupTexcoords.resize(m_mesh->vertexCount()); + memcpy(m_backupTexcoords.data(), m_mesh->texcoords(), m_mesh->vertexCount() * sizeof(Vector2)); + } + + void destroyUnifiedMesh() + { + if (m_unifiedMesh) { + m_unifiedMesh->~Mesh(); + XA_FREE(m_unifiedMesh); + m_unifiedMesh = nullptr; + } + if (m_unmodifiedUnifiedMesh) { + m_unmodifiedUnifiedMesh->~Mesh(); + XA_FREE(m_unmodifiedUnifiedMesh); + m_unmodifiedUnifiedMesh = nullptr; + } + // Don't need this when unified meshes are destroyed. + m_chartVertexToUnifiedVertexMap.destroy(); + } + Basis m_basis; Mesh *m_mesh; Mesh *m_unifiedMesh; @@ -6935,482 +7587,511 @@ private: uint32_t m_initialFaceCount; // Before fixing T-junctions and/or closing holes. uint32_t m_closedHolesCount, m_fixedTJunctionsCount; - // List of faces of the original mesh that belong to this chart. - Array m_faceArray; + // List of faces of the source mesh that belong to this chart. + Array m_faceToSourceFaceMap; - // Map vertices of the chart mesh to vertices of the original mesh. - Array m_chartToOriginalMap; + // Map vertices of the chart mesh to vertices of the source mesh. + Array m_vertexToSourceVertexMap; - Array m_chartToUnifiedMap; + Array m_chartVertexToUnifiedVertexMap; + + Array m_backupTexcoords; Quality m_quality; #if XA_DEBUG_EXPORT_OBJ_INVALID_PARAMETERIZATION Array m_paramFlippedFaces; #endif + bool m_isInvalid; }; -struct CreateChartTaskArgs +struct CreateAndParameterizeChartTaskArgs { - const Mesh *mesh; const Basis *basis; + ThreadLocal *boundaryGrid; + Chart *chart; // output + Array charts; // output (if more than one chart) + ThreadLocal *chartBuffers; + const Mesh *mesh; + const ParameterizeOptions *options; +#if XA_RECOMPUTE_CHARTS + ThreadLocal *pp; +#endif ConstArrayView faces; - uint32_t meshId; uint32_t chartGroupId; uint32_t chartId; - ThreadLocal *chartBuffers; - Chart **chart; }; -static void runCreateChartTask(void *userData) -{ - XA_PROFILE_START(createChartMeshesThread) - auto args = (CreateChartTaskArgs *)userData; - *(args->chart) = XA_NEW_ARGS(MemTag::Default, Chart, args->chartBuffers->get(), *(args->basis), args->faces, args->mesh, args->meshId, args->chartGroupId, args->chartId); - XA_PROFILE_END(createChartMeshesThread) -} - -struct ParameterizeChartTaskArgs +static void runCreateAndParameterizeChartTask(void *userData) { - Chart *chart; - ParameterizeFunc func; - ThreadLocal *boundaryGrid; -}; - -static void runParameterizeChartTask(void *userData) -{ - auto args = (ParameterizeChartTaskArgs *)userData; - Mesh *mesh = args->chart->unifiedMesh(); - XA_PROFILE_START(parameterizeChartsOrthogonal) - { - // Project vertices to plane. - const uint32_t vertexCount = mesh->vertexCount(); - const Basis &basis = args->chart->basis(); - for (uint32_t i = 0; i < vertexCount; i++) - mesh->texcoord(i) = Vector2(dot(basis.tangent, mesh->position(i)), dot(basis.bitangent, mesh->position(i))); - } - XA_PROFILE_END(parameterizeChartsOrthogonal) - // Computing charts checks for flipped triangles and boundary intersection. Don't need to do that again here if chart is planar. - if (args->chart->type() != ChartType::Planar) - args->chart->evaluateOrthoQuality(args->boundaryGrid->get()); - if (args->chart->type() == ChartType::LSCM) { - XA_PROFILE_START(parameterizeChartsLSCM) - if (args->func) - args->func(&mesh->position(0).x, &mesh->texcoord(0).x, mesh->vertexCount(), mesh->indices(), mesh->indexCount()); - else - computeLeastSquaresConformalMap(mesh); - XA_PROFILE_END(parameterizeChartsLSCM) - args->chart->evaluateQuality(args->boundaryGrid->get()); + auto args = (CreateAndParameterizeChartTaskArgs *)userData; + XA_PROFILE_START(createChartMesh) + args->chart = XA_NEW_ARGS(MemTag::Default, Chart, args->chartBuffers->get(), *args->options, *args->basis, args->faces, args->mesh, args->chartGroupId, args->chartId); + XA_PROFILE_END(createChartMesh) + args->chart->parameterize(*args->options, args->boundaryGrid->get()); +#if XA_RECOMPUTE_CHARTS + if (!args->chart->isInvalid()) + return; + // Recompute charts with invalid parameterizations. + XA_PROFILE_START(parameterizeChartsRecompute) + Chart *invalidChart = args->chart; + // Fixing t-junctions rewrites unified mesh faces, and we need to map faces back to input mesh. So use the unmodified unified mesh. + const Mesh *invalidMesh = invalidChart->unmodifiedUnifiedMesh(); + uint32_t faceCount = 0; + if (invalidMesh) { + faceCount = invalidMesh->faceCount(); + } else { + invalidMesh = invalidChart->unifiedMesh(); + faceCount = invalidChart->initialFaceCount(); // Not invalidMesh->faceCount(). Don't want faces added by hole closing. + } + PiecewiseParam &pp = args->pp->get(); + pp.reset(invalidMesh, faceCount); +#if XA_DEBUG_EXPORT_OBJ_RECOMPUTED_CHARTS + char filename[256]; + XA_SPRINTF(filename, sizeof(filename), "debug_mesh_%03u_chartgroup_%03u_chart_%03u_recomputed.obj", args->mesh->id(), args->chartGroupId, args->chartId); + FILE *file; + XA_FOPEN(file, filename, "w"); + uint32_t subChartIndex = 0; +#endif + for (;;) { + XA_PROFILE_START(parameterizeChartsPiecewise) + const bool facesRemaining = pp.computeChart(); + XA_PROFILE_END(parameterizeChartsPiecewise) + if (!facesRemaining) + break; + Chart *chart = XA_NEW_ARGS(MemTag::Default, Chart, args->chartBuffers->get(), invalidChart, invalidMesh, pp.chartFaces(), pp.texcoords(), args->mesh); + args->charts.push_back(chart); +#if XA_DEBUG_EXPORT_OBJ_RECOMPUTED_CHARTS + if (file) { + for (uint32_t j = 0; j < invalidMesh->vertexCount(); j++) { + fprintf(file, "v %g %g %g\n", invalidMesh->position(j).x, invalidMesh->position(j).y, invalidMesh->position(j).z); + fprintf(file, "vt %g %g\n", pp.texcoords()[j].x, pp.texcoords()[j].y); + } + fprintf(file, "o chart%03u\n", subChartIndex); + fprintf(file, "s off\n"); + for (uint32_t f = 0; f < pp.chartFaces().length; f++) { + fprintf(file, "f "); + const uint32_t face = pp.chartFaces()[f]; + for (uint32_t j = 0; j < 3; j++) { + const uint32_t index = invalidMesh->vertexCount() * subChartIndex + invalidMesh->vertexAt(face * 3 + j) + 1; // 1-indexed + fprintf(file, "%d/%d/%c", index, index, j == 2 ? '\n' : ' '); + } + } + } + subChartIndex++; +#endif } - // Transfer parameterization from unified mesh to chart mesh. - args->chart->transferParameterization(); +#if XA_DEBUG_EXPORT_OBJ_RECOMPUTED_CHARTS + if (file) + fclose(file); +#endif + XA_PROFILE_END(parameterizeChartsRecompute) +#endif // XA_RECOMPUTE_CHARTS } // Set of charts corresponding to mesh faces in the same face group. class ChartGroup { public: - ChartGroup(uint32_t id, const Mesh *sourceMesh, uint16_t faceGroup) : m_sourceId(sourceMesh->id()), m_id(id), m_isVertexMap(faceGroup == Mesh::kInvalidFaceGroup), m_paramAddedChartsCount(0), m_paramDeletedChartsCount(0) + ChartGroup(uint32_t id, const Mesh *sourceMesh, const MeshFaceGroups *sourceMeshFaceGroups, MeshFaceGroups::Handle faceGroup) : m_id(id), m_sourceMesh(sourceMesh), m_sourceMeshFaceGroups(sourceMeshFaceGroups), m_faceGroup(faceGroup), m_faceCount(0), m_paramAddedChartsCount(0), m_paramDeletedChartsCount(0) { - // Create new mesh from the source mesh, using faces that belong to this group. - const uint32_t sourceFaceCount = sourceMesh->faceCount(); - if (!m_isVertexMap) { - m_faceToSourceFaceMap.reserve(sourceMesh->faceGroupFaceCount(faceGroup)); - for (Mesh::GroupFaceIterator it(sourceMesh, faceGroup); !it.isDone(); it.advance()) - m_faceToSourceFaceMap.push_back(it.face()); - } else { - for (uint32_t f = 0; f < sourceFaceCount; f++) { - if (sourceMesh->faceGroupAt(f) == faceGroup) - m_faceToSourceFaceMap.push_back(f); - } - } - // Only initial meshes have face groups and ignored faces. The only flag we care about is HasNormals. - const uint32_t faceCount = m_faceToSourceFaceMap.size(); - XA_DEBUG_ASSERT(faceCount > 0); - const uint32_t approxVertexCount = faceCount * 3; - m_mesh = XA_NEW_ARGS(MemTag::Mesh, Mesh, sourceMesh->epsilon(), approxVertexCount, faceCount, sourceMesh->flags() & MeshFlags::HasNormals); - m_vertexToSourceVertexMap.reserve(approxVertexCount); - HashMap sourceVertexToVertexMap(MemTag::Mesh, approxVertexCount); - for (uint32_t f = 0; f < faceCount; f++) { - const uint32_t face = m_faceToSourceFaceMap[f]; - for (uint32_t i = 0; i < 3; i++) { - const uint32_t vertex = sourceMesh->vertexAt(face * 3 + i); - if (sourceVertexToVertexMap.get(vertex) == UINT32_MAX) { - sourceVertexToVertexMap.add(vertex); - m_vertexToSourceVertexMap.push_back(vertex); - Vector3 normal(0.0f); - if (sourceMesh->flags() & MeshFlags::HasNormals) - normal = sourceMesh->normal(vertex); - m_mesh->addVertex(sourceMesh->position(vertex), normal, sourceMesh->texcoord(vertex)); - } - } - } - // Add faces. - for (uint32_t f = 0; f < faceCount; f++) { - const uint32_t face = m_faceToSourceFaceMap[f]; - uint32_t indices[3]; - for (uint32_t i = 0; i < 3; i++) { - const uint32_t vertex = sourceMesh->vertexAt(face * 3 + i); - indices[i] = sourceVertexToVertexMap.get(vertex); - XA_DEBUG_ASSERT(indices[i] != UINT32_MAX); - } - // Don't copy flags, it doesn't matter if a face is ignored after this point. All ignored faces get their own vertex map (m_isVertexMap) ChartGroup. - // Don't hash edges if m_isVertexMap, they may be degenerate. - Mesh::AddFaceResult::Enum result = m_mesh->addFace(indices, false, !m_isVertexMap); - XA_UNUSED(result); - XA_DEBUG_ASSERT(result == Mesh::AddFaceResult::OK); - } - if (!m_isVertexMap) { - m_mesh->createColocals(); - m_mesh->createBoundaries(); - } -#if XA_DEBUG_EXPORT_OBJ_CHART_GROUPS - char filename[256]; - XA_SPRINTF(filename, sizeof(filename), "debug_mesh_%03u_chartgroup_%03u.obj", m_sourceId, m_id); - m_mesh->writeObjFile(filename); -#else - XA_UNUSED(m_id); -#endif } ~ChartGroup() { - m_mesh->~Mesh(); - XA_FREE(m_mesh); for (uint32_t i = 0; i < m_charts.size(); i++) { m_charts[i]->~Chart(); XA_FREE(m_charts[i]); } } + uint32_t segmentChartCount() const { return m_chartBasis.size(); } uint32_t chartCount() const { return m_charts.size(); } Chart *chartAt(uint32_t i) const { return m_charts[i]; } + uint32_t faceCount() const { return m_faceCount; } uint32_t paramAddedChartsCount() const { return m_paramAddedChartsCount; } uint32_t paramDeletedChartsCount() const { return m_paramDeletedChartsCount; } - bool isVertexMap() const { return m_isVertexMap; } - uint32_t mapFaceToSourceFace(uint32_t face) const { return m_faceToSourceFaceMap[face]; } - uint32_t mapVertexToSourceVertex(uint32_t i) const { return m_vertexToSourceVertexMap[i]; } - const Mesh *mesh() const { return m_mesh; } - /* - Compute charts using a simple segmentation algorithm. - - 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. - */ - void computeCharts(TaskScheduler *taskScheduler, const ChartOptions &options, segment::Atlas &atlas, ThreadLocal *chartBuffers) - { - m_chartOptions = options; - // This function may be called multiple times, so destroy existing charts. - for (uint32_t i = 0; i < m_charts.size(); i++) { - m_charts[i]->~Chart(); - XA_FREE(m_charts[i]); - } - m_charts.clear(); + void computeChartFaces(const ChartOptions &options, segment::Atlas &atlas) + { + // Create mesh from source mesh, using only the faces in this face group. + XA_PROFILE_START(createChartGroupMesh) + Mesh *mesh = createMesh(); + XA_PROFILE_END(createChartGroupMesh) + // Segment mesh into charts (arrays of faces). #if XA_DEBUG_SINGLE_CHART - Array chartFaces; - chartFaces.resize(m_mesh->faceCount()); - for (uint32_t i = 0; i < chartFaces.size(); i++) - chartFaces[i] = i; - Chart *chart = XA_NEW_ARGS(MemTag::Default, Chart, m_mesh, chartFaces, m_sourceId, m_id, 0); - m_charts.push_back(chart); + m_chartBasis.resize(1); + Fit::computeBasis(&mesh->position(0), mesh->vertexCount(), &m_chartBasis[0]); + m_chartFaces.resize(1 + mesh->faceCount()); + m_chartFaces[0] = mesh->faceCount(); + for (uint32_t i = 0; i < m_chartFaces.size(); i++) + m_chartFaces[i + 1] = i; #else XA_PROFILE_START(buildAtlas) - atlas.reset(m_sourceId, m_id, m_mesh, options); - buildAtlas(atlas, options); + atlas.reset(mesh, options); + atlas.compute(); XA_PROFILE_END(buildAtlas) - const uint32_t chartCount = atlas.chartCount(); - m_charts.resize(chartCount); - Array taskArgs; - taskArgs.resize(chartCount); - for (uint32_t i = 0; i < chartCount; i++) { - CreateChartTaskArgs &args = taskArgs[i]; - args.basis = &atlas.chartBasis(i); - args.faces = atlas.chartFaces(i); - args.mesh = m_mesh; - args.meshId = m_sourceId; - args.chartGroupId = m_id; - args.chartId = i; - args.chartBuffers = chartBuffers; - args.chart = &m_charts[i]; - } - XA_PROFILE_START(createChartMeshesReal) - TaskGroupHandle taskGroup = taskScheduler->createTaskGroup(chartCount); - for (uint32_t i = 0; i < chartCount; i++) { - Task task; - task.userData = &taskArgs[i]; - task.func = runCreateChartTask; - taskScheduler->run(taskGroup, task); - } - taskScheduler->wait(&taskGroup); - XA_PROFILE_END(createChartMeshesReal) -#endif #if XA_DEBUG_EXPORT_OBJ_CHARTS char filename[256]; - XA_SPRINTF(filename, sizeof(filename), "debug_mesh_%03u_chartgroup_%03u_charts.obj", m_sourceId, m_id); + XA_SPRINTF(filename, sizeof(filename), "debug_mesh_%03u_chartgroup_%03u_charts.obj", m_sourceMesh->id(), m_id); FILE *file; XA_FOPEN(file, filename, "w"); if (file) { - m_mesh->writeObjVertices(file); - for (uint32_t i = 0; i < chartCount; i++) { + mesh->writeObjVertices(file); + for (uint32_t i = 0; i < atlas.chartCount(); i++) { fprintf(file, "o chart_%04d\n", i); fprintf(file, "s off\n"); - const Array &faces = builder.chartFaces(i); - for (uint32_t f = 0; f < faces.size(); f++) - m_mesh->writeObjFace(file, faces[f]); + ConstArrayView faces = atlas.chartFaces(i); + for (uint32_t f = 0; f < faces.length; f++) + mesh->writeObjFace(file, faces[f]); } - m_mesh->writeObjBoundaryEges(file); - m_mesh->writeObjLinkedBoundaries(file); + mesh->writeObjBoundaryEges(file); + mesh->writeObjLinkedBoundaries(file); fclose(file); } +#endif + // Destroy mesh. + const uint32_t faceCount = mesh->faceCount(); + mesh->~Mesh(); + XA_FREE(mesh); + XA_PROFILE_START(copyChartFaces) + // Copy basis. + const uint32_t chartCount = atlas.chartCount(); + m_chartBasis.resize(chartCount); + for (uint32_t i = 0; i < chartCount; i++) + m_chartBasis[i] = atlas.chartBasis(i); + // Copy faces from segment::Atlas to m_chartFaces array with etc. encoding. + // segment::Atlas faces refer to the chart group mesh. Map them to the input mesh instead. + m_chartFaces.resize(chartCount + faceCount); + uint32_t offset = 0; + for (uint32_t i = 0; i < chartCount; i++) { + ConstArrayView faces = atlas.chartFaces(i); + m_chartFaces[offset++] = faces.length; + for (uint32_t j = 0; j < faces.length; j++) + m_chartFaces[offset++] = m_faceToSourceFaceMap[faces[j]]; + } + XA_PROFILE_END(copyChartFaces) #endif } #if XA_RECOMPUTE_CHARTS - void parameterizeCharts(TaskScheduler *taskScheduler, ParameterizeFunc func, ThreadLocal *boundaryGrid, ThreadLocal *chartBuffers, ThreadLocal *piecewiseParam) + void parameterizeCharts(TaskScheduler *taskScheduler, const ParameterizeOptions &options, ThreadLocal *boundaryGrid, ThreadLocal *chartBuffers, ThreadLocal *piecewiseParam) #else - void parameterizeCharts(TaskScheduler* taskScheduler, ParameterizeFunc func, ThreadLocal* boundaryGrid, ThreadLocal* /*chartBuffers*/) + void parameterizeCharts(TaskScheduler* taskScheduler, const ParameterizeOptions &options, ThreadLocal* boundaryGrid, ThreadLocal* chartBuffers) #endif { + // This function may be called multiple times, so destroy existing charts. + for (uint32_t i = 0; i < m_charts.size(); i++) { + m_charts[i]->~Chart(); + XA_FREE(m_charts[i]); + } m_paramAddedChartsCount = 0; - const uint32_t chartCount = m_charts.size(); - Array taskArgs; + const uint32_t chartCount = m_chartBasis.size(); + Array taskArgs; taskArgs.resize(chartCount); + taskArgs.runCtors(); // Has Array member. TaskGroupHandle taskGroup = taskScheduler->createTaskGroup(chartCount); + uint32_t offset = 0; for (uint32_t i = 0; i < chartCount; i++) { - ParameterizeChartTaskArgs &args = taskArgs[i]; - args.chart = m_charts[i]; - args.func = func; + CreateAndParameterizeChartTaskArgs &args = taskArgs[i]; + args.basis = &m_chartBasis[i]; args.boundaryGrid = boundaryGrid; + args.chart = nullptr; + args.chartGroupId = m_id; + args.chartId = i; + args.chartBuffers = chartBuffers; + const uint32_t faceCount = m_chartFaces[offset++]; + args.faces = ConstArrayView(&m_chartFaces[offset], faceCount); + offset += faceCount; + args.mesh = m_sourceMesh; + args.options = &options; +#if XA_RECOMPUTE_CHARTS + args.pp = piecewiseParam; +#endif Task task; task.userData = &args; - task.func = runParameterizeChartTask; + task.func = runCreateAndParameterizeChartTask; taskScheduler->run(taskGroup, task); } taskScheduler->wait(&taskGroup); #if XA_RECOMPUTE_CHARTS - // Find charts with invalid parameterizations. - Array invalidCharts; + // Count charts. Skip invalid ones and include new ones added by recomputing. + uint32_t newChartCount = 0; for (uint32_t i = 0; i < chartCount; i++) { - Chart *chart = m_charts[i]; - const Quality &quality = chart->quality(); - if (quality.boundaryIntersection || quality.flippedTriangleCount > 0) - invalidCharts.push_back(chart); + if (taskArgs[i].chart->isInvalid()) + newChartCount += taskArgs[i].charts.size(); + else + newChartCount++; } - if (invalidCharts.isEmpty()) - return; - // Recompute charts with invalid parameterizations. - PiecewiseParam &pp = piecewiseParam->get(); - for (uint32_t i = 0; i < invalidCharts.size(); i++) { - Chart *invalidChart = invalidCharts[i]; - // Fixing t-junctions rewrites unified mesh faces, and we need to map faces back to input mesh. So use the unmodified unified mesh. - const Mesh *invalidMesh = invalidChart->unmodifiedUnifiedMesh(); - uint32_t faceCount = 0; - if (invalidMesh) { - faceCount = invalidMesh->faceCount(); - } else { - invalidMesh = invalidChart->unifiedMesh(); - faceCount = invalidChart->initialFaceCount(); // Not invalidMesh->faceCount(). Don't want faces added by hole closing. - } - pp.reset(invalidMesh, faceCount); -#if XA_DEBUG_EXPORT_OBJ_RECOMPUTED_CHARTS - char filename[256]; - XA_SPRINTF(filename, sizeof(filename), "debug_mesh_%03u_chartgroup_%03u_recomputed_chart_%03u.obj", m_sourceId, m_id, m_paramAddedChartsCount); - FILE *file; - XA_FOPEN(file, filename, "w"); - uint32_t subChartIndex = 0; -#endif - for (;;) { - if (!pp.computeChart()) - break; - Chart *chart = XA_NEW_ARGS(MemTag::Default, Chart, chartBuffers->get(), invalidChart, invalidMesh, pp.chartFaces(), pp.texcoords(), m_mesh, m_sourceId, m_id, m_charts.size()); - m_charts.push_back(chart); -#if XA_DEBUG_EXPORT_OBJ_RECOMPUTED_CHARTS - if (file) { - for (uint32_t j = 0; j < invalidMesh->vertexCount(); j++) { - fprintf(file, "v %g %g %g\n", invalidMesh->position(j).x, invalidMesh->position(j).y, invalidMesh->position(j).z); - fprintf(file, "vt %g %g\n", pp.texcoords()[j].x, pp.texcoords()[j].y); - } - fprintf(file, "o chart%03u\n", subChartIndex); - fprintf(file, "s off\n"); - for (uint32_t f = 0; f < pp.chartFaces().length; f++) { - fprintf(file, "f "); - const uint32_t face = pp.chartFaces()[f]; - for (uint32_t j = 0; j < 3; j++) { - const uint32_t index = invalidMesh->vertexCount() * subChartIndex + invalidMesh->vertexAt(face * 3 + j) + 1; // 1-indexed - fprintf(file, "%d/%d/%c", index, index, j == 2 ? '\n' : ' '); - } - } - } - subChartIndex++; -#endif - m_paramAddedChartsCount++; + m_charts.resize(newChartCount); + // Add valid charts first. Destroy invalid ones. + uint32_t current = 0; + for (uint32_t i = 0; i < chartCount; i++) { + Chart *chart = taskArgs[i].chart; + if (chart->isInvalid()) { + chart->~Chart(); + XA_FREE(chart); + m_paramDeletedChartsCount++; + continue; } -#if XA_DEBUG_EXPORT_OBJ_RECOMPUTED_CHARTS - if (file) - fclose(file); -#endif + m_charts[current++] = chart; } - // Remove and delete the invalid charts. - for (uint32_t i = 0; i < invalidCharts.size(); i++) { - Chart *chart = invalidCharts[i]; - removeChart(chart); - chart->~Chart(); - XA_FREE(chart); - m_paramDeletedChartsCount++; + // Now add new charts. + for (uint32_t i = 0; i < chartCount; i++) { + CreateAndParameterizeChartTaskArgs &args = taskArgs[i]; + for (uint32_t j = 0; j < args.charts.size(); j++) { + m_charts[current++] = args.charts[j]; + m_paramAddedChartsCount++; + } } +#else // XA_RECOMPUTE_CHARTS + m_charts.resize(chartCount); + for (uint32_t i = 0; i < chartCount; i++) + m_charts[i] = taskArgs[i].chart; #endif // XA_RECOMPUTE_CHARTS + taskArgs.runDtors(); // Has Array member. } private: - void buildAtlas(segment::Atlas &atlas, const ChartOptions &options) + Mesh *createMesh() { - if (atlas.facesLeft() == 0) - return; - // Create initial charts greedely. - atlas.placeSeeds(options.maxThreshold * 0.5f); - if (options.maxIterations == 0) { - XA_DEBUG_ASSERT(atlas.facesLeft() == 0); - return; - } - atlas.relocateSeeds(); - atlas.resetCharts(); - // Restart process growing charts in parallel. - uint32_t iteration = 0; - for (;;) { - atlas.growCharts(options.maxThreshold); - // When charts cannot grow more: fill holes, merge charts, relocate seeds and start new iteration. - atlas.fillHoles(options.maxThreshold * 0.5f); -#if XA_MERGE_CHARTS - atlas.mergeCharts(); -#endif - if (++iteration == options.maxIterations) - break; - if (!atlas.relocateSeeds()) - break; - atlas.resetCharts(); + XA_DEBUG_ASSERT(m_faceGroup != MeshFaceGroups::kInvalid); + // Create new mesh from the source mesh, using faces that belong to this group. + m_faceToSourceFaceMap.reserve(m_sourceMeshFaceGroups->faceCount(m_faceGroup)); + for (MeshFaceGroups::Iterator it(m_sourceMeshFaceGroups, m_faceGroup); !it.isDone(); it.advance()) + m_faceToSourceFaceMap.push_back(it.face()); + // Only initial meshes has ignored faces. The only flag we care about is HasNormals. + const uint32_t faceCount = m_faceCount = m_faceToSourceFaceMap.size(); + XA_DEBUG_ASSERT(faceCount > 0); + const uint32_t approxVertexCount = min(faceCount * 3, m_sourceMesh->vertexCount()); + Mesh *mesh = XA_NEW_ARGS(MemTag::Mesh, Mesh, m_sourceMesh->epsilon(), approxVertexCount, faceCount, m_sourceMesh->flags() & MeshFlags::HasNormals); + HashMap> sourceVertexToVertexMap(MemTag::Mesh, approxVertexCount); + for (uint32_t f = 0; f < faceCount; f++) { + const uint32_t face = m_faceToSourceFaceMap[f]; + for (uint32_t i = 0; i < 3; i++) { + const uint32_t vertex = m_sourceMesh->vertexAt(face * 3 + i); + if (sourceVertexToVertexMap.get(vertex) == UINT32_MAX) { + sourceVertexToVertexMap.add(vertex); + Vector3 normal(0.0f); + if (m_sourceMesh->flags() & MeshFlags::HasNormals) + normal = m_sourceMesh->normal(vertex); + mesh->addVertex(m_sourceMesh->position(vertex), normal, m_sourceMesh->texcoord(vertex)); + } + } } - // Make sure no holes are left! - XA_DEBUG_ASSERT(atlas.facesLeft() == 0); - } - - void removeChart(const Chart *chart) - { - for (uint32_t i = 0; i < m_charts.size(); i++) { - if (m_charts[i] == chart) { - m_charts.removeAt(i); - return; + // Add faces. + for (uint32_t f = 0; f < faceCount; f++) { + const uint32_t face = m_faceToSourceFaceMap[f]; + XA_DEBUG_ASSERT(!m_sourceMesh->isFaceIgnored(face)); + uint32_t indices[3]; + for (uint32_t i = 0; i < 3; i++) { + const uint32_t vertex = m_sourceMesh->vertexAt(face * 3 + i); + indices[i] = sourceVertexToVertexMap.get(vertex); + XA_DEBUG_ASSERT(indices[i] != UINT32_MAX); } + // Don't copy flags - ignored faces aren't used by chart groups, they are handled by InvalidMeshGeometry. + Mesh::AddFaceResult::Enum result = mesh->addFace(indices); + XA_UNUSED(result); + XA_DEBUG_ASSERT(result == Mesh::AddFaceResult::OK); } + XA_PROFILE_START(createChartGroupMeshColocals) + mesh->createColocals(); + XA_PROFILE_END(createChartGroupMeshColocals) + XA_PROFILE_START(createChartGroupMeshBoundaries) + mesh->createBoundaries(); + mesh->destroyEdgeMap(); // Only needed it for createBoundaries. + XA_PROFILE_END(createChartGroupMeshBoundaries) +#if XA_DEBUG_EXPORT_OBJ_CHART_GROUPS + char filename[256]; + XA_SPRINTF(filename, sizeof(filename), "debug_mesh_%03u_chartgroup_%03u.obj", m_sourceMesh->id(), m_id); + mesh->writeObjFile(filename); +#endif + return mesh; } - uint32_t m_sourceId, m_id; - bool m_isVertexMap; - Mesh *m_mesh; + const uint32_t m_id; + const Mesh * const m_sourceMesh; + const MeshFaceGroups * const m_sourceMeshFaceGroups; + const MeshFaceGroups::Handle m_faceGroup; Array m_faceToSourceFaceMap; // List of faces of the source mesh that belong to this chart group. - Array m_vertexToSourceVertexMap; // Map vertices of the mesh to vertices of the source mesh. + Array m_chartBasis; // Copied from segment::Atlas. + Array m_chartFaces; // Copied from segment::Atlas. Encoding: etc. Array m_charts; - ChartOptions m_chartOptions; + uint32_t m_faceCount; // Set by createMesh(). Used for sorting. uint32_t m_paramAddedChartsCount; // Number of new charts added by recomputing charts with invalid parameterizations. uint32_t m_paramDeletedChartsCount; // Number of charts with invalid parameterizations that were deleted, after charts were recomputed. }; -struct CreateChartGroupTaskArgs +// References invalid faces and vertices in a mesh. +struct InvalidMeshGeometry { - uint16_t faceGroup; - uint32_t groupId; - const Mesh *mesh; - ChartGroup **chartGroup; + // Invalid faces have the face groups MeshFaceGroups::kInvalid. + void extract(const Mesh *mesh, const MeshFaceGroups *meshFaceGroups) + { + // Copy invalid faces. + m_faces.clear(); + const uint32_t meshFaceCount = mesh->faceCount(); + for (uint32_t f = 0; f < meshFaceCount; f++) { + if (meshFaceGroups->groupAt(f) == MeshFaceGroups::kInvalid) + m_faces.push_back(f); + } + // Create *unique* list of vertices of invalid faces. + const uint32_t faceCount = m_faces.size(); + m_indices.resize(faceCount * 3); + const uint32_t approxVertexCount = min(faceCount * 3, mesh->vertexCount()); + m_vertexToSourceVertexMap.clear(); + m_vertexToSourceVertexMap.reserve(approxVertexCount); + HashMap> sourceVertexToVertexMap(MemTag::Mesh, approxVertexCount); + for (uint32_t f = 0; f < faceCount; f++) { + const uint32_t face = m_faces[f]; + for (uint32_t i = 0; i < 3; i++) { + const uint32_t vertex = mesh->vertexAt(face * 3 + i); + uint32_t newVertex = sourceVertexToVertexMap.get(vertex); + if (newVertex == UINT32_MAX) { + newVertex = sourceVertexToVertexMap.add(vertex); + m_vertexToSourceVertexMap.push_back(vertex); + } + m_indices[f * 3 + i] = newVertex; + } + } + } + + ConstArrayView faces() const { return m_faces; } + ConstArrayView indices() const { return m_indices; } + ConstArrayView vertices() const { return m_vertexToSourceVertexMap; } + +private: + Array m_faces, m_indices; + Array m_vertexToSourceVertexMap; // Map face vertices to vertices of the source mesh. +}; + +struct ChartGroupComputeChartFacesTaskArgs +{ + ThreadLocal *atlas; + ChartGroup *chartGroup; + const ChartOptions *options; + Progress *progress; }; -static void runCreateChartGroupTask(void *userData) +static void runChartGroupComputeChartFacesJob(void *userData) { - XA_PROFILE_START(addMeshCreateChartGroupsThread) - auto args = (CreateChartGroupTaskArgs *)userData; - *(args->chartGroup) = XA_NEW_ARGS(MemTag::Default, ChartGroup, args->groupId, args->mesh, args->faceGroup); - XA_PROFILE_END(addMeshCreateChartGroupsThread) + auto args = (ChartGroupComputeChartFacesTaskArgs *)userData; + if (args->progress->cancel) + return; + XA_PROFILE_START(chartGroupComputeChartsThread) + args->chartGroup->computeChartFaces(*args->options, args->atlas->get()); + XA_PROFILE_END(chartGroupComputeChartsThread) } -struct ComputeChartsTaskArgs +struct MeshComputeChartFacesTaskArgs { - TaskScheduler *taskScheduler; - ChartGroup *chartGroup; + Array *chartGroups; // output + InvalidMeshGeometry *invalidMeshGeometry; // output ThreadLocal *atlas; - ThreadLocal *chartBuffers; const ChartOptions *options; Progress *progress; + const Mesh *sourceMesh; + TaskScheduler *taskScheduler; }; -static void runComputeChartsJob(void *userData) +#if XA_DEBUG_EXPORT_OBJ_FACE_GROUPS +static uint32_t s_faceGroupsCurrentVertex = 0; +#endif + +static void runMeshComputeChartFacesJob(void *userData) { - auto args = (ComputeChartsTaskArgs *)userData; + auto args = (MeshComputeChartFacesTaskArgs *)userData; if (args->progress->cancel) return; XA_PROFILE_START(computeChartsThread) - args->chartGroup->computeCharts(args->taskScheduler, *args->options, args->atlas->get(), args->chartBuffers); + // Create face groups. + XA_PROFILE_START(createFaceGroups) + MeshFaceGroups *meshFaceGroups = XA_NEW_ARGS(MemTag::Mesh, MeshFaceGroups, args->sourceMesh); + meshFaceGroups->compute(); + const uint32_t chartGroupCount = meshFaceGroups->groupCount(); + XA_PROFILE_END(createFaceGroups) + if (args->progress->cancel) + goto cleanup; +#if XA_DEBUG_EXPORT_OBJ_FACE_GROUPS + { + static std::mutex s_mutex; + std::lock_guard lock(s_mutex); + char filename[256]; + XA_SPRINTF(filename, sizeof(filename), "debug_face_groups.obj"); + FILE *file; + XA_FOPEN(file, filename, s_faceGroupsCurrentVertex == 0 ? "w" : "a"); + if (file) { + const Mesh *mesh = args->sourceMesh; + mesh->writeObjVertices(file); + // groups + uint32_t numGroups = 0; + for (uint32_t i = 0; i < mesh->faceCount(); i++) { + if (meshFaceGroups->groupAt(i) != MeshFaceGroups::kInvalid) + numGroups = max(numGroups, meshFaceGroups->groupAt(i) + 1); + } + for (uint32_t i = 0; i < numGroups; i++) { + fprintf(file, "o mesh_%03u_group_%04d\n", mesh->id(), i); + fprintf(file, "s off\n"); + for (uint32_t f = 0; f < mesh->faceCount(); f++) { + if (meshFaceGroups->groupAt(f) == i) + mesh->writeObjFace(file, f, s_faceGroupsCurrentVertex); + } + } + fprintf(file, "o mesh_%03u_group_ignored\n", mesh->id()); + fprintf(file, "s off\n"); + for (uint32_t f = 0; f < mesh->faceCount(); f++) { + if (meshFaceGroups->groupAt(f) == MeshFaceGroups::kInvalid) + mesh->writeObjFace(file, f, s_faceGroupsCurrentVertex); + } + mesh->writeObjBoundaryEges(file); + s_faceGroupsCurrentVertex += mesh->vertexCount(); + fclose(file); + } + } +#endif + // Create a chart group for each face group. + args->chartGroups->resize(chartGroupCount); + for (uint32_t i = 0; i < chartGroupCount; i++) + (*args->chartGroups)[i] = XA_NEW_ARGS(MemTag::Default, ChartGroup, i, args->sourceMesh, meshFaceGroups, MeshFaceGroups::Handle(i)); + // Extract invalid geometry via the invalid face group (MeshFaceGroups::kInvalid). + XA_PROFILE_START(extractInvalidMeshGeometry) + args->invalidMeshGeometry->extract(args->sourceMesh, meshFaceGroups); + XA_PROFILE_END(extractInvalidMeshGeometry) + // One task for each chart group - compute chart faces. + { + XA_PROFILE_START(chartGroupComputeChartsReal) + Array taskArgs; + taskArgs.resize(chartGroupCount); + for (uint32_t i = 0; i < chartGroupCount; i++) { + taskArgs[i].atlas = args->atlas; + taskArgs[i].chartGroup = (*args->chartGroups)[i]; + taskArgs[i].options = args->options; + taskArgs[i].progress = args->progress; + } + TaskGroupHandle taskGroup = args->taskScheduler->createTaskGroup(chartGroupCount); + for (uint32_t i = 0; i < chartGroupCount; i++) { + Task task; + task.userData = &taskArgs[i]; + task.func = runChartGroupComputeChartFacesJob; + args->taskScheduler->run(taskGroup, task); + } + args->taskScheduler->wait(&taskGroup); + XA_PROFILE_END(chartGroupComputeChartsReal) + } XA_PROFILE_END(computeChartsThread) args->progress->value++; args->progress->update(); +cleanup: + if (meshFaceGroups) { + meshFaceGroups->~MeshFaceGroups(); + XA_FREE(meshFaceGroups); + } } struct ParameterizeChartsTaskArgs { TaskScheduler *taskScheduler; ChartGroup *chartGroup; - ParameterizeFunc func; + const ParameterizeOptions *options; ThreadLocal *boundaryGrid; ThreadLocal *chartBuffers; #if XA_RECOMPUTE_CHARTS @@ -7426,9 +8107,9 @@ static void runParameterizeChartsJob(void *userData) return; XA_PROFILE_START(parameterizeChartsThread) #if XA_RECOMPUTE_CHARTS - args->chartGroup->parameterizeCharts(args->taskScheduler, args->func, args->boundaryGrid, args->chartBuffers, args->piecewiseParam); + args->chartGroup->parameterizeCharts(args->taskScheduler, *args->options, args->boundaryGrid, args->chartBuffers, args->piecewiseParam); #else - args->chartGroup->parameterizeCharts(args->taskScheduler, args->func, args->boundaryGrid, args->chartBuffers); + args->chartGroup->parameterizeCharts(args->taskScheduler, *args->options, args->boundaryGrid, args->chartBuffers); #endif XA_PROFILE_END(parameterizeChartsThread) args->progress->value++; @@ -7439,140 +8120,83 @@ static void runParameterizeChartsJob(void *userData) class Atlas { public: - Atlas() : m_meshCount(0), m_chartsComputed(false), m_chartsParameterized(false) {} + Atlas() : m_chartsComputed(false), m_chartsParameterized(false) {} ~Atlas() { - for (uint32_t i = 0; i < m_chartGroups.size(); i++) { - m_chartGroups[i]->~ChartGroup(); - XA_FREE(m_chartGroups[i]); + for (uint32_t i = 0; i < m_meshChartGroups.size(); i++) { + for (uint32_t j = 0; j < m_meshChartGroups[i].size(); j++) { + m_meshChartGroups[i][j]->~ChartGroup(); + XA_FREE(m_meshChartGroups[i][j]); + } } + m_meshChartGroups.runDtors(); + m_invalidMeshGeometry.runDtors(); } + uint32_t meshCount() const { return m_meshes.size(); } + const InvalidMeshGeometry &invalidMeshGeometry(uint32_t meshIndex) const { return m_invalidMeshGeometry[meshIndex]; } bool chartsComputed() const { return m_chartsComputed; } bool chartsParameterized() const { return m_chartsParameterized; } - uint32_t chartGroupCount() const { return m_chartGroups.size(); } - const ChartGroup *chartGroupAt(uint32_t index) const { return m_chartGroups[index]; } - - uint32_t chartGroupCount(uint32_t mesh) const - { - uint32_t count = 0; - for (uint32_t i = 0; i < m_chartGroups.size(); i++) { - if (m_chartGroupSourceMeshes[i] == mesh) - count++; - } - return count; - } + uint32_t chartGroupCount(uint32_t mesh) const { return m_meshChartGroups[mesh].size(); } + const ChartGroup *chartGroupAt(uint32_t mesh, uint32_t group) const { return m_meshChartGroups[mesh][group]; } - const ChartGroup *chartGroupAt(uint32_t mesh, uint32_t group) const + void addMesh(const Mesh *mesh) { - for (uint32_t c = 0; c < m_chartGroups.size(); c++) { - if (m_chartGroupSourceMeshes[c] != mesh) - continue; - if (group == 0) - return m_chartGroups[c]; - group--; - } - return nullptr; - } - - // This function is thread safe. - void addMesh(TaskScheduler *taskScheduler, const Mesh *mesh) - { - // Create one chart group per face group. - // If there's any ignored faces in the mesh, create an extra face group for that (vertex map). - // Chart group creation is slow since it copies a chunk of the source mesh, so use tasks. - Array chartGroups; - chartGroups.resize(mesh->faceGroupCount() + (mesh->ignoredFaceCount() > 0 ? 1 : 0)); - Array taskArgs; - taskArgs.resize(chartGroups.size()); - for (uint32_t g = 0; g < chartGroups.size(); g++) { - CreateChartGroupTaskArgs &args = taskArgs[g]; - args.chartGroup = &chartGroups[g]; - args.faceGroup = uint16_t(g < mesh->faceGroupCount() ? g : Mesh::kInvalidFaceGroup); - args.groupId = g; - args.mesh = mesh; - } - TaskGroupHandle taskGroup = taskScheduler->createTaskGroup(chartGroups.size()); - for (uint32_t g = 0; g < chartGroups.size(); g++) { - Task task; - task.userData = &taskArgs[g]; - task.func = runCreateChartGroupTask; - taskScheduler->run(taskGroup, task); - } - taskScheduler->wait(&taskGroup); - // Thread-safe append. - m_addMeshMutex.lock(); - for (uint32_t g = 0; g < chartGroups.size(); g++) { - m_chartGroups.push_back(chartGroups[g]); - m_chartGroupSourceMeshes.push_back(mesh->id()); - } - m_meshCount++; - m_addMeshMutex.unlock(); - } - - // Chart id/index is determined by depth-first hierarchy of mesh -> chart group -> chart. - // For chart index to be consistent here, chart groups needs to sorted by mesh index. Since addMesh is called by multithreaded tasks, order is indeterminate, so chart groups need to be explicitly sorted after all meshes are added. - void sortChartGroups() - { - Array oldChartGroups; - oldChartGroups.resize(m_chartGroups.size()); - memcpy(oldChartGroups.data(), m_chartGroups.data(), sizeof(ChartGroup *) * m_chartGroups.size()); - Array oldChartGroupSourceMeshes; - oldChartGroupSourceMeshes.resize(m_chartGroupSourceMeshes.size()); - memcpy(oldChartGroupSourceMeshes.data(), m_chartGroupSourceMeshes.data(), sizeof(uint32_t) * m_chartGroupSourceMeshes.size()); - uint32_t current = 0; - for (uint32_t i = 0; i < m_meshCount; i++) { - for (uint32_t j = 0; j < oldChartGroups.size(); j++) { - if (oldChartGroupSourceMeshes[j] == i) { - m_chartGroups[current] = oldChartGroups[j]; - m_chartGroupSourceMeshes[current] = oldChartGroupSourceMeshes[j]; - current++; - } - } - } + m_meshes.push_back(mesh); } bool computeCharts(TaskScheduler *taskScheduler, const ChartOptions &options, ProgressFunc progressFunc, void *progressUserData) { +#if XA_DEBUG_EXPORT_OBJ_PLANAR_REGIONS + segment::s_planarRegionsCurrentRegion = segment::s_planarRegionsCurrentVertex = 0; +#endif m_chartsComputed = false; m_chartsParameterized = false; - // Ignore vertex maps. - uint32_t chartGroupCount = 0; - for (uint32_t i = 0; i < m_chartGroups.size(); i++) { - if (!m_chartGroups[i]->isVertexMap()) - chartGroupCount++; + // Clear chart groups, since this function may be called multiple times. + if (!m_meshChartGroups.isEmpty()) { + for (uint32_t i = 0; i < m_meshChartGroups.size(); i++) { + for (uint32_t j = 0; j < m_meshChartGroups[i].size(); j++) { + m_meshChartGroups[i][j]->~ChartGroup(); + XA_FREE(m_meshChartGroups[i][j]); + } + m_meshChartGroups[i].clear(); + } + XA_ASSERT(m_meshChartGroups.size() == m_meshes.size()); // The number of meshes shouldn't have changed. } - Progress progress(ProgressCategory::ComputeCharts, progressFunc, progressUserData, chartGroupCount); + m_meshChartGroups.resize(m_meshes.size()); + m_meshChartGroups.runCtors(); + m_invalidMeshGeometry.resize(m_meshes.size()); + m_invalidMeshGeometry.runCtors(); + // One task per mesh. + const uint32_t meshCount = m_meshes.size(); + Progress progress(ProgressCategory::ComputeCharts, progressFunc, progressUserData, meshCount); ThreadLocal atlas; - ThreadLocal chartBuffers; - Array taskArgs; - taskArgs.reserve(chartGroupCount); - for (uint32_t i = 0; i < m_chartGroups.size(); i++) { - if (!m_chartGroups[i]->isVertexMap()) { - ComputeChartsTaskArgs args; - args.taskScheduler = taskScheduler; - args.chartGroup = m_chartGroups[i]; - args.atlas = &atlas; - args.chartBuffers = &chartBuffers; - args.options = &options; - args.progress = &progress; - taskArgs.push_back(args); - } - } - // Sort chart groups by mesh indexCount. - m_chartGroupsRadix = RadixSort(); - Array chartGroupSortData; - chartGroupSortData.resize(chartGroupCount); - for (uint32_t i = 0; i < chartGroupCount; i++) - chartGroupSortData[i] = (float)taskArgs[i].chartGroup->mesh()->indexCount(); - m_chartGroupsRadix.sort(chartGroupSortData); - // Larger chart group meshes are added first to reduce the chance of thread starvation. - TaskGroupHandle taskGroup = taskScheduler->createTaskGroup(chartGroupCount); - for (uint32_t i = 0; i < chartGroupCount; i++) { + Array taskArgs; + taskArgs.resize(meshCount); + for (uint32_t i = 0; i < meshCount; i++) { + MeshComputeChartFacesTaskArgs &args = taskArgs[i]; + args.atlas = &atlas; + args.chartGroups = &m_meshChartGroups[i]; + args.invalidMeshGeometry = &m_invalidMeshGeometry[i]; + args.options = &options; + args.progress = &progress; + args.sourceMesh = m_meshes[i]; + args.taskScheduler = taskScheduler; + } + // Sort meshes by indexCount. + Array meshSortData; + meshSortData.resize(meshCount); + for (uint32_t i = 0; i < meshCount; i++) + meshSortData[i] = (float)m_meshes[i]->indexCount(); + RadixSort meshSort; + meshSort.sort(meshSortData); + // Larger meshes are added first to reduce the chance of thread starvation. + TaskGroupHandle taskGroup = taskScheduler->createTaskGroup(meshCount); + for (uint32_t i = 0; i < meshCount; i++) { Task task; - task.userData = &taskArgs[m_chartGroupsRadix.ranks()[chartGroupCount - i - 1]]; - task.func = runComputeChartsJob; + task.userData = &taskArgs[meshSort.ranks()[meshCount - i - 1]]; + task.func = runMeshComputeChartFacesJob; taskScheduler->run(taskGroup, task); } taskScheduler->wait(&taskGroup); @@ -7582,15 +8206,12 @@ public: return true; } - bool parameterizeCharts(TaskScheduler *taskScheduler, ParameterizeFunc func, ProgressFunc progressFunc, void *progressUserData) + bool parameterizeCharts(TaskScheduler *taskScheduler, const ParameterizeOptions &options, ProgressFunc progressFunc, void *progressUserData) { m_chartsParameterized = false; - // Ignore vertex maps. uint32_t chartGroupCount = 0; - for (uint32_t i = 0; i < m_chartGroups.size(); i++) { - if (!m_chartGroups[i]->isVertexMap()) - chartGroupCount++; - } + for (uint32_t i = 0; i < m_meshChartGroups.size(); i++) + chartGroupCount += m_meshChartGroups[i].size(); Progress progress(ProgressCategory::ParameterizeCharts, progressFunc, progressUserData, chartGroupCount); ThreadLocal boundaryGrid; // For Quality boundary intersection. ThreadLocal chartBuffers; @@ -7598,27 +8219,45 @@ public: ThreadLocal piecewiseParam; #endif Array taskArgs; - taskArgs.reserve(chartGroupCount); - for (uint32_t i = 0; i < m_chartGroups.size(); i++) { - if (!m_chartGroups[i]->isVertexMap()) { - ParameterizeChartsTaskArgs args; - args.taskScheduler = taskScheduler; - args.chartGroup = m_chartGroups[i]; - args.func = func; - args.boundaryGrid = &boundaryGrid; - args.chartBuffers = &chartBuffers; + taskArgs.resize(chartGroupCount); + { + uint32_t k = 0; + for (uint32_t i = 0; i < m_meshChartGroups.size(); i++) { + const uint32_t count = m_meshChartGroups[i].size(); + for (uint32_t j = 0; j < count; j++) { + ParameterizeChartsTaskArgs &args = taskArgs[k]; + args.taskScheduler = taskScheduler; + args.chartGroup = m_meshChartGroups[i][j]; + args.options = &options; + args.boundaryGrid = &boundaryGrid; + args.chartBuffers = &chartBuffers; #if XA_RECOMPUTE_CHARTS - args.piecewiseParam = &piecewiseParam; + args.piecewiseParam = &piecewiseParam; #endif - args.progress = &progress; - taskArgs.push_back(args); + args.progress = &progress; + k++; + } } } - // Larger chart group meshes are added first to reduce the chance of thread starvation. + // Sort chart groups by face count. + Array chartGroupSortData; + chartGroupSortData.resize(chartGroupCount); + { + uint32_t k = 0; + for (uint32_t i = 0; i < m_meshChartGroups.size(); i++) { + const uint32_t count = m_meshChartGroups[i].size(); + for (uint32_t j = 0; j < count; j++) { + chartGroupSortData[k++] = (float)m_meshChartGroups[i][j]->faceCount(); + } + } + } + RadixSort chartGroupSort; + chartGroupSort.sort(chartGroupSortData); + // Larger chart groups are added first to reduce the chance of thread starvation. TaskGroupHandle taskGroup = taskScheduler->createTaskGroup(chartGroupCount); for (uint32_t i = 0; i < chartGroupCount; i++) { Task task; - task.userData = &taskArgs[m_chartGroupsRadix.ranks()[chartGroupCount - i - 1]]; + task.userData = &taskArgs[chartGroupSort.ranks()[chartGroupCount - i - 1]]; task.func = runParameterizeChartsJob; taskScheduler->run(taskGroup, task); } @@ -7630,13 +8269,11 @@ public: } private: - std::mutex m_addMeshMutex; - uint32_t m_meshCount; + Array m_meshes; + Array m_invalidMeshGeometry; // 1 per mesh. + Array > m_meshChartGroups; bool m_chartsComputed; bool m_chartsParameterized; - Array m_chartGroups; - RadixSort m_chartGroupsRadix; // By mesh indexCount. - Array m_chartGroupSourceMeshes; }; } // namespace param @@ -7777,7 +8414,7 @@ static void runAddChartTask(void *userData) auto args = (AddChartTaskArgs *)userData; param::Chart *paramChart = args->paramChart; XA_PROFILE_START(packChartsAddChartsRestoreTexcoords) - paramChart->transferParameterization(); + paramChart->restoreTexcoords(); XA_PROFILE_END(packChartsAddChartsRestoreTexcoords) Mesh *mesh = paramChart->mesh(); Chart *chart = args->chart = XA_NEW(MemTag::Default, Chart); @@ -7842,12 +8479,12 @@ struct Atlas { // Count charts. uint32_t chartCount = 0; - const uint32_t chartGroupsCount = paramAtlas->chartGroupCount(); - for (uint32_t i = 0; i < chartGroupsCount; i++) { - const param::ChartGroup *chartGroup = paramAtlas->chartGroupAt(i); - if (chartGroup->isVertexMap()) - continue; - chartCount += chartGroup->chartCount(); + for (uint32_t i = 0; i < paramAtlas->meshCount(); i++) { + const uint32_t chartGroupsCount = paramAtlas->chartGroupCount(i); + for (uint32_t j = 0; j < chartGroupsCount; j++) { + const param::ChartGroup *chartGroup = paramAtlas->chartGroupAt(i, j); + chartCount += chartGroup->chartCount(); + } } if (chartCount == 0) return; @@ -7857,20 +8494,21 @@ struct Atlas TaskGroupHandle taskGroup = taskScheduler->createTaskGroup(chartCount); uint32_t chartIndex = 0; ThreadLocal boundingBox; - for (uint32_t i = 0; i < chartGroupsCount; i++) { - const param::ChartGroup *chartGroup = paramAtlas->chartGroupAt(i); - if (chartGroup->isVertexMap()) - continue; - const uint32_t count = chartGroup->chartCount(); - for (uint32_t j = 0; j < count; j++) { - AddChartTaskArgs &args = taskArgs[chartIndex]; - args.boundingBox = &boundingBox; - args.paramChart = chartGroup->chartAt(j); - Task task; - task.userData = &taskArgs[chartIndex]; - task.func = runAddChartTask; - taskScheduler->run(taskGroup, task); - chartIndex++; + for (uint32_t i = 0; i < paramAtlas->meshCount(); i++) { + const uint32_t chartGroupsCount = paramAtlas->chartGroupCount(i); + for (uint32_t j = 0; j < chartGroupsCount; j++) { + const param::ChartGroup *chartGroup = paramAtlas->chartGroupAt(i, j); + const uint32_t count = chartGroup->chartCount(); + for (uint32_t k = 0; k < count; k++) { + AddChartTaskArgs &args = taskArgs[chartIndex]; + args.boundingBox = &boundingBox; + args.paramChart = chartGroup->chartAt(k); + Task task; + task.userData = &taskArgs[chartIndex]; + task.func = runAddChartTask; + taskScheduler->run(taskGroup, task); + chartIndex++; + } } } taskScheduler->wait(&taskGroup); @@ -8070,7 +8708,6 @@ struct Atlas maxChartPerimeter = max(maxChartPerimeter, chartOrderArray[c]); } // Sort charts by perimeter. - m_radix = RadixSort(); m_radix.sort(chartOrderArray); const uint32_t *ranks = m_radix.ranks(); // Divide chart perimeter range into buckets. @@ -8265,10 +8902,8 @@ struct Atlas } texcoord.x = best_x + t.x; texcoord.y = best_y + t.y; - if (!options.blockAlign) { - texcoord.x -= (float)options.padding; - texcoord.y -= (float)options.padding; - } + texcoord.x -= (float)options.padding; + texcoord.y -= (float)options.padding; XA_ASSERT(texcoord.x >= 0 && texcoord.y >= 0); XA_ASSERT(isFinite(texcoord.x) && isFinite(texcoord.y)); } @@ -8281,21 +8916,12 @@ struct Atlas } } } - if (options.blockAlign) { - if (maxResolution == 0) { - m_width = max(0, atlasSizes[0].x); - m_height = max(0, atlasSizes[0].y); - } else { - m_width = m_height = maxResolution; - } + // Remove padding from outer edges. + if (maxResolution == 0) { + m_width = max(0, atlasSizes[0].x - (int)options.padding * 2); + m_height = max(0, atlasSizes[0].y - (int)options.padding * 2); } else { - // Remove padding from outer edges. - if (maxResolution == 0) { - m_width = max(0, atlasSizes[0].x - (int)options.padding * 2); - m_height = max(0, atlasSizes[0].y - (int)options.padding * 2); - } else { - m_width = m_height = maxResolution - (int)options.padding * 2; - } + m_width = m_height = maxResolution - (int)options.padding * 2; } XA_PRINT(" %dx%d resolution\n", m_width, m_height); m_utilization.resize(m_bitImages.size()); @@ -8557,13 +9183,13 @@ private: struct Context { Atlas atlas; - uint32_t meshCount = 0; internal::Progress *addMeshProgress = nullptr; internal::TaskGroupHandle addMeshTaskGroup; internal::param::Atlas paramAtlas; ProgressFunc progressFunc = nullptr; void *progressUserData = nullptr; internal::TaskScheduler *taskScheduler; + internal::Array meshes; internal::Array uvMeshes; internal::Array uvMeshInstances; }; @@ -8582,19 +9208,19 @@ static void DestroyOutputMeshes(Context *ctx) return; for (int i = 0; i < (int)ctx->atlas.meshCount; i++) { Mesh &mesh = ctx->atlas.meshes[i]; - for (uint32_t j = 0; j < mesh.chartCount; j++) { - if (mesh.chartArray[j].faceArray) - XA_FREE(mesh.chartArray[j].faceArray); - } - if (mesh.chartArray) + if (mesh.chartArray) { + for (uint32_t j = 0; j < mesh.chartCount; j++) { + if (mesh.chartArray[j].faceArray) + XA_FREE(mesh.chartArray[j].faceArray); + } XA_FREE(mesh.chartArray); + } if (mesh.vertexArray) XA_FREE(mesh.vertexArray); if (mesh.indexArray) XA_FREE(mesh.indexArray); } - if (ctx->atlas.meshes) - XA_FREE(ctx->atlas.meshes); + XA_FREE(ctx->atlas.meshes); ctx->atlas.meshes = nullptr; } @@ -8613,6 +9239,11 @@ void Destroy(Atlas *atlas) } ctx->taskScheduler->~TaskScheduler(); XA_FREE(ctx->taskScheduler); + for (uint32_t i = 0; i < ctx->meshes.size(); i++) { + internal::Mesh *mesh = ctx->meshes[i]; + mesh->~Mesh(); + XA_FREE(mesh); + } for (uint32_t i = 0; i < ctx->uvMeshes.size(); i++) { internal::UvMesh *mesh = ctx->uvMeshes[i]; for (uint32_t j = 0; j < mesh->charts.size(); j++) { @@ -8653,58 +9284,11 @@ static void runAddMeshTask(void *userData) mesh->createColocals(); XA_PROFILE_END(addMeshCreateColocals) } - if (progress->cancel) - goto cleanup; - { - XA_PROFILE_START(addMeshCreateFaceGroups) - mesh->createFaceGroups(); - XA_PROFILE_END(addMeshCreateFaceGroups) - } - if (progress->cancel) - goto cleanup; -#if XA_DEBUG_EXPORT_OBJ_SOURCE_MESHES - char filename[256]; - XA_SPRINTF(filename, sizeof(filename), "debug_mesh_%03u.obj", mesh->id()); - FILE *file; - XA_FOPEN(file, filename, "w"); - if (file) { - mesh->writeObjVertices(file); - // groups - uint32_t numGroups = 0; - for (uint32_t i = 0; i < mesh->faceCount(); i++) { - if (mesh->faceGroupAt(i) != Mesh::kInvalidFaceGroup) - numGroups = internal::max(numGroups, mesh->faceGroupAt(i) + 1); - } - for (uint32_t i = 0; i < numGroups; i++) { - fprintf(file, "o group_%04d\n", i); - fprintf(file, "s off\n"); - for (uint32_t f = 0; f < mesh->faceCount(); f++) { - if (mesh->faceGroupAt(f) == i) - mesh->writeObjFace(file, f); - } - } - fprintf(file, "o group_ignored\n"); - fprintf(file, "s off\n"); - for (uint32_t f = 0; f < mesh->faceCount(); f++) { - if (mesh->faceGroupAt(f) == Mesh::kInvalidFaceGroup) - mesh->writeObjFace(file, f); - } - mesh->writeObjBoundaryEges(file); - fclose(file); - } -#endif - { - XA_PROFILE_START(addMeshCreateChartGroupsReal) - args->ctx->paramAtlas.addMesh(args->ctx->taskScheduler, mesh); // addMesh is thread safe - XA_PROFILE_END(addMeshCreateChartGroupsReal) - } if (progress->cancel) goto cleanup; progress->value++; progress->update(); cleanup: - mesh->~Mesh(); - XA_FREE(mesh); args->~AddMeshTaskArgs(); XA_FREE(args); XA_PROFILE_END(addMeshThread) @@ -8752,7 +9336,7 @@ AddMeshError::Enum AddMesh(Atlas *atlas, const MeshDecl &meshDecl, uint32_t mesh return AddMeshError::Error; } #if XA_PROFILE - if (ctx->meshCount == 0) + if (ctx->meshes.isEmpty()) internal::s_profile.addMeshReal = clock(); #endif // Don't know how many times AddMesh will be called, so progress needs to adjusted each time. @@ -8760,12 +9344,12 @@ AddMeshError::Enum AddMesh(Atlas *atlas, const MeshDecl &meshDecl, uint32_t mesh ctx->addMeshProgress = XA_NEW_ARGS(internal::MemTag::Default, internal::Progress, ProgressCategory::AddMesh, ctx->progressFunc, ctx->progressUserData, 1); } else { - ctx->addMeshProgress->setMaxValue(internal::max(ctx->meshCount + 1, meshCountHint)); + ctx->addMeshProgress->setMaxValue(internal::max(ctx->meshes.size() + 1, meshCountHint)); } XA_PROFILE_START(addMeshCopyData) const bool hasIndices = meshDecl.indexCount > 0; const uint32_t indexCount = hasIndices ? meshDecl.indexCount : meshDecl.vertexCount; - XA_PRINT("Adding mesh %d: %u vertices, %u triangles\n", ctx->meshCount, meshDecl.vertexCount, indexCount / 3); + XA_PRINT("Adding mesh %d: %u vertices, %u triangles\n", ctx->meshes.size(), meshDecl.vertexCount, indexCount / 3); // Expecting triangle faces. if ((indexCount % 3) != 0) return AddMeshError::InvalidIndexCount; @@ -8777,10 +9361,10 @@ AddMeshError::Enum AddMesh(Atlas *atlas, const MeshDecl &meshDecl, uint32_t mesh return AddMeshError::IndexOutOfRange; } } - uint32_t meshFlags = internal::MeshFlags::HasFaceGroups | internal::MeshFlags::HasIgnoredFaces; + uint32_t meshFlags = internal::MeshFlags::HasIgnoredFaces; if (meshDecl.vertexNormalData) meshFlags |= internal::MeshFlags::HasNormals; - internal::Mesh *mesh = XA_NEW_ARGS(internal::MemTag::Mesh, internal::Mesh, meshDecl.epsilon, meshDecl.vertexCount, indexCount / 3, meshFlags, ctx->meshCount); + internal::Mesh *mesh = XA_NEW_ARGS(internal::MemTag::Mesh, internal::Mesh, meshDecl.epsilon, meshDecl.vertexCount, indexCount / 3, meshFlags, ctx->meshes.size()); for (uint32_t i = 0; i < meshDecl.vertexCount; i++) { internal::Vector3 normal(0.0f); internal::Vector2 texcoord(0.0f); @@ -8790,6 +9374,8 @@ AddMeshError::Enum AddMesh(Atlas *atlas, const MeshDecl &meshDecl, uint32_t mesh texcoord = DecodeUv(meshDecl, i); mesh->addVertex(DecodePosition(meshDecl, i), normal, texcoord); } + const uint32_t kMaxWarnings = 50; + uint32_t warningCount = 0; for (uint32_t i = 0; i < indexCount / 3; i++) { uint32_t tri[3]; for (int j = 0; j < 3; j++) @@ -8801,14 +9387,16 @@ AddMeshError::Enum AddMesh(Atlas *atlas, const MeshDecl &meshDecl, uint32_t mesh const uint32_t index2 = tri[(j + 1) % 3]; if (index1 == index2) { ignore = true; - XA_PRINT(" Degenerate edge: index %d, index %d\n", index1, index2); + if (++warningCount <= kMaxWarnings) + XA_PRINT(" Degenerate edge: index %d, index %d\n", index1, index2); break; } const internal::Vector3 &pos1 = mesh->position(index1); const internal::Vector3 &pos2 = mesh->position(index2); if (internal::length(pos2 - pos1) <= 0.0f) { ignore = true; - XA_PRINT(" Zero length edge: index %d position (%g %g %g), index %d position (%g %g %g)\n", index1, pos1.x, pos1.y, pos1.z, index2, pos2.x, pos2.y, pos2.z); + if (++warningCount <= kMaxWarnings) + XA_PRINT(" Zero length edge: index %d position (%g %g %g), index %d position (%g %g %g)\n", index1, pos1.x, pos1.y, pos1.z, index2, pos2.x, pos2.y, pos2.z); break; } } @@ -8817,14 +9405,16 @@ AddMeshError::Enum AddMesh(Atlas *atlas, const MeshDecl &meshDecl, uint32_t mesh for (int j = 0; j < 3; j++) { const internal::Vector3 &pos = mesh->position(tri[j]); if (internal::isNan(pos.x) || internal::isNan(pos.y) || internal::isNan(pos.z)) { - XA_PRINT(" NAN position in face: %d\n", i); + if (++warningCount <= kMaxWarnings) + XA_PRINT(" NAN position in face: %d\n", i); ignore = true; break; } if (meshDecl.vertexNormalData) { const internal::Vector3 &normal = mesh->normal(tri[j]); if (internal::isNan(normal.x) || internal::isNan(normal.y) || internal::isNan(normal.z)) { - XA_PRINT(" NAN normal in face: %d\n", i); + if (++warningCount <= kMaxWarnings) + XA_PRINT(" NAN normal in face: %d\n", i); ignore = true; break; } @@ -8832,7 +9422,8 @@ AddMeshError::Enum AddMesh(Atlas *atlas, const MeshDecl &meshDecl, uint32_t mesh if (meshDecl.vertexUvData) { const internal::Vector2 &uv = mesh->texcoord(tri[j]); if (internal::isNan(uv.x) || internal::isNan(uv.y)) { - XA_PRINT(" NAN texture coordinate in face: %d\n", i); + if (++warningCount <= kMaxWarnings) + XA_PRINT(" NAN texture coordinate in face: %d\n", i); ignore = true; break; } @@ -8848,20 +9439,26 @@ AddMeshError::Enum AddMesh(Atlas *atlas, const MeshDecl &meshDecl, uint32_t mesh area = internal::length(internal::cross(b - a, c - a)) * 0.5f; if (area <= internal::kAreaEpsilon) { ignore = true; - XA_PRINT(" Zero area face: %d, indices (%d %d %d), area is %f\n", i, tri[0], tri[1], tri[2], area); + if (++warningCount <= kMaxWarnings) + XA_PRINT(" Zero area face: %d, indices (%d %d %d), area is %f\n", i, tri[0], tri[1], tri[2], area); } } if (!ignore) { if (internal::equal(a, b, meshDecl.epsilon) || internal::equal(a, c, meshDecl.epsilon) || internal::equal(b, c, meshDecl.epsilon)) { ignore = true; - XA_PRINT(" Degenerate face: %d, area is %f\n", i, area); + if (++warningCount <= kMaxWarnings) + XA_PRINT(" Degenerate face: %d, area is %f\n", i, area); } } if (meshDecl.faceIgnoreData && meshDecl.faceIgnoreData[i]) ignore = true; mesh->addFace(tri[0], tri[1], tri[2], ignore); } + if (warningCount > kMaxWarnings) + XA_PRINT(" %u additional warnings truncated\n", warningCount - kMaxWarnings); XA_PROFILE_END(addMeshCopyData) + ctx->meshes.push_back(mesh); + ctx->paramAtlas.addMesh(mesh); if (ctx->addMeshTaskGroup.value == UINT32_MAX) ctx->addMeshTaskGroup = ctx->taskScheduler->createTaskGroup(); AddMeshTaskArgs *taskArgs = XA_NEW(internal::MemTag::Default, AddMeshTaskArgs); // The task frees this. @@ -8871,7 +9468,6 @@ AddMeshError::Enum AddMesh(Atlas *atlas, const MeshDecl &meshDecl, uint32_t mesh task.userData = taskArgs; task.func = runAddMeshTask; ctx->taskScheduler->run(ctx->addMeshTaskGroup, task); - ctx->meshCount++; return AddMeshError::Success; } @@ -8889,19 +9485,21 @@ void AddMeshJoin(Atlas *atlas) ctx->addMeshProgress->~Progress(); XA_FREE(ctx->addMeshProgress); ctx->addMeshProgress = nullptr; - ctx->paramAtlas.sortChartGroups(); #if XA_PROFILE - XA_PRINT("Added %u meshes\n", ctx->meshCount); + XA_PRINT("Added %u meshes\n", ctx->meshes.size()); internal::s_profile.addMeshReal = clock() - internal::s_profile.addMeshReal; #endif XA_PROFILE_PRINT_AND_RESET(" Total (real): ", addMeshReal) XA_PROFILE_PRINT_AND_RESET(" Copy data: ", addMeshCopyData) XA_PROFILE_PRINT_AND_RESET(" Total (thread): ", addMeshThread) XA_PROFILE_PRINT_AND_RESET(" Create colocals: ", addMeshCreateColocals) - XA_PROFILE_PRINT_AND_RESET(" Create face groups: ", addMeshCreateFaceGroups) - XA_PROFILE_PRINT_AND_RESET(" Create chart groups (real): ", addMeshCreateChartGroupsReal) - XA_PROFILE_PRINT_AND_RESET(" Create chart groups (thread): ", addMeshCreateChartGroupsThread) +#if XA_PROFILE_ALLOC + XA_PROFILE_PRINT_AND_RESET(" Alloc: ", alloc) +#endif XA_PRINT_MEM_USAGE +#if XA_DEBUG_EXPORT_OBJ_FACE_GROUPS + internal::param::s_faceGroupsCurrentVertex = 0; +#endif } struct EdgeKey @@ -8923,7 +9521,7 @@ AddMeshError::Enum AddUvMesh(Atlas *atlas, const UvMeshDecl &decl) return AddMeshError::Error; } Context *ctx = (Context *)atlas; - if (ctx->meshCount > 0) { + if (!ctx->meshes.isEmpty()) { XA_PRINT_WARNING("AddUvMesh: Meshes and UV meshes cannot be added to the same atlas.\n"); return AddMeshError::Error; } @@ -9026,7 +9624,7 @@ AddMeshError::Enum AddUvMesh(Atlas *atlas, const UvMeshDecl &decl) return AddMeshError::Success; } -void ComputeCharts(Atlas *atlas, ChartOptions chartOptions) +void ComputeCharts(Atlas *atlas, ChartOptions options) { if (!atlas) { XA_PRINT_WARNING("ComputeCharts: atlas is null.\n"); @@ -9038,69 +9636,65 @@ void ComputeCharts(Atlas *atlas, ChartOptions chartOptions) return; } AddMeshJoin(atlas); - if (ctx->meshCount == 0) { + if (ctx->meshes.isEmpty()) { XA_PRINT_WARNING("ComputeCharts: No meshes. Call AddMesh first.\n"); return; } XA_PRINT("Computing charts\n"); - uint32_t chartCount = 0, chartsWithHolesCount = 0, holesCount = 0, chartsWithTJunctionsCount = 0, tJunctionsCount = 0; XA_PROFILE_START(computeChartsReal) - if (!ctx->paramAtlas.computeCharts(ctx->taskScheduler, chartOptions, ctx->progressFunc, ctx->progressUserData)) { + if (!ctx->paramAtlas.computeCharts(ctx->taskScheduler, options, ctx->progressFunc, ctx->progressUserData)) { XA_PRINT(" Cancelled by user\n"); return; } XA_PROFILE_END(computeChartsReal) - // Count charts and print warnings. - for (uint32_t i = 0; i < ctx->meshCount; i++) { + // Count charts. + uint32_t chartCount = 0; + const uint32_t meshCount = ctx->meshes.size(); + for (uint32_t i = 0; i < meshCount; i++) { for (uint32_t j = 0; j < ctx->paramAtlas.chartGroupCount(i); j++) { const internal::param::ChartGroup *chartGroup = ctx->paramAtlas.chartGroupAt(i, j); - if (chartGroup->isVertexMap()) - continue; - for (uint32_t k = 0; k < chartGroup->chartCount(); k++) { - const internal::param::Chart *chart = chartGroup->chartAt(k); -#if XA_PRINT_CHART_WARNINGS - if (chart->warningFlags() & internal::param::ChartWarningFlags::CloseHolesFailed) - XA_PRINT_WARNING(" Chart %u (mesh %u, group %u, id %u): failed to close holes\n", chartCount, i, j, k); - if (chart->warningFlags() & internal::param::ChartWarningFlags::FixTJunctionsDuplicatedEdge) - XA_PRINT_WARNING(" Chart %u (mesh %u, group %u, id %u): fixing t-junctions created non-manifold geometry\n", chartCount, i, j, k); - if (chart->warningFlags() & internal::param::ChartWarningFlags::FixTJunctionsFailed) - XA_PRINT_WARNING(" Chart %u (mesh %u, group %u, id %u): fixing t-junctions failed\n", chartCount, i, j, k); - if (chart->warningFlags() & internal::param::ChartWarningFlags::TriangulateDuplicatedEdge) - XA_PRINT_WARNING(" Chart %u (mesh %u, group %u, id %u): triangulation created non-manifold geometry\n", chartCount, i, j, k); -#endif - holesCount += chart->closedHolesCount(); - if (chart->closedHolesCount() > 0) - chartsWithHolesCount++; - tJunctionsCount += chart->fixedTJunctionsCount(); - if (chart->fixedTJunctionsCount() > 0) - chartsWithTJunctionsCount++; - chartCount++; - } + chartCount += chartGroup->segmentChartCount(); } } - if (holesCount > 0) - XA_PRINT(" Closed %u holes in %u charts\n", holesCount, chartsWithHolesCount); - if (tJunctionsCount > 0) - XA_PRINT(" Fixed %u t-junctions in %u charts\n", tJunctionsCount, chartsWithTJunctionsCount); XA_PRINT(" %u charts\n", chartCount); +#if XA_PROFILE + XA_PRINT(" Chart groups\n"); + uint32_t chartGroupCount = 0; + for (uint32_t i = 0; i < meshCount; i++) { + XA_PRINT(" Mesh %u: %u chart groups\n", i, ctx->paramAtlas.chartGroupCount(i)); + chartGroupCount += ctx->paramAtlas.chartGroupCount(i); + } + XA_PRINT(" %u total\n", chartGroupCount); +#endif XA_PROFILE_PRINT_AND_RESET(" Total (real): ", computeChartsReal) XA_PROFILE_PRINT_AND_RESET(" Total (thread): ", computeChartsThread) - XA_PROFILE_PRINT_AND_RESET(" Build atlas: ", buildAtlas) - XA_PROFILE_PRINT_AND_RESET(" Init: ", buildAtlasInit) - XA_PROFILE_PRINT_AND_RESET(" Place seeds: ", buildAtlasPlaceSeeds) - XA_PROFILE_PRINT_AND_RESET(" Relocate seeds: ", buildAtlasRelocateSeeds) - XA_PROFILE_PRINT_AND_RESET(" Reset charts: ", buildAtlasResetCharts) - XA_PROFILE_PRINT_AND_RESET(" Grow charts: ", buildAtlasGrowCharts) - XA_PROFILE_PRINT_AND_RESET(" Merge charts: ", buildAtlasMergeCharts) - XA_PROFILE_PRINT_AND_RESET(" Fill holes: ", buildAtlasFillHoles) - XA_PROFILE_PRINT_AND_RESET(" Create chart meshes (real): ", createChartMeshesReal) - XA_PROFILE_PRINT_AND_RESET(" Create chart meshes (thread): ", createChartMeshesThread) - XA_PROFILE_PRINT_AND_RESET(" Fix t-junctions: ", fixChartMeshTJunctions) - XA_PROFILE_PRINT_AND_RESET(" Close holes: ", closeChartMeshHoles) + XA_PROFILE_PRINT_AND_RESET(" Create face groups: ", createFaceGroups) + XA_PROFILE_PRINT_AND_RESET(" Extract invalid mesh geometry: ", extractInvalidMeshGeometry) + XA_PROFILE_PRINT_AND_RESET(" Chart group compute charts (real): ", chartGroupComputeChartsReal) + XA_PROFILE_PRINT_AND_RESET(" Chart group compute charts (thread): ", chartGroupComputeChartsThread) + XA_PROFILE_PRINT_AND_RESET(" Create chart group mesh: ", createChartGroupMesh) + XA_PROFILE_PRINT_AND_RESET(" Create colocals: ", createChartGroupMeshColocals) + XA_PROFILE_PRINT_AND_RESET(" Create boundaries: ", createChartGroupMeshBoundaries) + XA_PROFILE_PRINT_AND_RESET(" Build atlas: ", buildAtlas) + XA_PROFILE_PRINT_AND_RESET(" Init: ", buildAtlasInit) + XA_PROFILE_PRINT_AND_RESET(" Planar charts: ", planarCharts) + XA_PROFILE_PRINT_AND_RESET(" Clustered charts: ", clusteredCharts) + XA_PROFILE_PRINT_AND_RESET(" Place seeds: ", clusteredChartsPlaceSeeds) + XA_PROFILE_PRINT_AND_RESET(" Boundary intersection: ", clusteredChartsPlaceSeedsBoundaryIntersection) + XA_PROFILE_PRINT_AND_RESET(" Relocate seeds: ", clusteredChartsRelocateSeeds) + XA_PROFILE_PRINT_AND_RESET(" Reset: ", clusteredChartsReset) + XA_PROFILE_PRINT_AND_RESET(" Grow: ", clusteredChartsGrow) + XA_PROFILE_PRINT_AND_RESET(" Boundary intersection: ", clusteredChartsGrowBoundaryIntersection) + XA_PROFILE_PRINT_AND_RESET(" Merge: ", clusteredChartsMerge) + XA_PROFILE_PRINT_AND_RESET(" Fill holes: ", clusteredChartsFillHoles) + XA_PROFILE_PRINT_AND_RESET(" Copy chart faces: ", copyChartFaces) +#if XA_PROFILE_ALLOC + XA_PROFILE_PRINT_AND_RESET(" Alloc: ", alloc) +#endif XA_PRINT_MEM_USAGE } -void ParameterizeCharts(Atlas *atlas, ParameterizeFunc func) +void ParameterizeCharts(Atlas *atlas, ParameterizeOptions options) { if (!atlas) { XA_PRINT_WARNING("ParameterizeCharts: atlas is null.\n"); @@ -9130,19 +9724,34 @@ void ParameterizeCharts(Atlas *atlas, ParameterizeFunc func) DestroyOutputMeshes(ctx); XA_PRINT("Parameterizing charts\n"); XA_PROFILE_START(parameterizeChartsReal) - if (!ctx->paramAtlas.parameterizeCharts(ctx->taskScheduler, func, ctx->progressFunc, ctx->progressUserData)) { + if (!ctx->paramAtlas.parameterizeCharts(ctx->taskScheduler, options, ctx->progressFunc, ctx->progressUserData)) { XA_PRINT(" Cancelled by user\n"); return; } XA_PROFILE_END(parameterizeChartsReal) - uint32_t chartCount = 0, orthoChartsCount = 0, planarChartsCount = 0, lscmChartsCount = 0, piecewiseChartsCount = 0, chartsAddedCount = 0, chartsDeletedCount = 0; - for (uint32_t i = 0; i < ctx->meshCount; i++) { + const uint32_t meshCount = ctx->meshes.size(); + uint32_t chartCount = 0, chartsWithHolesCount = 0, holesCount = 0, chartsWithTJunctionsCount = 0, tJunctionsCount = 0, orthoChartsCount = 0, planarChartsCount = 0, lscmChartsCount = 0, piecewiseChartsCount = 0, chartsAddedCount = 0, chartsDeletedCount = 0; + for (uint32_t i = 0; i < meshCount; i++) { for (uint32_t j = 0; j < ctx->paramAtlas.chartGroupCount(i); j++) { const internal::param::ChartGroup *chartGroup = ctx->paramAtlas.chartGroupAt(i, j); - if (chartGroup->isVertexMap()) - continue; for (uint32_t k = 0; k < chartGroup->chartCount(); k++) { const internal::param::Chart *chart = chartGroup->chartAt(k); +#if XA_PRINT_CHART_WARNINGS + if (chart->warningFlags() & internal::param::ChartWarningFlags::CloseHolesFailed) + XA_PRINT_WARNING(" Chart %u (mesh %u, group %u, id %u): failed to close holes\n", chartCount, i, j, k); + if (chart->warningFlags() & internal::param::ChartWarningFlags::FixTJunctionsDuplicatedEdge) + XA_PRINT_WARNING(" Chart %u (mesh %u, group %u, id %u): fixing t-junctions created non-manifold geometry\n", chartCount, i, j, k); + if (chart->warningFlags() & internal::param::ChartWarningFlags::FixTJunctionsFailed) + XA_PRINT_WARNING(" Chart %u (mesh %u, group %u, id %u): fixing t-junctions failed\n", chartCount, i, j, k); + if (chart->warningFlags() & internal::param::ChartWarningFlags::TriangulateDuplicatedEdge) + XA_PRINT_WARNING(" Chart %u (mesh %u, group %u, id %u): triangulation created non-manifold geometry\n", chartCount, i, j, k); +#endif + holesCount += chart->closedHolesCount(); + if (chart->closedHolesCount() > 0) + chartsWithHolesCount++; + tJunctionsCount += chart->fixedTJunctionsCount(); + if (chart->fixedTJunctionsCount() > 0) + chartsWithTJunctionsCount++; if (chart->type() == ChartType::Planar) planarChartsCount++; else if (chart->type() == ChartType::Ortho) @@ -9157,19 +9766,21 @@ void ParameterizeCharts(Atlas *atlas, ParameterizeFunc func) chartsDeletedCount += chartGroup->paramDeletedChartsCount(); } } + if (holesCount > 0) + XA_PRINT(" %u holes closed in %u charts\n", holesCount, chartsWithHolesCount); + if (tJunctionsCount > 0) + XA_PRINT(" %u t-junctions fixed in %u charts\n", tJunctionsCount, chartsWithTJunctionsCount); XA_PRINT(" %u planar charts, %u ortho charts, %u LSCM charts, %u piecewise charts\n", planarChartsCount, orthoChartsCount, lscmChartsCount, piecewiseChartsCount); if (chartsDeletedCount > 0) { XA_PRINT(" %u charts with invalid parameterizations replaced with %u new charts\n", chartsDeletedCount, chartsAddedCount); XA_PRINT(" %u charts\n", chartCount); } uint32_t chartIndex = 0, invalidParamCount = 0; - for (uint32_t i = 0; i < ctx->meshCount; i++) { + for (uint32_t i = 0; i < meshCount; i++) { for (uint32_t j = 0; j < ctx->paramAtlas.chartGroupCount(i); j++) { const internal::param::ChartGroup *chartGroup = ctx->paramAtlas.chartGroupAt(i, j); - if (chartGroup->isVertexMap()) - continue; for (uint32_t k = 0; k < chartGroup->chartCount(); k++) { - const internal::param::Chart *chart = chartGroup->chartAt(k); + internal::param::Chart *chart = chartGroup->chartAt(k); const internal::param::Quality &quality = chart->quality(); #if XA_DEBUG_EXPORT_OBJ_CHARTS_AFTER_PARAMETERIZATION { @@ -9178,7 +9789,6 @@ void ParameterizeCharts(Atlas *atlas, ParameterizeFunc func) chart->unifiedMesh()->writeObjFile(filename); } #endif - bool invalid = false; const char *type = "LSCM"; if (chart->type() == ChartType::Planar) type = "planar"; @@ -9186,18 +9796,15 @@ void ParameterizeCharts(Atlas *atlas, ParameterizeFunc func) type = "ortho"; else if (chart->type() == ChartType::Piecewise) type = "piecewise"; - if (quality.boundaryIntersection) { - invalid = true; - XA_PRINT_WARNING(" Chart %u (mesh %u, group %u, id %u) (%s): invalid parameterization, self-intersecting boundary.\n", chartIndex, i, j, k, type); - } - if (quality.flippedTriangleCount > 0) { - invalid = true; - XA_PRINT_WARNING(" Chart %u (mesh %u, group %u, id %u) (%s): invalid parameterization, %u / %u flipped triangles.\n", chartIndex, i, j, k, type, quality.flippedTriangleCount, quality.totalTriangleCount); - } - if (invalid) + if (chart->isInvalid()) { + if (quality.boundaryIntersection) { + XA_PRINT_WARNING(" Chart %u (mesh %u, group %u, id %u) (%s): invalid parameterization, self-intersecting boundary.\n", chartIndex, i, j, k, type); + } + if (quality.flippedTriangleCount > 0) { + XA_PRINT_WARNING(" Chart %u (mesh %u, group %u, id %u) (%s): invalid parameterization, %u / %u flipped triangles.\n", chartIndex, i, j, k, type, quality.flippedTriangleCount, quality.totalTriangleCount); + } invalidParamCount++; #if XA_DEBUG_EXPORT_OBJ_INVALID_PARAMETERIZATION - if (invalid) { char filename[256]; XA_SPRINTF(filename, sizeof(filename), "debug_chart_%03u_invalid_parameterization.obj", chartIndex); const internal::Mesh *mesh = chart->unifiedMesh(); @@ -9218,8 +9825,8 @@ void ParameterizeCharts(Atlas *atlas, ParameterizeFunc func) mesh->writeObjLinkedBoundaries(file); fclose(file); } - } #endif + } chartIndex++; } } @@ -9228,9 +9835,18 @@ void ParameterizeCharts(Atlas *atlas, ParameterizeFunc func) XA_PRINT_WARNING(" %u charts with invalid parameterizations\n", invalidParamCount); XA_PROFILE_PRINT_AND_RESET(" Total (real): ", parameterizeChartsReal) XA_PROFILE_PRINT_AND_RESET(" Total (thread): ", parameterizeChartsThread) + XA_PROFILE_PRINT_AND_RESET(" Create chart mesh: ", createChartMesh) + XA_PROFILE_PRINT_AND_RESET(" Fix t-junctions: ", fixChartMeshTJunctions) + XA_PROFILE_PRINT_AND_RESET(" Close holes: ", closeChartMeshHoles) XA_PROFILE_PRINT_AND_RESET(" Orthogonal: ", parameterizeChartsOrthogonal) XA_PROFILE_PRINT_AND_RESET(" LSCM: ", parameterizeChartsLSCM) + XA_PROFILE_PRINT_AND_RESET(" Recompute: ", parameterizeChartsRecompute) + XA_PROFILE_PRINT_AND_RESET(" Piecewise: ", parameterizeChartsPiecewise) + XA_PROFILE_PRINT_AND_RESET(" Boundary intersection: ", parameterizeChartsPiecewiseBoundaryIntersection) XA_PROFILE_PRINT_AND_RESET(" Evaluate quality: ", parameterizeChartsEvaluateQuality) +#if XA_PROFILE_ALLOC + XA_PROFILE_PRINT_AND_RESET(" Alloc: ", alloc) +#endif XA_PRINT_MEM_USAGE } @@ -9242,7 +9858,7 @@ void PackCharts(Atlas *atlas, PackOptions packOptions) return; } Context *ctx = (Context *)atlas; - if (ctx->meshCount == 0 && ctx->uvMeshInstances.isEmpty()) { + if (ctx->meshes.isEmpty() && ctx->uvMeshInstances.isEmpty()) { XA_PRINT_WARNING("PackCharts: No meshes. Call AddMesh or AddUvMesh first.\n"); return; } @@ -9299,7 +9915,7 @@ void PackCharts(Atlas *atlas, PackOptions packOptions) if (packOptions.createImage) { atlas->image = XA_ALLOC_ARRAY(internal::MemTag::Default, uint32_t, atlas->atlasCount * atlas->width * atlas->height); for (uint32_t i = 0; i < atlas->atlasCount; i++) - packAtlas.getImages()[i]->copyTo(&atlas->image[atlas->width * atlas->height * i], atlas->width, atlas->height, packOptions.blockAlign ? 0 : packOptions.padding); + packAtlas.getImages()[i]->copyTo(&atlas->image[atlas->width * atlas->height * i], atlas->width, atlas->height, packOptions.padding); } XA_PROFILE_PRINT_AND_RESET(" Total: ", packCharts) XA_PROFILE_PRINT_AND_RESET(" Add charts (real): ", packChartsAddCharts) @@ -9309,6 +9925,9 @@ void PackCharts(Atlas *atlas, PackOptions packOptions) XA_PROFILE_PRINT_AND_RESET(" Dilate (padding): ", packChartsDilate) XA_PROFILE_PRINT_AND_RESET(" Find location: ", packChartsFindLocation) XA_PROFILE_PRINT_AND_RESET(" Blit: ", packChartsBlit) +#if XA_PROFILE_ALLOC + XA_PROFILE_PRINT_AND_RESET(" Alloc: ", alloc) +#endif XA_PRINT_MEM_USAGE XA_PRINT("Building output meshes\n"); XA_PROFILE_START(buildOutputMeshes) @@ -9318,91 +9937,92 @@ void PackCharts(Atlas *atlas, PackOptions packOptions) return; } if (ctx->uvMeshInstances.isEmpty()) - atlas->meshCount = ctx->meshCount; + atlas->meshCount = ctx->meshes.size(); else atlas->meshCount = ctx->uvMeshInstances.size(); atlas->meshes = XA_ALLOC_ARRAY(internal::MemTag::Default, Mesh, atlas->meshCount); memset(atlas->meshes, 0, sizeof(Mesh) * atlas->meshCount); if (ctx->uvMeshInstances.isEmpty()) { uint32_t chartIndex = 0; - for (uint32_t i = 0; i < ctx->meshCount; i++) { + for (uint32_t i = 0; i < atlas->meshCount; i++) { Mesh &outputMesh = atlas->meshes[i]; - // Count and alloc arrays. Ignore vertex mapped chart groups in Mesh::chartCount, since they're ignored faces. + // Count and alloc arrays. + const internal::param::InvalidMeshGeometry &invalid = ctx->paramAtlas.invalidMeshGeometry(i); + outputMesh.vertexCount += invalid.vertices().length; + outputMesh.indexCount += invalid.faces().length * 3; for (uint32_t cg = 0; cg < ctx->paramAtlas.chartGroupCount(i); cg++) { const internal::param::ChartGroup *chartGroup = ctx->paramAtlas.chartGroupAt(i, cg); - if (chartGroup->isVertexMap()) { - outputMesh.vertexCount += chartGroup->mesh()->vertexCount(); - outputMesh.indexCount += chartGroup->mesh()->faceCount() * 3; - } else { - for (uint32_t c = 0; c < chartGroup->chartCount(); c++) { - const internal::param::Chart *chart = chartGroup->chartAt(c); - outputMesh.vertexCount += chart->mesh()->vertexCount(); - outputMesh.indexCount += chart->mesh()->faceCount() * 3; - outputMesh.chartCount++; - } + for (uint32_t c = 0; c < chartGroup->chartCount(); c++) { + const internal::param::Chart *chart = chartGroup->chartAt(c); + outputMesh.vertexCount += chart->mesh()->vertexCount(); + outputMesh.indexCount += chart->mesh()->faceCount() * 3; + outputMesh.chartCount++; } } outputMesh.vertexArray = XA_ALLOC_ARRAY(internal::MemTag::Default, Vertex, outputMesh.vertexCount); outputMesh.indexArray = XA_ALLOC_ARRAY(internal::MemTag::Default, uint32_t, outputMesh.indexCount); outputMesh.chartArray = XA_ALLOC_ARRAY(internal::MemTag::Default, Chart, outputMesh.chartCount); - XA_PRINT(" mesh %u: %u vertices, %u triangles, %u charts\n", i, outputMesh.vertexCount, outputMesh.indexCount / 3, outputMesh.chartCount); + XA_PRINT(" Mesh %u: %u vertices, %u triangles, %u charts\n", i, outputMesh.vertexCount, outputMesh.indexCount / 3, outputMesh.chartCount); // Copy mesh data. - uint32_t firstVertex = 0, meshChartIndex = 0; + uint32_t firstVertex = 0; + { + const internal::param::InvalidMeshGeometry &mesh = ctx->paramAtlas.invalidMeshGeometry(i); + internal::ConstArrayView faces = mesh.faces(); + internal::ConstArrayView indices = mesh.indices(); + internal::ConstArrayView vertices = mesh.vertices(); + // Vertices. + for (uint32_t v = 0; v < vertices.length; v++) { + Vertex &vertex = outputMesh.vertexArray[v]; + vertex.atlasIndex = -1; + vertex.chartIndex = -1; + vertex.uv[0] = vertex.uv[1] = 0.0f; + vertex.xref = vertices[v]; + } + // Indices. + for (uint32_t f = 0; f < faces.length; f++) { + const uint32_t indexOffset = faces[f] * 3; + for (uint32_t j = 0; j < 3; j++) + outputMesh.indexArray[indexOffset + j] = indices[f * 3 + j]; + } + firstVertex = vertices.length; + } + uint32_t meshChartIndex = 0; for (uint32_t cg = 0; cg < ctx->paramAtlas.chartGroupCount(i); cg++) { const internal::param::ChartGroup *chartGroup = ctx->paramAtlas.chartGroupAt(i, cg); - if (chartGroup->isVertexMap()) { - const internal::Mesh *mesh = chartGroup->mesh(); + for (uint32_t c = 0; c < chartGroup->chartCount(); c++) { + const internal::param::Chart *chart = chartGroup->chartAt(c); + const internal::Mesh *mesh = chart->mesh(); // Vertices. for (uint32_t v = 0; v < mesh->vertexCount(); v++) { Vertex &vertex = outputMesh.vertexArray[firstVertex + v]; - vertex.atlasIndex = -1; - vertex.chartIndex = -1; - vertex.uv[0] = vertex.uv[1] = 0.0f; - vertex.xref = chartGroup->mapVertexToSourceVertex(v); + vertex.atlasIndex = packAtlas.getChart(chartIndex)->atlasIndex; + XA_DEBUG_ASSERT(vertex.atlasIndex >= 0); + vertex.chartIndex = (int32_t)chartIndex; + const internal::Vector2 &uv = mesh->texcoord(v); + vertex.uv[0] = internal::max(0.0f, uv.x); + vertex.uv[1] = internal::max(0.0f, uv.y); + vertex.xref = chart->mapChartVertexToSourceVertex(v); } // Indices. for (uint32_t f = 0; f < mesh->faceCount(); f++) { - const uint32_t indexOffset = chartGroup->mapFaceToSourceFace(f) * 3; + const uint32_t indexOffset = chart->mapFaceToSourceFace(f) * 3; for (uint32_t j = 0; j < 3; j++) outputMesh.indexArray[indexOffset + j] = firstVertex + mesh->vertexAt(f * 3 + j); } + // Charts. + Chart *outputChart = &outputMesh.chartArray[meshChartIndex]; + const int32_t atlasIndex = packAtlas.getChart(chartIndex)->atlasIndex; + XA_DEBUG_ASSERT(atlasIndex >= 0); + outputChart->atlasIndex = (uint32_t)atlasIndex; + outputChart->type = chart->isInvalid() ? ChartType::Invalid : chart->type(); + outputChart->faceCount = mesh->faceCount(); + outputChart->faceArray = XA_ALLOC_ARRAY(internal::MemTag::Default, uint32_t, outputChart->faceCount); + for (uint32_t f = 0; f < outputChart->faceCount; f++) + outputChart->faceArray[f] = chart->mapFaceToSourceFace(f); + outputChart->material = 0; + meshChartIndex++; + chartIndex++; firstVertex += mesh->vertexCount(); - } else { - for (uint32_t c = 0; c < chartGroup->chartCount(); c++) { - const internal::param::Chart *chart = chartGroup->chartAt(c); - const internal::Mesh *mesh = chart->mesh(); - // Vertices. - for (uint32_t v = 0; v < mesh->vertexCount(); v++) { - Vertex &vertex = outputMesh.vertexArray[firstVertex + v]; - vertex.atlasIndex = packAtlas.getChart(chartIndex)->atlasIndex; - XA_DEBUG_ASSERT(vertex.atlasIndex >= 0); - vertex.chartIndex = (int32_t)chartIndex; - const internal::Vector2 &uv = mesh->texcoord(v); - vertex.uv[0] = internal::max(0.0f, uv.x); - vertex.uv[1] = internal::max(0.0f, uv.y); - vertex.xref = chartGroup->mapVertexToSourceVertex(chart->mapChartVertexToOriginalVertex(v)); - } - // Indices. - for (uint32_t f = 0; f < mesh->faceCount(); f++) { - const uint32_t indexOffset = chartGroup->mapFaceToSourceFace(chart->mapFaceToSourceFace(f)) * 3; - for (uint32_t j = 0; j < 3; j++) - outputMesh.indexArray[indexOffset + j] = firstVertex + mesh->vertexAt(f * 3 + j); - } - // Charts. - Chart *outputChart = &outputMesh.chartArray[meshChartIndex]; - const int32_t atlasIndex = packAtlas.getChart(chartIndex)->atlasIndex; - XA_DEBUG_ASSERT(atlasIndex >= 0); - outputChart->atlasIndex = (uint32_t)atlasIndex; - outputChart->type = chart->type(); - outputChart->faceCount = mesh->faceCount(); - outputChart->faceArray = XA_ALLOC_ARRAY(internal::MemTag::Default, uint32_t, outputChart->faceCount); - for (uint32_t f = 0; f < outputChart->faceCount; f++) - outputChart->faceArray[f] = chartGroup->mapFaceToSourceFace(chart->mapFaceToSourceFace(f)); - outputChart->material = 0; - meshChartIndex++; - chartIndex++; - firstVertex += mesh->vertexCount(); - } } } XA_DEBUG_ASSERT(outputMesh.vertexCount == firstVertex); @@ -9476,10 +10096,13 @@ void PackCharts(Atlas *atlas, PackOptions packOptions) ctx->progressFunc(ProgressCategory::BuildOutputMeshes, 100, ctx->progressUserData); XA_PROFILE_END(buildOutputMeshes) XA_PROFILE_PRINT_AND_RESET(" Total: ", buildOutputMeshes) +#if XA_PROFILE_ALLOC + XA_PROFILE_PRINT_AND_RESET(" Alloc: ", alloc) +#endif XA_PRINT_MEM_USAGE } -void Generate(Atlas *atlas, ChartOptions chartOptions, ParameterizeFunc paramFunc, PackOptions packOptions) +void Generate(Atlas *atlas, ChartOptions chartOptions, ParameterizeOptions parameterizeOptions, PackOptions packOptions) { if (!atlas) { XA_PRINT_WARNING("Generate: atlas is null.\n"); @@ -9490,12 +10113,12 @@ void Generate(Atlas *atlas, ChartOptions chartOptions, ParameterizeFunc paramFun XA_PRINT_WARNING("Generate: This function should not be called with UV meshes.\n"); return; } - if (ctx->meshCount == 0) { + if (ctx->meshes.isEmpty()) { XA_PRINT_WARNING("Generate: No meshes. Call AddMesh first.\n"); return; } ComputeCharts(atlas, chartOptions); - ParameterizeCharts(atlas, paramFunc); + ParameterizeCharts(atlas, parameterizeOptions); PackCharts(atlas, packOptions); } diff --git a/thirdparty/xatlas/xatlas.h b/thirdparty/xatlas/xatlas.h index e59f493287..cc47f4837e 100644 --- a/thirdparty/xatlas/xatlas.h +++ b/thirdparty/xatlas/xatlas.h @@ -1,7 +1,7 @@ /* MIT License -Copyright (c) 2018-2019 Jonathan Young +Copyright (c) 2018-2020 Jonathan Young Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal @@ -42,18 +42,19 @@ struct ChartType Planar, Ortho, LSCM, - Piecewise + Piecewise, + Invalid }; }; // A group of connected faces, belonging to a single atlas. struct Chart { - uint32_t atlasIndex; // Sub-atlas index. uint32_t *faceArray; + uint32_t atlasIndex; // Sub-atlas index. uint32_t faceCount; - uint32_t material; ChartType::Enum type; + uint32_t material; }; // Output vertex. @@ -69,10 +70,10 @@ struct Vertex struct Mesh { Chart *chartArray; - uint32_t chartCount; uint32_t *indexArray; - uint32_t indexCount; Vertex *vertexArray; + uint32_t chartCount; + uint32_t indexCount; uint32_t vertexCount; }; @@ -84,15 +85,15 @@ static const uint32_t kImageIsPaddingBit = 0x20000000; // Empty on creation. Populated after charts are packed. struct Atlas { + uint32_t *image; + Mesh *meshes; // The output meshes, corresponding to each AddMesh call. uint32_t width; // Atlas width in texels. uint32_t height; // Atlas height in texels. uint32_t atlasCount; // Number of sub-atlases. Equal to 0 unless PackOptions resolution is changed from default (0). uint32_t chartCount; // Total number of charts in all meshes. uint32_t meshCount; // Number of output meshes. Equal to the number of times AddMesh was called. - Mesh *meshes; // The output meshes, corresponding to each AddMesh call. float *utilization; // Normalized atlas texel utilization array. E.g. a value of 0.8 means 20% empty space. atlasCount in length. float texelsPerUnit; // Equal to PackOptions texelsPerUnit if texelsPerUnit > 0, otherwise an estimated value to match PackOptions resolution. - uint32_t *image; }; // Create an empty atlas. @@ -112,22 +113,23 @@ struct IndexFormat // Input mesh declaration. struct MeshDecl { - uint32_t vertexCount = 0; const void *vertexPositionData = nullptr; - uint32_t vertexPositionStride = 0; const void *vertexNormalData = nullptr; // optional - uint32_t vertexNormalStride = 0; // optional const void *vertexUvData = nullptr; // optional. The input UVs are provided as a hint to the chart generator. - uint32_t vertexUvStride = 0; // optional - uint32_t indexCount = 0; const void *indexData = nullptr; // optional - int32_t indexOffset = 0; // optional. Add this offset to all indices. - IndexFormat::Enum indexFormat = IndexFormat::UInt16; // Optional. indexCount / 3 (triangle count) in length. // Don't atlas faces set to true. Ignored faces still exist in the output meshes, Vertex uv is set to (0, 0) and Vertex atlasIndex to -1. const bool *faceIgnoreData = nullptr; + uint32_t vertexCount = 0; + uint32_t vertexPositionStride = 0; + uint32_t vertexNormalStride = 0; // optional + uint32_t vertexUvStride = 0; // optional + uint32_t indexCount = 0; + int32_t indexOffset = 0; // optional. Add this offset to all indices. + IndexFormat::Enum indexFormat = IndexFormat::UInt16; + // Vertex positions within epsilon distance of each other are considered colocal. float epsilon = 1.192092896e-07F; }; @@ -151,14 +153,14 @@ void AddMeshJoin(Atlas *atlas); struct UvMeshDecl { + const void *vertexUvData = nullptr; + const void *indexData = nullptr; // optional + const uint32_t *faceMaterialData = nullptr; // Optional. Faces with different materials won't be assigned to the same chart. Must be indexCount / 3 in length. uint32_t vertexCount = 0; uint32_t vertexStride = 0; - const void *vertexUvData = nullptr; uint32_t indexCount = 0; - const void *indexData = nullptr; // optional int32_t indexOffset = 0; // optional. Add this offset to all indices. IndexFormat::Enum indexFormat = IndexFormat::UInt16; - const uint32_t *faceMaterialData = nullptr; // Optional. Faces with different materials won't be assigned to the same chart. Must be indexCount / 3 in length. bool rotateCharts = true; }; @@ -170,24 +172,31 @@ struct ChartOptions float maxBoundaryLength = 0.0f; // Don't grow charts to have a longer boundary than this. 0 means no limit. // Weights determine chart growth. Higher weights mean higher cost for that metric. - float proxyFitMetricWeight = 2.0f; // Angle between face and average chart normal. - float roundnessMetricWeight = 0.01f; - float straightnessMetricWeight = 6.0f; - float normalSeamMetricWeight = 4.0f; // If > 1000, normal seams are fully respected. - float textureSeamMetricWeight = 0.5f; + float normalDeviationWeight = 2.0f; // Angle between face and average chart normal. + float roundnessWeight = 0.01f; + float straightnessWeight = 6.0f; + float normalSeamWeight = 4.0f; // If > 1000, normal seams are fully respected. + float textureSeamWeight = 0.5f; - float maxThreshold = 2.0f; // If total of all metrics * weights > maxThreshold, don't grow chart. Lower values result in more charts. + float maxCost = 2.0f; // If total of all metrics * weights > maxCost, don't grow chart. Lower values result in more charts. uint32_t maxIterations = 1; // Number of iterations of the chart growing and seeding phases. Higher values result in better charts. }; // Call after all AddMesh calls. Can be called multiple times to recompute charts with different options. -void ComputeCharts(Atlas *atlas, ChartOptions chartOptions = ChartOptions()); +void ComputeCharts(Atlas *atlas, ChartOptions options = ChartOptions()); // Custom parameterization function. texcoords initial values are an orthogonal parameterization. typedef void (*ParameterizeFunc)(const float *positions, float *texcoords, uint32_t vertexCount, const uint32_t *indices, uint32_t indexCount); +struct ParameterizeOptions +{ + ParameterizeFunc func = nullptr; + bool closeHoles = true; // If the custom parameterization function works with multiple boundaries, this can be set to false to improve performance. + bool fixTJunctions = true; // If meshes don't have T-junctions, this can be set to false to improve performance. +}; + // Call after ComputeCharts. Can be called multiple times to re-parameterize charts with a different ParameterizeFunc. -void ParameterizeCharts(Atlas *atlas, ParameterizeFunc func = nullptr); +void ParameterizeCharts(Atlas *atlas, ParameterizeOptions options = ParameterizeOptions()); struct PackOptions { @@ -224,7 +233,7 @@ struct PackOptions void PackCharts(Atlas *atlas, PackOptions packOptions = PackOptions()); // Equivalent to calling ComputeCharts, ParameterizeCharts and PackCharts in sequence. Can be called multiple times to regenerate with different options. -void Generate(Atlas *atlas, ChartOptions chartOptions = ChartOptions(), ParameterizeFunc paramFunc = nullptr, PackOptions packOptions = PackOptions()); +void Generate(Atlas *atlas, ChartOptions chartOptions = ChartOptions(), ParameterizeOptions parameterizeOptions = ParameterizeOptions(), PackOptions packOptions = PackOptions()); // Progress tracking. struct ProgressCategory -- cgit v1.2.3