refactor: segment processing and mesh deduplication

.github/workflows/build.yml

@@ -89,7 +89,7 @@ jobs:
       - name: 🏗️ Compile (other)
         if: ${{ matrix.os != 'windows-latest' }}
         run: |
-          cmake -DVIENNALS_BUILD_TESTS=ON -B build
+          cmake -DVIENNALS_BUILD_TESTS=ON -DCMAKE_BUILD_TYPE=${{ matrix.config }} -B build
           cmake --build build --config ${{ matrix.config }}
 
       - name: 🧪 Test

CMakeLists.txt

@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.20 FATAL_ERROR)
 project(
   ViennaLS
   LANGUAGES CXX
-  VERSION 5.8.5)
+  VERSION 5.9.0)
 
 # --------------------------------------------------------------------------------------------------------
 # Library options
@@ -521,11 +521,16 @@ endif()
 # Install Target
 # --------------------------------------------------------------------------------------------------------
 
+set(VIENNALS_DEPENDENCIES "ViennaHRLE;ViennaCore")
+if(VIENNALS_USE_VTK)
+  list(APPEND VIENNALS_DEPENDENCIES "VTK")
+endif()
+
 packageProject(
   NAME ${PROJECT_NAME} NAMESPACE ViennaTools
   VERSION ${PROJECT_VERSION}
   BINARY_DIR ${PROJECT_BINARY_DIR}
   INCLUDE_DIR ${PROJECT_SOURCE_DIR}/include/viennals
   INCLUDE_DESTINATION include/viennals-${PROJECT_VERSION}
   COMPATIBILITY SameMajorVersion
-  DEPENDENCIES "VTK;ViennaHRLE;ViennaCore")
+  DEPENDENCIES ${VIENNALS_DEPENDENCIES})

README.md

@@ -90,9 +90,6 @@ This will install the necessary headers and CMake files to the specified path. I
 
 In order to install ViennaLS without VTK, run:
 ```bash
-git clone https://github.com/ViennaTools/ViennaLS.git
-cd ViennaLS
-
 cmake -B build -D CMAKE_INSTALL_PREFIX=/path/to/your/custom/install/ -D VIENNALS_USE_VTK=OFF
 cmake --install build
 ```
@@ -152,9 +149,6 @@ ViennaLS uses CTest to run its tests.
 In order to check whether ViennaLS runs without issues on your system, you can run:
 
 ```bash
-git clone https://github.com/ViennaTools/ViennaLS.git
-cd ViennaLS
-
 cmake -B build -DVIENNALS_BUILD_TESTS=ON
 cmake --build build
 ctest -E "Benchmark|Performance" --test-dir build
@@ -175,7 +169,9 @@ We recommend using [CPM.cmake](https://github.com/cpm-cmake/CPM.cmake) to consum
 
 * Installation with CPM
   ```cmake
-  CPMAddPackage("gh:viennatools/viennals@5.8.5")
+  CPMAddPackage("gh:viennatools/viennals@5.9.0")
+
+  target_link_libraries(${PROJECT_NAME} PUBLIC ViennaTools::ViennaLS)
   ```
 
 * With a local installation
@@ -216,10 +212,6 @@ cmake --build build --target format
 
 ## Authors
 
-Current contributors: Tobias Reiter, Roman Kostal, Lado Filipovic
-
-Founder and initial developer: Otmar Ertl
-
 Contact us via: viennatools@iue.tuwien.ac.at
 
 ViennaLS was developed under the aegis of the 'Institute for Microelectronics' at the 'TU Wien'.  

include/viennals/lsAdvect.hpp

@@ -98,13 +98,10 @@ template <class T, int D> class Advect {
 
     VectorType<T, D> finalAlphas{};
 
-#pragma omp parallel num_threads((levelSets.back())->getNumberOfSegments())
-    {
+#pragma omp parallel for
+    for (unsigned p = 0; p < levelSets.back()->getNumberOfSegments(); ++p) {
       VectorType<T, D> localAlphas{};
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+
       viennahrle::Index<D> startVector =
           (p == 0) ? grid.getMinGridPoint()
                    : topDomain.getSegmentation()[p - 1];
@@ -274,20 +271,16 @@ template <class T, int D> class Advect {
     }
 #endif
 
-#pragma omp parallel num_threads(newDomain.getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < newDomain.getNumberOfSegments(); ++p) {
       auto &domainSegment = newDomain.getDomainSegment(p);
 
       viennahrle::Index<D> startVector =
           (p == 0) ? grid.getMinGridPoint()
                    : newDomain.getSegmentation()[p - 1];
 
       viennahrle::Index<D> endVector =
-          (p != static_cast<int>(newDomain.getNumberOfSegments() - 1))
+          (p != newDomain.getNumberOfSegments() - 1)
               ? newDomain.getSegmentation()[p]
               : grid.incrementIndices(grid.getMaxGridPoint());
 
@@ -461,18 +454,15 @@ template <class T, int D> class Advect {
 
     storedRates.resize(topDomain.getNumberOfSegments());
 
-#pragma omp parallel num_threads(topDomain.getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < topDomain.getNumberOfSegments(); ++p) {
+
       viennahrle::Index<D> startVector =
           (p == 0) ? grid.getMinGridPoint()
                    : topDomain.getSegmentation()[p - 1];
 
       viennahrle::Index<D> endVector =
-          (p != static_cast<int>(topDomain.getNumberOfSegments() - 1))
+          (p != topDomain.getNumberOfSegments() - 1)
               ? topDomain.getSegmentation()[p]
               : grid.incrementIndices(grid.getMaxGridPoint());
 
@@ -707,12 +697,9 @@ template <class T, int D> class Advect {
 
     const bool checkDiss = checkDissipation;
 
-#pragma omp parallel num_threads(topDomain.getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < topDomain.getNumberOfSegments(); ++p) {
+
       auto itRS = storedRates[p].cbegin();
       auto &segment = topDomain.getDomainSegment(p);
       const unsigned maxId = segment.getNumberOfPoints();
@@ -738,8 +725,8 @@ template <class T, int D> class Advect {
         T velocity = gradient - dissipation;
         // check if dissipation is too high and would cause a change in
         // direction of the velocity
-        if (checkDiss && (gradient < 0 && velocity > 0) ||
-            (gradient > 0 && velocity < 0)) {
+        if (checkDiss && ((gradient < 0 && velocity > 0) ||
+                          (gradient > 0 && velocity < 0))) {
           velocity = 0;
         }
 
@@ -751,8 +738,8 @@ template <class T, int D> class Advect {
 
           // recalculate velocity and rate
           velocity = itRS->first.first - itRS->first.second;
-          if (checkDiss && (itRS->first.first < 0 && velocity > 0) ||
-              (itRS->first.first > 0 && velocity < 0)) {
+          if (checkDiss && ((itRS->first.first < 0 && velocity > 0) ||
+                            (itRS->first.first > 0 && velocity < 0))) {
             velocity = 0;
           }
           rate = time * velocity;

include/viennals/lsBooleanOperation.hpp

@@ -74,12 +74,8 @@ template <class T, int D> class BooleanOperation {
       newDataLS.resize(newDataSourceIds.size());
     }
 
-#pragma omp parallel num_threads(newDomain.getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < newDomain.getNumberOfSegments(); ++p) {
 
       auto &domainSegment = newDomain.getDomainSegment(p);
 
@@ -214,12 +210,8 @@ template <class T, int D> class BooleanOperation {
 
   void invert() {
     auto &hrleDomain = levelSetA->getDomain();
-#pragma omp parallel num_threads(hrleDomain.getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < hrleDomain.getNumberOfSegments(); ++p) {
       auto &domainSegment = hrleDomain.getDomainSegment(p);
 
       // change all defined values

include/viennals/lsCalculateCurvatures.hpp

@@ -89,12 +89,8 @@ template <class T, int D> class CalculateCurvatures {
         (type == CurvatureEnum::MEAN_AND_GAUSSIAN_CURVATURE);
 
     //! Calculate Curvatures
-#pragma omp parallel num_threads(levelSet->getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < levelSet->getNumberOfSegments(); ++p) {
 
       auto &meanCurvatures = meanCurvaturesVector[p];
       auto &gaussCurvatures = gaussCurvaturesVector[p];

include/viennals/lsCalculateNormalVectors.hpp

@@ -131,12 +131,8 @@ template <class T, int D> class CalculateNormalVectors {
     auto grid = levelSet->getGrid();
 
     // Calculate Normalvectors
-#pragma omp parallel num_threads(levelSet->getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < levelSet->getNumberOfSegments(); ++p) {
 
       auto &normalVectors = normalVectorsVector[p];
       normalVectors.reserve(pointsPerSegment);
@@ -209,12 +205,8 @@ template <class T, int D> class CalculateNormalVectors {
     // points.
     std::vector<Vec3D<T>> normalVectors(levelSet->getNumberOfPoints());
 
-#pragma omp parallel num_threads(levelSet->getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < levelSet->getNumberOfSegments(); ++p) {
 
       viennahrle::Index<D> startVector =
           (p == 0) ? grid.getMinGridPoint()

include/viennals/lsCalculateVisibilities.hpp

@@ -65,12 +65,8 @@ template <class T, int D> class CalculateVisibilities {
       }
       //****************************
 
-#pragma omp parallel num_threads(levelSet->getNumberOfSegments())
-      {
-        int p = 0;
-#ifdef _OPENMP
-        p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+      for (unsigned p = 0; p < domain.getNumberOfSegments(); ++p) {
 
         const viennahrle::Index<D> startVector =
             (p == 0) ? grid.getMinGridPoint() : domain.getSegmentation()[p - 1];

include/viennals/lsDetectFeatures.hpp

@@ -94,12 +94,8 @@ template <class T, int D> class DetectFeatures {
     typename Domain<T, D>::DomainType &domain = levelSet->getDomain();
     std::vector<std::vector<T>> flagsReserve(levelSet->getNumberOfSegments());
 
-#pragma omp parallel num_threads((levelSet)->getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < levelSet->getNumberOfSegments(); ++p) {
 
       auto &flagsSegment = flagsReserve[p];
       flagsSegment.reserve(
@@ -169,12 +165,8 @@ template <class T, int D> class DetectFeatures {
     std::vector<std::vector<T>> flagsReserve(levelSet->getNumberOfSegments());
 
     // Compare angles between normal vectors
-#pragma omp parallel num_threads((levelSet)->getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < levelSet->getNumberOfSegments(); ++p) {
 
       Vec3D<T> zeroVector{};
 

include/viennals/lsExpand.hpp

@@ -81,12 +81,8 @@ template <class T, int D> class Expand {
       if (updateData)
         newDataSourceIds.resize(newDomain.getNumberOfSegments());
 
-#pragma omp parallel num_threads(newDomain.getNumberOfSegments())
-      {
-        int p = 0;
-#ifdef _OPENMP
-        p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+      for (unsigned p = 0; p < newDomain.getNumberOfSegments(); ++p) {
 
         auto &domainSegment = newDomain.getDomainSegment(p);
 

include/viennals/lsGeometricAdvect.hpp

@@ -297,12 +297,8 @@ template <class T, int D> class GeometricAdvect {
     }
 #endif
 // set up multithreading
-#pragma omp parallel num_threads(domain.getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
+#pragma omp parallel for
+    for (unsigned p = 0; p < domain.getNumberOfSegments(); ++p) {
 
       viennahrle::Index<D> startVector;
       if (p == 0) {

include/viennals/lsInterior.hpp

@@ -73,14 +73,8 @@ template <class T, int D> class Interior {
       if (updateData)
         newDataSourceIds.resize(newDomain.getNumberOfSegments());
 
-#pragma omp parallel num_threads(newDomain.getNumberOfSegments())              \
-    reduction(+ : addedPoints)
-      {
-        int p = 0;
-#ifdef _OPENMP
-        p = omp_get_thread_num();
-#endif
-
+#pragma omp parallel for reduction(+ : addedPoints)
+      for (unsigned p = 0; p < newDomain.getNumberOfSegments(); ++p) {
         auto &domainSegment = newDomain.getDomainSegment(p);
 
         viennahrle::Index<D> const startVector =

include/viennals/lsMaterialMap.hpp

@@ -45,6 +45,9 @@ class MaterialMap {
   void setMaterialId(const std::size_t index, const int materialId) {
     if (index >= materialMap.size()) {
       materialMap.resize(index + 1, -1); // Initialize new elements with -1
+    } else {
+      auto oldMaterialId = materialMap[index];
+      materials.erase(oldMaterialId); // Remove old material ID if it exists
     }
     materialMap[index] = materialId;
     materials.insert(materialId);