Merge pull request #1234 from p12tic/fix-openmp-atomics

Fix broken openmp atomic usage

INSTALL.md

@@ -28,7 +28,7 @@ Dependencies
 AliceVision depends on external libraries:
 
 * [Assimp >= 5.0.0](https://github.com/assimp/assimp)
-* [Boost >= 1.72.0](https://www.boost.org)
+* [Boost >= 1.73.0](https://www.boost.org)
 * [Ceres >= 1.10.0](https://github.com/ceres-solver/ceres-solver)
 * [Eigen >= 3.3.4](https://gitlab.com/libeigen/eigen)
 * [Geogram >= 1.7.5](https://gforge.inria.fr/frs/?group_id=5833)

src/aliceVision/alicevision_omp.hpp

@@ -25,15 +25,3 @@ inline void omp_destroy_lock(omp_lock_t *lock) {}
 inline void omp_set_lock(omp_lock_t *lock) {}
 inline void omp_unset_lock(omp_lock_t *lock) {}
 #endif
-
-// OpenMP >= 3.1 for advanced atomic clauses (https://software.intel.com/en-us/node/608160)
-// OpenMP preprocessor version: https://github.com/jeffhammond/HPCInfo/wiki/Preprocessor-Macros
-#if defined _OPENMP && _OPENMP >= 201107
-#define OMP_ATOMIC_UPDATE _Pragma("omp atomic update")
-#define OMP_ATOMIC_WRITE  _Pragma("omp atomic write")
-#define OMP_HAVE_MIN_MAX_REDUCTION
-#else
-#define OMP_ATOMIC_UPDATE _Pragma("omp atomic")
-#define OMP_ATOMIC_WRITE  _Pragma("omp atomic")
-#endif
-

src/aliceVision/fuseCut/DelaunayGraphCut.cpp

@@ -36,6 +36,7 @@
 #include <boost/math/constants/constants.hpp>
 #include <boost/accumulators/accumulators.hpp>
 #include <boost/accumulators/statistics.hpp>
+#include <boost/atomic/atomic_ref.hpp>
 
 namespace aliceVision {
 namespace fuseCut {
@@ -683,11 +684,8 @@ StaticVector<int> DelaunayGraphCut::getIsUsedPerCamera() const
         for(int c = 0; c < v.cams.size(); ++c)
         {
             const int obsCam = v.cams[c];
-
-//#pragma OMP_ATOMIC_WRITE
-            {
-                cams[obsCam] = 1;
-            }
+            // boost::atomic_ref<int>(cams[obsCam]) = 1;
+            cams[obsCam] = 1;
         }
     }
 
@@ -1892,6 +1890,9 @@ void DelaunayGraphCut::fillGraph(double nPixelSizeBehind, bool labatutWeights, b
         const int vertexIndex = verticesRandIds[i];
         const GC_vertexInfo& v = _verticesAttr[vertexIndex];
 
+        GeometriesCount subTotalGeometriesIntersectedFrontCount;
+        GeometriesCount subTotalGeometriesIntersectedBehindCount;
+
         if(v.isReal())
         {
             ++totalIsRealNrc;
@@ -1917,22 +1918,25 @@ void DelaunayGraphCut::fillGraph(double nPixelSizeBehind, bool labatutWeights, b
                 totalStepsBehind += stepsBehind;
                 totalRayBehind += 1;
 
-#pragma OMP_ATOMIC_UPDATE
-                totalGeometriesIntersectedFrontCount.facets += geometriesIntersectedFrontCount.facets;
-#pragma OMP_ATOMIC_UPDATE
-                totalGeometriesIntersectedFrontCount.vertices += geometriesIntersectedFrontCount.vertices;
-#pragma OMP_ATOMIC_UPDATE
-                totalGeometriesIntersectedFrontCount.edges += geometriesIntersectedFrontCount.edges;
-#pragma OMP_ATOMIC_UPDATE
-                totalGeometriesIntersectedBehindCount.facets += geometriesIntersectedBehindCount.facets;
-#pragma OMP_ATOMIC_UPDATE
-                totalGeometriesIntersectedBehindCount.vertices += geometriesIntersectedBehindCount.vertices;
-#pragma OMP_ATOMIC_UPDATE
-                totalGeometriesIntersectedBehindCount.edges += geometriesIntersectedBehindCount.edges;
+                subTotalGeometriesIntersectedFrontCount += geometriesIntersectedFrontCount;
+                subTotalGeometriesIntersectedBehindCount += geometriesIntersectedBehindCount;
             } // for c
 
             totalCamHaveVisibilityOnVertex += v.cams.size();
             totalOfVertex += 1;
+
+            boost::atomic_ref<std::size_t>{totalGeometriesIntersectedFrontCount.facets} +=
+                    subTotalGeometriesIntersectedFrontCount.facets;
+            boost::atomic_ref<std::size_t>{totalGeometriesIntersectedFrontCount.vertices} +=
+                    subTotalGeometriesIntersectedFrontCount.vertices;
+            boost::atomic_ref<std::size_t>{totalGeometriesIntersectedFrontCount.edges} +=
+                    subTotalGeometriesIntersectedFrontCount.edges;
+            boost::atomic_ref<std::size_t>{totalGeometriesIntersectedBehindCount.facets} +=
+                    subTotalGeometriesIntersectedBehindCount.facets;
+            boost::atomic_ref<std::size_t>{totalGeometriesIntersectedBehindCount.vertices} +=
+                    subTotalGeometriesIntersectedBehindCount.vertices;
+            boost::atomic_ref<std::size_t>{totalGeometriesIntersectedBehindCount.edges} +=
+                    subTotalGeometriesIntersectedBehindCount.edges;
         }
     }
 
@@ -2032,15 +2036,11 @@ void DelaunayGraphCut::fillGraphPartPtRc(
             if (geometry.type == EGeometryType::Facet)
             {
                 ++outFrontCount.facets;
-                {
-#pragma OMP_ATOMIC_UPDATE
-                    _cellsAttr[geometry.facet.cellIndex].emptinessScore += weight;
-                }
+                boost::atomic_ref<float>{_cellsAttr[geometry.facet.cellIndex].emptinessScore} += weight;
 
                 {
                     const float dist = distFcn(maxDist, (originPt - lastIntersectPt).size(), distFcnHeight);
-#pragma OMP_ATOMIC_UPDATE
-                    _cellsAttr[geometry.facet.cellIndex].gEdgeVisWeight[geometry.facet.localVertexIndex] += weight * dist;
+                    boost::atomic_ref<float>{_cellsAttr[geometry.facet.cellIndex].gEdgeVisWeight[geometry.facet.localVertexIndex]} += weight * dist;
                 }
 
                 // Take the mirror facet to iterate over the next cell
@@ -2067,8 +2067,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(
                 // These geometries do not have a cellIndex, so we use the previousGeometry to retrieve the cell between the previous geometry and the current one.
                 if (previousGeometry.type == EGeometryType::Facet)
                 {
-#pragma OMP_ATOMIC_UPDATE
-                    _cellsAttr[previousGeometry.facet.cellIndex].emptinessScore += weight;
+                    boost::atomic_ref<float>{_cellsAttr[previousGeometry.facet.cellIndex].emptinessScore} += weight;
                 }
 
                 if (geometry.type == EGeometryType::Vertex)
@@ -2096,8 +2095,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(
             if (lastIntersectedFacet.cellIndex != GEO::NO_CELL &&
                 (_mp.CArr[cam] - intersectPt).size() < 0.2 * pointCamDistance)
             {
-    #pragma OMP_ATOMIC_WRITE
-                _cellsAttr[lastIntersectedFacet.cellIndex].cellSWeight = (float)maxint;
+                boost::atomic_ref<float>{_cellsAttr[lastIntersectedFacet.cellIndex].cellSWeight} = (float)maxint;
             }
         }
 
@@ -2109,8 +2107,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(
                 // lastGeoIsVertex is supposed to be positive in almost all cases.
                 // If we do not reach the camera, we still vote on the last tetrehedra.
                 // Possible reaisons: the camera is not part of the vertices or we encounter a numerical error in intersectNextGeom
-#pragma OMP_ATOMIC_WRITE
-                _cellsAttr[lastIntersectedFacet.cellIndex].cellSWeight = (float)maxint;
+                boost::atomic_ref<float>{_cellsAttr[lastIntersectedFacet.cellIndex].cellSWeight} = (float)maxint;
             }
             // else
             // {
@@ -2183,15 +2180,11 @@ void DelaunayGraphCut::fillGraphPartPtRc(
                 // Vote for the first cell found (only once)
                 if (firstIteration)
                 {
-#pragma OMP_ATOMIC_UPDATE
-                    _cellsAttr[geometry.facet.cellIndex].on += fWeight;
+                    boost::atomic_ref<float>{_cellsAttr[geometry.facet.cellIndex].on} += fWeight;
                     firstIteration = false;
                 }
 
-                {
-#pragma OMP_ATOMIC_UPDATE
-                    _cellsAttr[geometry.facet.cellIndex].fullnessScore += fWeight;
-                }
+                boost::atomic_ref<float>{_cellsAttr[geometry.facet.cellIndex].fullnessScore} += fWeight;
 
                 // Take the mirror facet to iterate over the next cell
                 const Facet mFacet = mirrorFacet(geometry.facet);
@@ -2205,8 +2198,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(
 
                 {
                     const float dist = distFcn(maxDist, (originPt - lastIntersectPt).size(), distFcnHeight);
-#pragma OMP_ATOMIC_UPDATE
-                    _cellsAttr[geometry.facet.cellIndex].gEdgeVisWeight[geometry.facet.localVertexIndex] +=
+                    boost::atomic_ref<float>{_cellsAttr[geometry.facet.cellIndex].gEdgeVisWeight[geometry.facet.localVertexIndex]} +=
                         fWeight * dist;
                 }
                 if(previousGeometry.type == EGeometryType::Facet && outBehindCount.facets > 1000)
@@ -2235,8 +2227,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(
 
                     for (const CellIndex& ci : neighboringCells)
                     {
-#pragma OMP_ATOMIC_UPDATE
-                        _cellsAttr[neighboringCells[0]].on += fWeight;
+                        boost::atomic_ref<float>{_cellsAttr[neighboringCells[0]].on} += fWeight;
                     }
                     firstIteration = false;
                 }
@@ -2245,8 +2236,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(
                 // These geometries do not have a cellIndex, so we use the previousGeometry to retrieve the cell between the previous geometry and the current one.
                 if (previousGeometry.type == EGeometryType::Facet)
                 {
-#pragma OMP_ATOMIC_UPDATE
-                    _cellsAttr[previousGeometry.facet.cellIndex].fullnessScore += fWeight;
+                    boost::atomic_ref<float>{_cellsAttr[previousGeometry.facet.cellIndex].fullnessScore} += fWeight;
                 }
 
                 if (geometry.type == EGeometryType::Vertex)
@@ -2274,8 +2264,7 @@ void DelaunayGraphCut::fillGraphPartPtRc(
         // Vote for the last intersected facet (farthest from the camera)
         if (lastIntersectedFacet.cellIndex != GEO::NO_CELL)
         {
-#pragma OMP_ATOMIC_UPDATE
-            _cellsAttr[lastIntersectedFacet.cellIndex].cellTWeight += fWeight;
+            boost::atomic_ref<float>{_cellsAttr[lastIntersectedFacet.cellIndex].cellTWeight} += fWeight;
         }
     }
 }
@@ -2460,12 +2449,12 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(float nPixelSizeBehind)
                     }
                 }
                 ++totalRayFront;
-#pragma OMP_ATOMIC_UPDATE
-                totalGeometriesIntersectedFrontCount.facets += geometriesIntersectedFrontCount.facets;
-#pragma OMP_ATOMIC_UPDATE
-                totalGeometriesIntersectedFrontCount.vertices += geometriesIntersectedFrontCount.vertices;
-#pragma OMP_ATOMIC_UPDATE
-                totalGeometriesIntersectedFrontCount.edges += geometriesIntersectedFrontCount.edges;
+                boost::atomic_ref<std::size_t>{totalGeometriesIntersectedFrontCount.facets} +=
+                        geometriesIntersectedFrontCount.facets;
+                boost::atomic_ref<std::size_t>{totalGeometriesIntersectedFrontCount.vertices} +=
+                        geometriesIntersectedFrontCount.vertices;
+                boost::atomic_ref<std::size_t>{totalGeometriesIntersectedFrontCount.edges} +=
+                        geometriesIntersectedFrontCount.edges;
             }
             {
                 // Initialisation
@@ -2613,17 +2602,13 @@ void DelaunayGraphCut::forceTedgesByGradientIJCV(float nPixelSizeBehind)
                        (maxSilent < maxSilentPartRange)) // g < k_outl                  //// k_outl=100  // 400 in the paper
                         //(maxSilent-minSilent<maxSilentPartRange))
                     {
-#pragma OMP_ATOMIC_UPDATE
-                        _cellsAttr[lastIntersectedFacet.cellIndex].on += (maxJump - midSilent);
+                        boost::atomic_ref<float>{_cellsAttr[lastIntersectedFacet.cellIndex].on} += (maxJump - midSilent);
                     }
                 }
                 ++totalRayBehind;
-#pragma OMP_ATOMIC_UPDATE
-                totalGeometriesIntersectedBehindCount.facets += totalGeometriesIntersectedBehindCount.facets;
-#pragma OMP_ATOMIC_UPDATE
-                totalGeometriesIntersectedBehindCount.vertices += totalGeometriesIntersectedBehindCount.vertices;
-#pragma OMP_ATOMIC_UPDATE
-                totalGeometriesIntersectedBehindCount.edges += totalGeometriesIntersectedBehindCount.edges;
+                boost::atomic_ref<std::size_t>{totalGeometriesIntersectedBehindCount.facets} += totalGeometriesIntersectedBehindCount.facets;
+                boost::atomic_ref<std::size_t>{totalGeometriesIntersectedBehindCount.vertices} += totalGeometriesIntersectedBehindCount.vertices;
+                boost::atomic_ref<std::size_t>{totalGeometriesIntersectedBehindCount.edges} += totalGeometriesIntersectedBehindCount.edges;
             }
         }
         totalCamHaveVisibilityOnVertex += v.cams.size();
@@ -2928,7 +2913,7 @@ void DelaunayGraphCut::cellsStatusFilteringBySolidAngleRatio(int nbSolidAngleFil
     // using solid angle ratio between full/empty parts.
     for(int i = 0; i < nbSolidAngleFilteringIterations; ++i)
     {
-        std::vector<bool> cellsInvertStatus(_cellIsFull.size(), false);
+        std::vector<std::uint8_t> cellsInvertStatus(_cellIsFull.size(), false);
         int toInvertCount = 0;
 
         std::vector<bool> vertexIsOnSurface;
@@ -3025,8 +3010,7 @@ void DelaunayGraphCut::cellsStatusFilteringBySolidAngleRatio(int nbSolidAngleFil
             {
                 if(_cellIsFull[ci] == invertFull)
                 {
-#pragma OMP_ATOMIC_WRITE
-                    cellsInvertStatus[ci] = true;
+                    boost::atomic_ref<std::uint8_t>{cellsInvertStatus[ci]} = true;
                     ++toInvertCount;
                 }
             }

src/aliceVision/mesh/Mesh.cpp

@@ -13,6 +13,7 @@
 
 #include <geogram/points/kd_tree.h>
 
+#include <boost/atomic/atomic_ref.hpp>
 #include <boost/filesystem.hpp>
 #include <boost/algorithm/string/case_conv.hpp> 
 
@@ -2651,14 +2652,12 @@ bool Mesh::lockSurfaceBoundaries(int neighbourIterations, StaticVectorBool& out_
 
             if(boundariesVertices[edge.first])
             {
-                #pragma OMP_ATOMIC_WRITE
-                boundariesVerticesCurrent[edge.second] = true;
+                boost::atomic_ref<char>{boundariesVerticesCurrent[edge.second]} = true;
             }
 
             if(boundariesVertices[edge.second])
             {
-                #pragma OMP_ATOMIC_WRITE
-                boundariesVerticesCurrent[edge.first] = true;
+                boost::atomic_ref<char>{boundariesVerticesCurrent[edge.first]} = true;
             }
         }
         std::swap(boundariesVertices, boundariesVerticesCurrent);
@@ -2778,8 +2777,7 @@ bool Mesh::getSurfaceBoundaries(StaticVectorBool& out_trisToConsider, bool inver
 
         if(boundariesEdges[i] == !invert)
         {
-            #pragma OMP_ATOMIC_WRITE
-            out_trisToConsider[edge.triId] = true;
+            boost::atomic_ref<char>{out_trisToConsider[edge.triId]} = true;
         }
     }
 

src/software/pipeline/main_cameraInit.cpp

@@ -14,6 +14,7 @@
 #include <aliceVision/system/cmdline.hpp>
 #include <aliceVision/image/io.cpp>
 
+#include <boost/atomic/atomic_ref.hpp>
 #include <boost/program_options.hpp>
 #include <boost/filesystem.hpp>
 #include <boost/algorithm/string.hpp>
@@ -481,8 +482,7 @@ int aliceVision_main(int argc, char **argv)
         if(intrinsic->getFocalLengthPixX() > 0)
         {
           // the view intrinsic is initialized
-          #pragma omp atomic
-          ++completeViewCount;
+          boost::atomic_ref<std::size_t>(completeViewCount)++;
 
           // don't need to build a new intrinsic
           continue;
@@ -616,8 +616,7 @@ int aliceVision_main(int argc, char **argv)
     if(intrinsic && intrinsic->isValid())
     {
       // the view intrinsic is initialized
-      #pragma omp atomic
-      ++completeViewCount;
+      boost::atomic_ref<std::size_t>(completeViewCount)++;
     }
 
     // Create serial number if not already filled