maths: be explicit in the computation type of Vec3::getNorm2 + friends

The result and computation types of Vec3<T>::getNorm2 were assumed to
be T.  This works fine for floating point types, but leads to overflow
when the underlying vector is Vec3<int32_t>.  This commit requires the
caller to supply the result of the computation type as a template
parameter to getNorm2.

tmc3/PCCMath.h

@@ -79,12 +79,17 @@ public:
   const T& y() const { return data[1]; }
   const T& z() const { return data[2]; }
 
-  T getNorm() const { return static_cast<T>(sqrt(getNorm2())); }
-  T getNorm2() const { return (*this) * (*this); }
+  template<typename ResultT>
+  ResultT getNorm2() const
+  {
+    return Vec3<ResultT>(*this) * Vec3<ResultT>(*this);
+  }
+
   T getNormInf() const
   {
     return std::max(data[2], std::max(abs(data[0]), abs(data[1])));
   }
+
   Vec3& operator=(const Vec3& rhs)
   {
     memcpy(data, rhs.data, sizeof(data));
@@ -394,11 +399,13 @@ struct Box3 {
       && max.z() >= box.min.z() && min.z() <= box.max.z();
   }
 
-  T getDist2(const Vec3<T>& point) const
+  template<typename SquaredT>
+  SquaredT getDist2(const Vec3<T>& point) const
   {
-    const T dx = std::max(std::max(min[0] - point[0], T()), point[0] - max[0]);
-    const T dy = std::max(std::max(min[1] - point[1], T()), point[1] - max[1]);
-    const T dz = std::max(std::max(min[2] - point[2], T()), point[2] - max[2]);
+    using U = SquaredT;
+    U dx = U(std::max(std::max(min[0] - point[0], T()), point[0] - max[0]));
+    U dy = U(std::max(std::max(min[1] - point[1], T()), point[1] - max[1]));
+    U dz = U(std::max(std::max(min[2] - point[2], T()), point[2] - max[2]));
     return dx * dx + dy * dy + dz * dz;
   }
 

tmc3/PCCTMC3Common.h

@@ -490,7 +490,7 @@ FindNeighborWithinDistance(
     const int32_t index1 = retained[j];
     const int32_t pointIndex1 = packedVoxel[index1].index;
     const auto& point1 = pointCloud[pointIndex1];
-    const auto d2 = (point1 - point).getNorm2();
+    const auto d2 = (point1 - point).getNorm2<double>();
     if (d2 <= radius2) {
       return index1;
     }
@@ -666,7 +666,9 @@ computeNearestNeighbors(
         aps.scalable_lifting_enabled_flag, nodeSizeLog2,
         pointCloud[pointIndex1]);
 
-      double norm2 = times(point - point1, aps.lod_neigh_bias).getNorm2();
+      double norm2 =
+        times(point - point1, aps.lod_neigh_bias).getNorm2<double>();
+
       if (nodeSizeLog2 > 0 && point == point1) {
         norm2 = double(1 << (nodeSizeLog2 - 1));
         norm2 = norm2 * norm2;
@@ -686,7 +688,7 @@ computeNearestNeighbors(
         }
         if (
           predictor.neighborCount < aps.num_pred_nearest_neighbours
-          || bBoxes[bucketIndex1].getDist2(point)
+          || bBoxes[bucketIndex1].getDist2<int64_t>(point)
             <= predictor.neighbors[index0].weight) {
           const int32_t k0 = std::max(bucketIndex1 * bucketSize, j0);
           const int32_t k1 = std::min((bucketIndex1 + 1) * bucketSize, j1);
@@ -697,7 +699,8 @@ computeNearestNeighbors(
               pointCloud[pointIndex1]);
 
             double norm2 =
-              times(point - point1, aps.lod_neigh_bias).getNorm2();
+              times(point - point1, aps.lod_neigh_bias).getNorm2<double>();
+
             if (nodeSizeLog2 > 0 && point == point1) {
               norm2 = (double)(1 << (nodeSizeLog2 - 1));
               norm2 = norm2 * norm2;
@@ -720,9 +723,9 @@ computeNearestNeighbors(
       for (int32_t k = k0; k < k1; ++k) {
         const int32_t pointIndex1 = packedVoxel[indexes[startIndex + k]].index;
         const auto& point1 = pointCloud[pointIndex1];
+        auto d2 = times(point - point1, aps.lod_neigh_bias).getNorm2<double>();
         predictor.insertNeighbor(
-          pointIndex1, times(point - point1, aps.lod_neigh_bias).getNorm2(),
-          aps.num_pred_nearest_neighbours,
+          pointIndex1, d2, aps.num_pred_nearest_neighbours,
           startIndex + k - i + 2 * aps.search_range);
       }
 
@@ -734,7 +737,7 @@ computeNearestNeighbors(
 
         if (
           predictor.neighborCount < aps.num_pred_nearest_neighbours
-          || bBoxesI[bucketIndex1].getDist2(point)
+          || bBoxesI[bucketIndex1].getDist2<int64_t>(point)
             <= predictor.neighbors[index0].weight) {
           const int32_t k0 = bucketIndex1 * bucketSize;
           const int32_t k1 = std::min((bucketIndex1 + 1) * bucketSize, j1);
@@ -742,10 +745,11 @@ computeNearestNeighbors(
             const int32_t pointIndex1 =
               packedVoxel[indexes[startIndex + k]].index;
             const auto& point1 = pointCloud[pointIndex1];
+            auto d2 =
+              times(point - point1, aps.lod_neigh_bias).getNorm2<double>();
+
             predictor.insertNeighbor(
-              pointIndex1,
-              times(point - point1, aps.lod_neigh_bias).getNorm2(),
-              aps.num_pred_nearest_neighbours,
+              pointIndex1, d2, aps.num_pred_nearest_neighbours,
               startIndex + k - i + 2 * aps.search_range);
           }
         }
@@ -776,9 +780,11 @@ subsampleByDistance(
         continue;
       }
       const auto& point = pointCloud[packedVoxel[index].index];
+      const auto& lastRetained =
+        pointCloud[packedVoxel[retained.back()].index];
+      auto distanceToLastRetained = (lastRetained - point).getNorm2<double>();
       if (
-        (pointCloud[packedVoxel[retained.back()].index] - point).getNorm2()
-          <= radius2
+        distanceToLastRetained <= radius2
         || FindNeighborWithinDistance(
              pointCloud, packedVoxel, index, radius2, searchRange, retained)
           != PCC_UNDEFINED_INDEX) {

tmc3/pointset_processing.cpp

@@ -294,7 +294,7 @@ recolourColour(
       double maxAttributeDist2 = std::numeric_limits<double>::min();
       for (int i = 0; i < nNN; ++i) {
         for (int j = 0; j < nNN; ++j) {
-          const double dist2 = (colors[i] - colors[j]).getNorm2();
+          const double dist2 = (colors[i] - colors[j]).getNorm2<double>();
           if (dist2 > maxAttributeDist2) {
             maxAttributeDist2 = dist2;
           }
@@ -416,7 +416,7 @@ recolourColour(
           double maxAttributeDist2 = std::numeric_limits<double>::min();
           for (int i = 0; i < nNN; ++i) {
             for (int j = 0; j < nNN; ++j) {
-              const double dist2 = (colors[i] - colors[j]).getNorm2();
+              const double dist2 = (colors[i] - colors[j]).getNorm2<double>();
               if (dist2 > maxAttributeDist2) {
                 maxAttributeDist2 = dist2;
               }
@@ -462,7 +462,7 @@ recolourColour(
       }
     }
     const double r = double(pointCountTarget) / double(pointCountSource);
-    const double delta2 = (centroid2 - centroid1).getNorm2();
+    const double delta2 = (centroid2 - centroid1).getNorm2<double>();
     const double eps = 0.000001;
 
     const bool fixWeight = 1;  // m42538