attr/m47352: permit partial geometry and attribute reconstruction

This provides a means to achieve spatial scalability for a G-PCC
bitstream.  For bitstreams encoded with aps_scalable_enabled_flag set,
the decoder may partially decode the geometry octree
(minGeomNodeSizeLog2) and the corresponding point attributes.  The
functionality is achieved by constraining the LoD generation process to
align with partial geometry decoding.

doc/README.options.md

@@ -72,8 +72,11 @@ conversion process is lossless.
 
 Decoder-specific options
 ========================
-There are no decoder specific options at this time.
 
+### `--minGeomNodeSizeLog2=INT-VALUE`
+The option indicates the number of skipped lod layers from leaf lod.
+If aps.scalable_enable_flag is 1, the option is valid.
+Otherwise, the option is ignored.
 
 Encoder-specific options
 ========================
@@ -287,6 +290,11 @@ Controls the level-of-detail generation method:
 Part of LoD attribute coding.  Permits (1) points to be predicted from
 previously reconstructed points within the same LoD.
 
+### `--aps_scalable_enabled_flag=0|1`
+Enable spatially scalable attribute encoding.
+The option is only valid when `transformType=2`, `lodDecimation=0`,
+and `trisoup_node_size_log2=0`.
+
 ### `--levelOfDetailCount=INT-VALUE`
 Attribute's number of levels of detail.
 

tmc3/AttributeDecoder.cpp

@@ -179,6 +179,8 @@ AttributeDecoder::decode(
   const SequenceParameterSet& sps,
   const AttributeDescription& attr_desc,
   const AttributeParameterSet& attr_aps,
+  int geom_num_points,
+  int minGeomNodeSizeLog2,
   const PayloadBuffer& payload,
   PCCPointSet3& pointCloud)
 {
@@ -204,7 +206,8 @@ AttributeDecoder::decode(
 
     case AttributeEncoding::kLiftingTransform:
       decodeReflectancesLift(
-        attr_desc, attr_aps, quantLayers, decoder, pointCloud);
+        attr_desc, attr_aps, quantLayers, geom_num_points, minGeomNodeSizeLog2,
+        decoder, pointCloud);
       break;
     }
   } else if (attr_desc.attr_num_dimensions == 3) {
@@ -218,7 +221,9 @@ AttributeDecoder::decode(
       break;
 
     case AttributeEncoding::kLiftingTransform:
-      decodeColorsLift(attr_desc, attr_aps, quantLayers, decoder, pointCloud);
+      decodeColorsLift(
+        attr_desc, attr_aps, quantLayers, geom_num_points, minGeomNodeSizeLog2,
+        decoder, pointCloud);
       break;
     }
   } else {
@@ -281,9 +286,10 @@ AttributeDecoder::decodeReflectancesPred(
 
   buildPredictorsFast(
     pointCloud, aps.lod_decimation_enabled_flag,
-    aps.intra_lod_prediction_enabled_flag, aps.dist2, aps.num_detail_levels,
-    aps.num_pred_nearest_neighbours, aps.search_range, aps.search_range,
-    predictors, numberOfPointsPerLOD, indexesLOD);
+    aps.scalable_lifting_enabled_flag, aps.intra_lod_prediction_enabled_flag,
+    aps.dist2, aps.num_detail_levels, aps.num_pred_nearest_neighbours,
+    aps.search_range, aps.search_range, 0, predictors, numberOfPointsPerLOD,
+    indexesLOD);
 
   const int64_t maxReflectance = (1ll << desc.attr_bitdepth) - 1;
   int zero_cnt = decoder.decodeZeroCnt(pointCount);
@@ -372,9 +378,10 @@ AttributeDecoder::decodeColorsPred(
 
   buildPredictorsFast(
     pointCloud, aps.lod_decimation_enabled_flag,
-    aps.intra_lod_prediction_enabled_flag, aps.dist2, aps.num_detail_levels,
-    aps.num_pred_nearest_neighbours, aps.search_range, aps.search_range,
-    predictors, numberOfPointsPerLOD, indexesLOD);
+    aps.scalable_lifting_enabled_flag, aps.intra_lod_prediction_enabled_flag,
+    aps.dist2, aps.num_detail_levels, aps.num_pred_nearest_neighbours,
+    aps.search_range, aps.search_range, 0, predictors, numberOfPointsPerLOD,
+    indexesLOD);
 
   uint32_t values[3];
   int zero_cnt = decoder.decodeZeroCnt(pointCount);
@@ -552,6 +559,8 @@ AttributeDecoder::decodeColorsLift(
   const AttributeDescription& desc,
   const AttributeParameterSet& aps,
   const std::vector<Quantizers>& quantLayers,
+  int geom_num_points,
+  int minGeomNodeSizeLog2,
   PCCResidualsDecoder& decoder,
   PCCPointSet3& pointCloud)
 {
@@ -560,24 +569,41 @@ AttributeDecoder::decodeColorsLift(
   std::vector<uint32_t> numberOfPointsPerLOD;
   std::vector<uint32_t> indexesLOD;
 
+  if (minGeomNodeSizeLog2 > 0)
+    assert(aps.scalable_lifting_enabled_flag);
+
   buildPredictorsFast(
     pointCloud, aps.lod_decimation_enabled_flag,
-    aps.intra_lod_prediction_enabled_flag, aps.dist2, aps.num_detail_levels,
-    aps.num_pred_nearest_neighbours, aps.search_range, aps.search_range,
-    predictors, numberOfPointsPerLOD, indexesLOD);
+    aps.scalable_lifting_enabled_flag, aps.intra_lod_prediction_enabled_flag,
+    aps.dist2, aps.num_detail_levels, aps.num_pred_nearest_neighbours,
+    aps.search_range, aps.search_range, minGeomNodeSizeLog2, predictors,
+    numberOfPointsPerLOD, indexesLOD);
 
   for (size_t predictorIndex = 0; predictorIndex < pointCount;
        ++predictorIndex) {
     predictors[predictorIndex].computeWeights();
   }
   std::vector<uint64_t> weights;
+  if (!aps.scalable_lifting_enabled_flag) {
     PCCComputeQuantizationWeights(predictors, weights);
+  } else {
+    computeQuantizationWeightsScalable(
+      predictors, numberOfPointsPerLOD, geom_num_points, minGeomNodeSizeLog2,
+      weights);
+  }
+
   const size_t lodCount = numberOfPointsPerLOD.size();
   std::vector<Vec3<int64_t>> colors;
   colors.resize(pointCount);
 
+  // NB: when partially decoding, the truncated unary limit for zero_run
+  // must be the original value.  geom_num_points may be the case.  However,
+  // the encoder does lie sometimes, and there are actually more points.
+  int zeroCntLimit = std::max(geom_num_points, int(pointCount));
+
   // decompress
-  int zero_cnt = decoder.decodeZeroCnt(pointCount);
+  int zero_cnt = decoder.decodeZeroCnt(zeroCntLimit);
+  std::cout << zero_cnt << '\n';
   int quantLayer = 0;
   for (size_t predictorIndex = 0; predictorIndex < pointCount;
        ++predictorIndex) {
@@ -592,7 +618,8 @@ AttributeDecoder::decodeColorsLift(
       zero_cnt--;
     } else {
       decoder.decode(values);
-      zero_cnt = decoder.decodeZeroCnt(pointCount);
+      zero_cnt = decoder.decodeZeroCnt(zeroCntLimit);
+  std::cout << zero_cnt << '\n';
     }
 
     const int64_t quantWeight = weights[predictorIndex];
@@ -634,6 +661,8 @@ AttributeDecoder::decodeReflectancesLift(
   const AttributeDescription& desc,
   const AttributeParameterSet& aps,
   const std::vector<Quantizers>& quantLayers,
+  int geom_num_points,
+  int minGeomNodeSizeLog2,
   PCCResidualsDecoder& decoder,
   PCCPointSet3& pointCloud)
 {
@@ -642,11 +671,15 @@ AttributeDecoder::decodeReflectancesLift(
   std::vector<uint32_t> numberOfPointsPerLOD;
   std::vector<uint32_t> indexesLOD;
 
+  if (minGeomNodeSizeLog2 > 0)
+    assert(aps.scalable_lifting_enabled_flag);
+
   buildPredictorsFast(
     pointCloud, aps.lod_decimation_enabled_flag,
-    aps.intra_lod_prediction_enabled_flag, aps.dist2, aps.num_detail_levels,
-    aps.num_pred_nearest_neighbours, aps.search_range, aps.search_range,
-    predictors, numberOfPointsPerLOD, indexesLOD);
+    aps.scalable_lifting_enabled_flag, aps.intra_lod_prediction_enabled_flag,
+    aps.dist2, aps.num_detail_levels, aps.num_pred_nearest_neighbours,
+    aps.search_range, aps.search_range, minGeomNodeSizeLog2, predictors,
+    numberOfPointsPerLOD, indexesLOD);
 
   for (size_t predictorIndex = 0; predictorIndex < pointCount;
        ++predictorIndex) {
@@ -654,13 +687,25 @@ AttributeDecoder::decodeReflectancesLift(
   }
 
   std::vector<uint64_t> weights;
+  if (!aps.scalable_lifting_enabled_flag) {
     PCCComputeQuantizationWeights(predictors, weights);
+  } else {
+    computeQuantizationWeightsScalable(
+      predictors, numberOfPointsPerLOD, geom_num_points, minGeomNodeSizeLog2,
+      weights);
+  }
+
   const size_t lodCount = numberOfPointsPerLOD.size();
   std::vector<int64_t> reflectances;
   reflectances.resize(pointCount);
 
+  // NB: when partially decoding, the truncated unary limit for zero_run
+  // must be the original value.  geom_num_points may be the case.  However,
+  // the encoder does lie sometimes, and there are actually more points.
+  int zeroCntLimit = std::max(geom_num_points, int(pointCount));
+
   // decompress
-  int zero_cnt = decoder.decodeZeroCnt(pointCount);
+  int zero_cnt = decoder.decodeZeroCnt(zeroCntLimit);
   int quantLayer = 0;
   for (size_t predictorIndex = 0; predictorIndex < pointCount;
        ++predictorIndex) {
@@ -674,7 +719,7 @@ AttributeDecoder::decodeReflectancesLift(
       zero_cnt--;
     } else {
       detail = decoder.decode();
-      zero_cnt = decoder.decodeZeroCnt(pointCount);
+      zero_cnt = decoder.decodeZeroCnt(zeroCntLimit);
     }
     const int64_t quantWeight = weights[predictorIndex];
     auto& reflectance = reflectances[predictorIndex];

tmc3/AttributeDecoder.h

@@ -56,6 +56,8 @@ public:
     const SequenceParameterSet& sps,
     const AttributeDescription& desc,
     const AttributeParameterSet& aps,
+    int geom_num_points,
+    int minGeomNodeSizeLog2,
     const PayloadBuffer&,
     PCCPointSet3& pointCloud);
 
@@ -66,6 +68,8 @@ protected:
     const AttributeDescription& desc,
     const AttributeParameterSet& aps,
     const std::vector<Quantizers>& quantLayers,
+    int geom_num_points,
+    int minGeomNodeSizeLog2,
     PCCResidualsDecoder& decoder,
     PCCPointSet3& pointCloud);
 
@@ -73,6 +77,8 @@ protected:
     const AttributeDescription& desc,
     const AttributeParameterSet& aps,
     const std::vector<Quantizers>& quantLayers,
+    int geom_num_points,
+    int minGeomNodeSizeLog2,
     PCCResidualsDecoder& decoder,
     PCCPointSet3& pointCloud);
 

tmc3/AttributeEncoder.cpp

@@ -481,9 +481,10 @@ AttributeEncoder::encodeReflectancesPred(
 
   buildPredictorsFast(
     pointCloud, aps.lod_decimation_enabled_flag,
-    aps.intra_lod_prediction_enabled_flag, aps.dist2, aps.num_detail_levels,
-    aps.num_pred_nearest_neighbours, aps.search_range, aps.search_range,
-    predictors, numberOfPointsPerLOD, indexesLOD);
+    aps.scalable_lifting_enabled_flag, aps.intra_lod_prediction_enabled_flag,
+    aps.dist2, aps.num_detail_levels, aps.num_pred_nearest_neighbours,
+    aps.search_range, aps.search_range, 0, predictors, numberOfPointsPerLOD,
+    indexesLOD);
 
   const int64_t clipMax = (1ll << desc.attr_bitdepth) - 1;
   PCCResidualsEntropyEstimator context;
@@ -671,9 +672,10 @@ AttributeEncoder::encodeColorsPred(
 
   buildPredictorsFast(
     pointCloud, aps.lod_decimation_enabled_flag,
-    aps.intra_lod_prediction_enabled_flag, aps.dist2, aps.num_detail_levels,
-    aps.num_pred_nearest_neighbours, aps.search_range, aps.search_range,
-    predictors, numberOfPointsPerLOD, indexesLOD);
+    aps.scalable_lifting_enabled_flag, aps.intra_lod_prediction_enabled_flag,
+    aps.dist2, aps.num_detail_levels, aps.num_pred_nearest_neighbours,
+    aps.search_range, aps.search_range, 0, predictors, numberOfPointsPerLOD,
+    indexesLOD);
 
   const int64_t clipMax = (1ll << desc.attr_bitdepth) - 1;
   uint32_t values[3];
@@ -924,16 +926,22 @@ AttributeEncoder::encodeColorsLift(
 
   buildPredictorsFast(
     pointCloud, aps.lod_decimation_enabled_flag,
-    aps.intra_lod_prediction_enabled_flag, aps.dist2, aps.num_detail_levels,
-    aps.num_pred_nearest_neighbours, aps.search_range, aps.search_range,
-    predictors, numberOfPointsPerLOD, indexesLOD);
+    aps.scalable_lifting_enabled_flag, aps.intra_lod_prediction_enabled_flag,
+    aps.dist2, aps.num_detail_levels, aps.num_pred_nearest_neighbours,
+    aps.search_range, aps.search_range, 0, predictors, numberOfPointsPerLOD,
+    indexesLOD);
 
   for (size_t predictorIndex = 0; predictorIndex < pointCount;
        ++predictorIndex) {
     predictors[predictorIndex].computeWeights();
   }
   std::vector<uint64_t> weights;
+  if (!aps.scalable_lifting_enabled_flag) {
     PCCComputeQuantizationWeights(predictors, weights);
+  } else {
+    computeQuantizationWeightsScalable(
+      predictors, numberOfPointsPerLOD, pointCount, 0, weights);
+  }
   const size_t lodCount = numberOfPointsPerLOD.size();
   std::vector<Vec3<int64_t>> colors;
   colors.resize(pointCount);
@@ -1027,16 +1035,22 @@ AttributeEncoder::encodeReflectancesLift(
 
   buildPredictorsFast(
     pointCloud, aps.lod_decimation_enabled_flag,
-    aps.intra_lod_prediction_enabled_flag, aps.dist2, aps.num_detail_levels,
-    aps.num_pred_nearest_neighbours, aps.search_range, aps.search_range,
-    predictors, numberOfPointsPerLOD, indexesLOD);
+    aps.scalable_lifting_enabled_flag, aps.intra_lod_prediction_enabled_flag,
+    aps.dist2, aps.num_detail_levels, aps.num_pred_nearest_neighbours,
+    aps.search_range, aps.search_range, 0, predictors, numberOfPointsPerLOD,
+    indexesLOD);
 
   for (size_t predictorIndex = 0; predictorIndex < pointCount;
        ++predictorIndex) {
     predictors[predictorIndex].computeWeights();
   }
   std::vector<uint64_t> weights;
+  if (!aps.scalable_lifting_enabled_flag) {
     PCCComputeQuantizationWeights(predictors, weights);
+  } else {
+    computeQuantizationWeightsScalable(
+      predictors, numberOfPointsPerLOD, pointCount, 0, weights);
+  }
 
   const size_t lodCount = numberOfPointsPerLOD.size();
   std::vector<int64_t> reflectances;

tmc3/PCCPointSet.h

@@ -371,6 +371,7 @@ public:
       frameidx.resize(size);
     }
   }
+
   void reserve(const size_t size)
   {
     positions.reserve(size);
@@ -392,6 +393,23 @@ public:
     frameidx.clear();
   }
 
+  size_t removeDuplicatePointInQuantizedPoint(int minGeomNodeSizeLog2)
+  {
+    for (int i = 0; i < positions.size(); i++) {
+      Vec3<double> newPoint = positions[i];
+      if (minGeomNodeSizeLog2 > 0) {
+        uint32_t mask = ((uint32_t)-1) << minGeomNodeSizeLog2;
+        positions[i].x() = ((int32_t)(positions[i].x()) & mask);
+        positions[i].y() = ((int32_t)(positions[i].y()) & mask);
+        positions[i].z() = ((int32_t)(positions[i].z()) & mask);
+      }
+    }
+    positions.erase(
+      std::unique(positions.begin(), positions.end()), positions.end());
+
+    return positions.size();
+  }
+
   void append(const PCCPointSet3& src)
   {
     if (!getPointCount())

tmc3/PCCTMC3Common.h

@@ -431,6 +431,47 @@ PCCComputeQuantizationWeights(
 
 //---------------------------------------------------------------------------
 
+inline void
+computeQuantizationWeightsScalable(
+  const std::vector<PCCPredictor>& predictors,
+  const std::vector<uint32_t>& numberOfPointsPerLOD,
+  size_t geom_num_points,
+  int32_t minGeomNodeSizeLog2,
+  std::vector<uint64_t>& quantizationWeights)
+{
+  const size_t pointCount = predictors.size();
+  quantizationWeights.resize(pointCount);
+  for (size_t i = 0; i < pointCount; ++i) {
+    quantizationWeights[i] = (1 << kFixedPointWeightShift);
+  }
+
+  const size_t lodCount = numberOfPointsPerLOD.size();
+  for (size_t lodIndex = 0; lodIndex < lodCount; ++lodIndex) {
+    const size_t startIndex =
+      (lodIndex == 0) ? 0 : numberOfPointsPerLOD[lodIndex - 1];
+    const size_t endIndex = numberOfPointsPerLOD[lodIndex];
+
+    const size_t predictorCount = endIndex - startIndex;
+    for (size_t index = 0; index < predictorCount; ++index) {
+      const size_t predictorIndex = index + startIndex;
+      const double currentQuantWeight = (geom_num_points / predictorCount);
+
+      if (!minGeomNodeSizeLog2 && (lodIndex == lodCount - 1)) {
+        quantizationWeights[predictorIndex] = (1 << kFixedPointWeightShift);
+      } else {
+        quantizationWeights[predictorIndex] =
+          currentQuantWeight * (1 << kFixedPointWeightShift);
+      }
+    }
+  }
+
+  for (auto& w : quantizationWeights) {
+    w = isqrt(w);
+  }
+}
+
+//---------------------------------------------------------------------------
+
 inline uint32_t
 FindNeighborWithinDistance(
   const PCCPointSet3& pointCloud,
@@ -459,6 +500,80 @@ FindNeighborWithinDistance(
 
 //---------------------------------------------------------------------------
 
+inline bool
+checkDistance(
+  const Vec3<double>& point,
+  const Vec3<double>& refpoint,
+  int32_t nodeSizeLog2)
+{
+  uint32_t mask = uint32_t(-1) << nodeSizeLog2;
+
+  int32_t minX = (int32_t(refpoint.x()) & mask) * 2 - 1;
+  int32_t minY = (int32_t(refpoint.y()) & mask) * 2 - 1;
+  int32_t minZ = (int32_t(refpoint.z()) & mask) * 2 - 1;
+
+  int32_t maxX = minX + (2 << nodeSizeLog2);
+  int32_t maxY = minY + (2 << nodeSizeLog2);
+  int32_t maxZ = minZ + (2 << nodeSizeLog2);
+
+  int32_t x = int32_t(point.x()) * 2;
+  int32_t y = int32_t(point.y()) * 2;
+  int32_t z = int32_t(point.z()) * 2;
+
+  return minX < x && maxX > x && minY < y && maxY > y && minZ < z && maxZ > z;
+}
+
+//---------------------------------------------------------------------------
+
+inline uint32_t
+findNeighborWithinVoxel(
+  const PCCPointSet3& pointCloud,
+  const std::vector<MortonCodeWithIndex>& packedVoxel,
+  int32_t index,
+  int32_t octreeNodeSizeLog2,
+  int32_t searchRange,
+  std::vector<uint32_t>& retained)
+{
+  const auto& point = pointCloud[packedVoxel[index].index];
+  int32_t retainedSize = retained.size();
+  int32_t j = retainedSize - 1;
+  int32_t k = 0;
+  while (j >= 0 && ++k < searchRange) {
+    int32_t index1 = retained[j];
+    int32_t pointIndex1 = packedVoxel[index1].index;
+    const auto& point1 = pointCloud[pointIndex1];
+
+    if (checkDistance(point, point1, octreeNodeSizeLog2))
+      return index1;
+    --j;
+  }
+  return PCC_UNDEFINED_INDEX;
+}
+
+//---------------------------------------------------------------------------
+
+inline Vec3<double>
+clacIntermediatePosition(
+  bool enabled, int32_t nodeSizeLog2, const Vec3<double>& point)
+{
+  if (!enabled || !nodeSizeLog2)
+    return point;
+
+  uint32_t mask = (uint32_t(-1)) << nodeSizeLog2;
+  int32_t centerX = (int32_t(point.x()) & mask) + (1 << (nodeSizeLog2 - 1));
+  int32_t centerY = (int32_t(point.y()) & mask) + (1 << (nodeSizeLog2 - 1));
+  int32_t centerZ = (int32_t(point.z()) & mask) + (1 << (nodeSizeLog2 - 1));
+
+  Vec3<double> newPoint = point;
+  newPoint.x() = centerX;
+  newPoint.y() = centerY;
+  newPoint.z() = centerZ;
+
+  return newPoint;
+}
+
+//---------------------------------------------------------------------------
+
 inline int
 indexTieBreaker(int a, int b)
 {
@@ -476,22 +591,27 @@ computeNearestNeighbors(
   const int32_t endIndex,
   const int32_t searchRange,
   const int32_t numberOfNearestNeighborsInPrediction,
+  const bool intraLodPredictionEnabled,
+  const bool scalableEnableFlag,
+  const int32_t nodeSizeLog2,
   std::vector<uint32_t>& indexes,
   std::vector<PCCPredictor>& predictors,
   std::vector<uint32_t>& pointIndexToPredictorIndex,
   int32_t& predIndex,
-  std::vector<Box3<double>>& bBoxes,
-  const bool intraLodPredictionEnabled)
+  std::vector<Box3<double>>& bBoxes)
 {
   const int32_t retainedSize = retained.size();
   const int32_t bucketSize = 8;
   bBoxes.resize((retainedSize + bucketSize - 1) / bucketSize);
   for (int32_t i = 0, b = 0; i < retainedSize; ++b) {
     auto& bBox = bBoxes[b];
-    bBox.min = bBox.max = pointCloud[packedVoxel[retained[i++]].index];
+    bBox.min = bBox.max = clacIntermediatePosition(
+      scalableEnableFlag, nodeSizeLog2,
+      pointCloud[packedVoxel[retained[i++]].index]);
     for (int32_t k = 1; k < bucketSize && i < retainedSize; ++k, ++i) {
       const int32_t pointIndex = packedVoxel[retained[i]].index;
-      const auto& point = pointCloud[pointIndex];
+      const auto point = clacIntermediatePosition(
+        scalableEnableFlag, nodeSizeLog2, pointCloud[pointIndex]);
       for (int32_t p = 0; p < 3; ++p) {
         bBox.min[p] = std::min(bBox.min[p], point[p]);
         bBox.max[p] = std::max(bBox.max[p], point[p]);
@@ -523,7 +643,8 @@ computeNearestNeighbors(
     const int32_t index = indexes[i];
     const int64_t mortonCode = packedVoxel[index].mortonCode;
     const int32_t pointIndex = packedVoxel[index].index;
-    const auto& point = pointCloud[pointIndex];
+    const auto point = clacIntermediatePosition(
+      scalableEnableFlag, nodeSizeLog2, pointCloud[pointIndex]);
     indexes[i] = pointIndex;
     while (j < retainedSize - 1
            && mortonCode >= packedVoxel[retained[j]].mortonCode)
@@ -542,10 +663,17 @@ computeNearestNeighbors(
 
     for (int32_t k = k0; k < k1; ++k) {
       const int32_t pointIndex1 = packedVoxel[retained[k]].index;
-      const auto& point1 = pointCloud[pointIndex1];
+      const auto point1 = clacIntermediatePosition(
+        scalableEnableFlag, nodeSizeLog2, pointCloud[pointIndex1]);
+
+      double norm2 = (point - point1).getNorm2();
+      if (nodeSizeLog2 > 0 && point == point1) {
+        norm2 = double(1 << (nodeSizeLog2 - 1));
+        norm2 = norm2 * norm2;
+      }
       predictor.insertNeighbor(
-        pointIndex1, (point - point1).getNorm2(),
-        numberOfNearestNeighborsInPrediction, indexTieBreaker(k, j));
+        pointIndex1, norm2, numberOfNearestNeighborsInPrediction,
+        indexTieBreaker(k, j));
     }
 
     for (int32_t s0 = 1, sr = (1 + searchRange / bucketSize); s0 < sr; ++s0) {
@@ -563,10 +691,17 @@ computeNearestNeighbors(
           const int32_t k1 = std::min((bucketIndex1 + 1) * bucketSize, j1);
           for (int32_t k = k0; k < k1; ++k) {
             const int32_t pointIndex1 = packedVoxel[retained[k]].index;
-            const auto& point1 = pointCloud[pointIndex1];
+            const auto point1 = clacIntermediatePosition(
+              scalableEnableFlag, nodeSizeLog2, pointCloud[pointIndex1]);
+
+            double norm2 = (point - point1).getNorm2();
+            if (nodeSizeLog2 > 0 && point == point1) {
+              norm2 = (double)(1 << (nodeSizeLog2 - 1));
+              norm2 = norm2 * norm2;
+            }
             predictor.insertNeighbor(
-              pointIndex1, (point - point1).getNorm2(),
-              numberOfNearestNeighborsInPrediction, indexTieBreaker(k, j));
+              pointIndex1, norm2, numberOfNearestNeighborsInPrediction,
+              indexTieBreaker(k, j));
           }
         }
       }
@@ -653,6 +788,46 @@ subsampleByDistance(
 
 //---------------------------------------------------------------------------
 
+inline void
+subsampleByOctree(
+  const PCCPointSet3& pointCloud,
+  const std::vector<MortonCodeWithIndex>& packedVoxel,
+  const std::vector<uint32_t>& input,
+  int32_t octreeNodeSizeLog2,
+  std::vector<uint32_t>& retained,
+  std::vector<uint32_t>& indexes)
+{
+  int32_t nodeSizeLog2OfUpperLayer = octreeNodeSizeLog2 + 1;
+  int32_t searchRange = 1;
+
+  if (input.size() == 1) {
+    indexes.push_back(input[0]);
+  } else {
+    for (const auto index : input) {
+      if (retained.empty()) {
+        retained.push_back(index);
+        continue;
+      }
+      const auto& point = pointCloud[packedVoxel[index].index];
+      const auto& retainedPoint =
+        pointCloud[packedVoxel[retained.back()].index];
+
+      if (
+        checkDistance(point, retainedPoint, nodeSizeLog2OfUpperLayer)
+        || findNeighborWithinVoxel(
+             pointCloud, packedVoxel, index, nodeSizeLog2OfUpperLayer,
+             searchRange, retained)
+          != PCC_UNDEFINED_INDEX) {
+        indexes.push_back(index);
+      } else {
+        retained.push_back(index);
+      }
+    }
+  }
+}
+
+//---------------------------------------------------------------------------
+
 inline void
 subsampleByDecimation(
   const std::vector<uint32_t>& input,
@@ -674,19 +849,25 @@ subsampleByDecimation(
 inline void