Merge pull request #1254 from p12tic/fix-int-division-ceil

Fix wrong integer division where rounding up was intended
alicevision · Oct 6, 2022 · ec91fab · ec91fab
2 parents fbc7ec5 + 2caea25
commit ec91fab
Show file tree

Hide file tree

Showing 3 changed files with 81 additions and 5 deletions.
diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -958,8 +958,8 @@ void DelaunayGraphCut::addMaskHelperPoints(const Point3d voxel[8], const StaticV
                 }
             }
 
-            int syMax = std::ceil(height / step);
-            int sxMax = std::ceil(width / step);
+            int syMax = divideRoundUp(height, step);
+            int sxMax = divideRoundUp(width, step);
 
             for(int sy = 0; sy < syMax; ++sy)
             {
@@ -1063,7 +1063,8 @@ void DelaunayGraphCut::fuseFromDepthMaps(const StaticVector<int>& cams, const Po
     {
         const auto& imgParams = _mp.getImageParams(i);
         startIndex[i] = realMaxVertices;
-        realMaxVertices += std::ceil(imgParams.width / step) * std::ceil(imgParams.height / step);
+        realMaxVertices += divideRoundUp(imgParams.width, step) *
+                           divideRoundUp(imgParams.height, step);
     }
     std::vector<Point3d> verticesCoordsPrepare(realMaxVertices);
     std::vector<double> pixSizePrepare(realMaxVertices);
@@ -1135,8 +1136,8 @@ void DelaunayGraphCut::fuseFromDepthMaps(const StaticVector<int>& cams, const Po
                 }
             }
 
-            int syMax = std::ceil(height/step);
-            int sxMax = std::ceil(width/step);
+            int syMax = divideRoundUp(height, step);
+            int sxMax = divideRoundUp(width, step);
             #pragma omp parallel for
             for(int sy = 0; sy < syMax; ++sy)
             {

diff --git a/src/aliceVision/numeric/numeric.hpp b/src/aliceVision/numeric/numeric.hpp
@@ -570,6 +570,20 @@ void SplitRange(const T range_start, const T range_end, const int nb_split,
   }
 }
 
+template<class T>
+constexpr T divideRoundUp(T x, T y)
+{
+    static_assert(std::is_integral<T>::value, "divideRoundUp only works with integer arguments");
+    const auto xPos = x >= 0;
+    const auto yPos = y >= 0;
+    if (xPos == yPos) {
+        return x / y + T((x % y) != 0);
+    } else {
+        // negative result, rounds towards zero anyways
+        return x / y;
+    }
+}
+
 /**
  * This function initializes the global state of random number generators that e.g. our tests
  * depend on. This makes it possible to have exactly reproducible program runtime behavior

diff --git a/src/aliceVision/numeric/numeric_test.cpp b/src/aliceVision/numeric/numeric_test.cpp
@@ -113,3 +113,64 @@ BOOST_AUTO_TEST_CASE(Numeric_MeanAndVarianceAlongRows) {
   BOOST_CHECK_SMALL(0.25-variance(0), 1e-8);
   BOOST_CHECK_SMALL(1.25-variance(1), 1e-8);
 }
+
+BOOST_AUTO_TEST_CASE(Numeric_divideRoundUp)
+{
+    BOOST_CHECK_EQUAL(divideRoundUp(0, 1), 0);
+    BOOST_CHECK_EQUAL(divideRoundUp(1, 1), 1);
+    BOOST_CHECK_EQUAL(divideRoundUp(2, 1), 2);
+    BOOST_CHECK_EQUAL(divideRoundUp(0, 2), 0);
+    BOOST_CHECK_EQUAL(divideRoundUp(1, 2), 1);
+    BOOST_CHECK_EQUAL(divideRoundUp(2, 2), 1);
+    BOOST_CHECK_EQUAL(divideRoundUp(3, 2), 2);
+    BOOST_CHECK_EQUAL(divideRoundUp(4, 2), 2);
+    BOOST_CHECK_EQUAL(divideRoundUp(999, 1000), 1);
+    BOOST_CHECK_EQUAL(divideRoundUp(1000, 1000), 1);
+    BOOST_CHECK_EQUAL(divideRoundUp(1001, 1000), 2);
+    BOOST_CHECK_EQUAL(divideRoundUp(1000999, 1000), 1001);
+    BOOST_CHECK_EQUAL(divideRoundUp(1001000, 1000), 1001);
+    BOOST_CHECK_EQUAL(divideRoundUp(1001001, 1000), 1002);
+
+    BOOST_CHECK_EQUAL(divideRoundUp(-1, 1), -1);
+    BOOST_CHECK_EQUAL(divideRoundUp(-2, 1), -2);
+    BOOST_CHECK_EQUAL(divideRoundUp(-0, 2), 0);
+    BOOST_CHECK_EQUAL(divideRoundUp(-1, 2), 0);
+    BOOST_CHECK_EQUAL(divideRoundUp(-2, 2), -1);
+    BOOST_CHECK_EQUAL(divideRoundUp(-3, 2), -1);
+    BOOST_CHECK_EQUAL(divideRoundUp(-4, 2), -2);
+    BOOST_CHECK_EQUAL(divideRoundUp(-999, 1000), 0);
+    BOOST_CHECK_EQUAL(divideRoundUp(-1000, 1000), -1);
+    BOOST_CHECK_EQUAL(divideRoundUp(-1001, 1000), -1);
+    BOOST_CHECK_EQUAL(divideRoundUp(-1000999, 1000), -1000);
+    BOOST_CHECK_EQUAL(divideRoundUp(-1001000, 1000), -1001);
+    BOOST_CHECK_EQUAL(divideRoundUp(-1001001, 1000), -1001);
+
+    BOOST_CHECK_EQUAL(divideRoundUp(0, -1), 0);
+    BOOST_CHECK_EQUAL(divideRoundUp(1, -1), -1);
+    BOOST_CHECK_EQUAL(divideRoundUp(2, -1), -2);
+    BOOST_CHECK_EQUAL(divideRoundUp(0, -2), 0);
+    BOOST_CHECK_EQUAL(divideRoundUp(1, -2), 0);
+    BOOST_CHECK_EQUAL(divideRoundUp(2, -2), -1);
+    BOOST_CHECK_EQUAL(divideRoundUp(3, -2), -1);
+    BOOST_CHECK_EQUAL(divideRoundUp(4, -2), -2);
+    BOOST_CHECK_EQUAL(divideRoundUp(999, -1000), 0);
+    BOOST_CHECK_EQUAL(divideRoundUp(1000, -1000), -1);
+    BOOST_CHECK_EQUAL(divideRoundUp(1001, -1000), -1);
+    BOOST_CHECK_EQUAL(divideRoundUp(1000999, -1000), -1000);
+    BOOST_CHECK_EQUAL(divideRoundUp(1001000, -1000), -1001);
+    BOOST_CHECK_EQUAL(divideRoundUp(1001001, -1000), -1001);
+
+    BOOST_CHECK_EQUAL(divideRoundUp(-1, -1), 1);
+    BOOST_CHECK_EQUAL(divideRoundUp(-2, -1), 2);
+    BOOST_CHECK_EQUAL(divideRoundUp(0, -2), 0);
+    BOOST_CHECK_EQUAL(divideRoundUp(-1, -2), 1);
+    BOOST_CHECK_EQUAL(divideRoundUp(-2, -2), 1);
+    BOOST_CHECK_EQUAL(divideRoundUp(-3, -2), 2);
+    BOOST_CHECK_EQUAL(divideRoundUp(-4, -2), 2);
+    BOOST_CHECK_EQUAL(divideRoundUp(-999, -1000), 1);
+    BOOST_CHECK_EQUAL(divideRoundUp(-1000, -1000), 1);
+    BOOST_CHECK_EQUAL(divideRoundUp(-1001, -1000), 2);
+    BOOST_CHECK_EQUAL(divideRoundUp(-1000999, -1000), 1001);
+    BOOST_CHECK_EQUAL(divideRoundUp(-1001000, -1000), 1001);
+    BOOST_CHECK_EQUAL(divideRoundUp(-1001001, -1000), 1002);
+}