AttributeEncoder.cpp 35.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
/* The copyright in this software is being made available under the BSD
 * Licence, included below.  This software may be subject to other third
 * party and contributor rights, including patent rights, and no such
 * rights are granted under this licence.
 *
 * Copyright (c) 2017-2018, ISO/IEC
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * * Redistributions of source code must retain the above copyright
 *   notice, this list of conditions and the following disclaimer.
 *
 * * Redistributions in binary form must reproduce the above copyright
 *   notice, this list of conditions and the following disclaimer in the
 *   documentation and/or other materials provided with the distribution.
 *
 * * Neither the name of the ISO/IEC nor the names of its contributors
 *   may be used to endorse or promote products derived from this
 *   software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */

#include "AttributeEncoder.h"
37

David Flynn's avatar
David Flynn committed
38
#include "ArithmeticCodec.h"
39
#include "DualLutCoder.h"
40
#include "constants.h"
41
#include "entropy.h"
42
#include "quantization.h"
43
#include "RAHT.h"
44
#include "FixedPoint.h"
45

46
// todo(df): promote to per-attribute encoder parameter
47
48
static const double kAttrPredLambdaR = 0.01;
static const double kAttrPredLambdaC = 0.01;
49

50
51
52
53
54
namespace pcc {
//============================================================================
// An encapsulation of the entropy coding methods used in attribute coding

struct PCCResidualsEncoder {
55
56
57
58
  EntropyEncoder arithmeticEncoder;
  StaticBitModel binaryModel0;
  AdaptiveBitModel binaryModelDiff[7];
  AdaptiveBitModel binaryModelIsZero[7];
59
  AdaptiveBitModel ctxPredMode[2];
60
61
  AdaptiveBitModel ctxZeroCnt[3];
  AdaptiveBitModel binaryModelIsOne[7];
62
  DualLutCoder<false> symbolCoder[2];
63

64
  void start(const SequenceParameterSet& sps, int numPoints);
65
  int stop();
66
  void encodePredMode(int value, int max);
67
  void encodeZeroCnt(int value, int max);
68
69
70
  void encodeSymbol(uint32_t value, int k1, int k2);
  void encode(uint32_t value0, uint32_t value1, uint32_t value2);
  void encode(uint32_t value);
71
72
73
74
};

//----------------------------------------------------------------------------

75
void
76
PCCResidualsEncoder::start(const SequenceParameterSet& sps, int pointCount)
77
{
78
79
  // todo(df): remove estimate when arithmetic codec is replaced
  int maxAcBufLen = pointCount * 3 * 2 + 1024;
80
  arithmeticEncoder.setBuffer(maxAcBufLen, nullptr);
81
  arithmeticEncoder.enableBypassStream(sps.cabac_bypass_stream_enabled_flag);
82
  arithmeticEncoder.start();
83
84
85
86
}

//----------------------------------------------------------------------------

87
int
88
89
PCCResidualsEncoder::stop()
{
90
  return arithmeticEncoder.stop();
91
92
93
94
}

//----------------------------------------------------------------------------

95
void
96
PCCResidualsEncoder::encodePredMode(int mode, int maxMode)
97
{
98
99
100
101
102
103
104
105
106
107
108
109
110
  // max = 0 => no direct predictors are used
  if (maxMode == 0)
    return;

  int ctxIdx = 0;
  for (int i = 0; i < mode; i++) {
    arithmeticEncoder.encode(1, ctxPredMode[ctxIdx]);
    ctxIdx = 1;
  }

  // Truncated unary
  if (mode != maxMode)
    arithmeticEncoder.encode(0, ctxPredMode[ctxIdx]);
111
112
113
114
}

//----------------------------------------------------------------------------

115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
void
PCCResidualsEncoder::encodeZeroCnt(int mode, int maxMode)
{
  // max = 0 => no direct predictors are used
  if (maxMode == 0)
    return;

  int ctxIdx = 0;
  for (int i = 0; i < mode; i++) {
    arithmeticEncoder.encode(1, ctxZeroCnt[ctxIdx]);
    ctxIdx = (ctxIdx == 0 ? 1 : 2);
  }

  // Truncated unary
  if (mode != maxMode)
    arithmeticEncoder.encode(0, ctxZeroCnt[ctxIdx]);
}

//----------------------------------------------------------------------------

135
136
void
PCCResidualsEncoder::encodeSymbol(uint32_t value, int k1, int k2)
137
{
138
139
140
141
142
  bool isZero = value == 0;
  arithmeticEncoder.encode(isZero, binaryModelIsZero[k1]);
  if (isZero) {
    return;
  }
143
144
145
146
147
148
149
150

  bool isOne = value == 1;
  arithmeticEncoder.encode(isOne, binaryModelIsOne[k1]);
  if (isOne) {
    return;
  }
  value -= 2;

151
  if (value < kAttributeResidualAlphabetSize) {
152
    symbolCoder[k2].encode(value, &arithmeticEncoder);
153
  } else {
154
    int alphabetSize = kAttributeResidualAlphabetSize;
155
    symbolCoder[k2].encode(alphabetSize, &arithmeticEncoder);
156
    arithmeticEncoder.encodeExpGolomb(
157
      value - alphabetSize, 0, binaryModel0, binaryModelDiff[k1]);
158
159
160
  }
}

161
162
163
//----------------------------------------------------------------------------

void
164
PCCResidualsEncoder::encode(uint32_t value0, uint32_t value1, uint32_t value2)
165
{
166
167
168
169
170
171
  if (value0 == value1 && value0 == value2) {
    value0--;
    value1--;
    value2--;
  }

172
173
174
175
176
177
178
179
180
181
182
183
  int b0 = value0 == 0;
  int b1 = value1 == 0;
  encodeSymbol(value0, 0, 0);
  encodeSymbol(value1, 1 + b0, 1);
  encodeSymbol(value2, 3 + (b0 << 1) + b1, 1);
}

//----------------------------------------------------------------------------

void
PCCResidualsEncoder::encode(uint32_t value)
{
184
  encodeSymbol(value - 1, 0, 0);
185
186
}

187
//============================================================================
188
189
190
// An encapsulation of the entropy coding methods used in attribute coding

struct PCCResidualsEntropyEstimator {
191
192
  size_t freq0[kAttributeResidualAlphabetSize + 1];
  size_t freq1[kAttributeResidualAlphabetSize + 1];
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
  size_t symbolCount0;
  size_t symbolCount1;
  size_t isZero0Count;
  size_t isZero1Count;
  PCCResidualsEntropyEstimator() { init(); }
  void init();
  double bitsDetail(
    const uint32_t detail,
    const size_t symbolCount,
    const size_t* const freq) const;
  double bits(const uint32_t value0) const;
  void update(const uint32_t value0);
  double bits(
    const uint32_t value0, const uint32_t value1, const uint32_t value2) const;
  void
  update(const uint32_t value0, const uint32_t value1, const uint32_t value2);
};

//----------------------------------------------------------------------------
212

213
void
214
PCCResidualsEntropyEstimator::init()
215
{
216
  for (size_t i = 0; i <= kAttributeResidualAlphabetSize; ++i) {
217
218
219
    freq0[i] = 1;
    freq1[i] = 1;
  }
220
221
  symbolCount0 = kAttributeResidualAlphabetSize + 1;
  symbolCount1 = kAttributeResidualAlphabetSize + 1;
222
223
  isZero1Count = isZero0Count = symbolCount0 / 2;
}
224

225
226
227
228
229
230
231
232
233
//----------------------------------------------------------------------------

double
PCCResidualsEntropyEstimator::bitsDetail(
  const uint32_t detail,
  const size_t symbolCount,
  const size_t* const freq) const
{
  const uint32_t detailClipped =
234
    std::min(detail, uint32_t(kAttributeResidualAlphabetSize));
235
236
237
  const double pDetail =
    PCCClip(double(freq[detailClipped]) / symbolCount, 0.001, 0.999);
  double bits = -log2(pDetail);
238
239
  if (detail >= kAttributeResidualAlphabetSize) {
    const double x = double(detail) - double(kAttributeResidualAlphabetSize);
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
    bits += 2.0 * std::floor(log2(x + 1.0)) + 1.0;
  }
  return bits;
}

//----------------------------------------------------------------------------

double
PCCResidualsEntropyEstimator::bits(const uint32_t value0) const
{
  const bool isZero0 = value0 == 0;
  const double pIsZero0 = isZero0
    ? double(isZero0Count) / symbolCount0
    : double(symbolCount0 - isZero0Count) / symbolCount0;
  double bits = -log2(PCCClip(pIsZero0, 0.001, 0.999));
  if (!isZero0) {
    bits += bitsDetail(value0 - 1, symbolCount0, freq0);
257
  }
258
  return bits;
259
260
261
262
}

//----------------------------------------------------------------------------

263
264
265
266
267
268
void
PCCResidualsEntropyEstimator::update(const uint32_t value0)
{
  const bool isZero0 = value0 == 0;
  ++symbolCount0;
  if (!isZero0) {
269
    ++freq0[std::min(value0 - 1, uint32_t(kAttributeResidualAlphabetSize))];
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
  } else {
    ++isZero0Count;
  }
}

//----------------------------------------------------------------------------

double
PCCResidualsEntropyEstimator::bits(
  const uint32_t value0, const uint32_t value1, const uint32_t value2) const
{
  const bool isZero0 = value0 == 0;
  const double pIsZero0 = isZero0
    ? double(isZero0Count) / symbolCount0
    : double(symbolCount0 - isZero0Count) / symbolCount0;
  double bits = -log2(PCCClip(pIsZero0, 0.001, 0.999));
  if (!isZero0) {
    bits += bitsDetail(value0 - 1, symbolCount0, freq0);
  }

  const bool isZero1 = value1 == 0 && value2 == 0;
  const double pIsZero1 = isZero1
    ? double(isZero1Count) / symbolCount0
    : double(symbolCount0 - isZero1Count) / symbolCount0;
  bits -= log2(PCCClip(pIsZero1, 0.001, 0.999));
  if (!isZero1) {
    bits += bitsDetail(value1, symbolCount1, freq1);
    bits += bitsDetail(value2, symbolCount1, freq1);
  }
  return bits;
}

//----------------------------------------------------------------------------

void
PCCResidualsEntropyEstimator::update(
  const uint32_t value0, const uint32_t value1, const uint32_t value2)
{
  const bool isZero0 = value0 == 0;
  ++symbolCount0;
  if (!isZero0) {
311
    ++freq0[std::min(value0 - 1, uint32_t(kAttributeResidualAlphabetSize))];
312
313
314
315
316
317
318
  } else {
    ++isZero0Count;
  }

  const bool isZero1 = value1 == 0 && value2 == 0;
  symbolCount1 += 2;
  if (!isZero1) {
319
320
    ++freq1[std::min(value1, uint32_t(kAttributeResidualAlphabetSize))];
    ++freq1[std::min(value2, uint32_t(kAttributeResidualAlphabetSize))];
321
322
323
324
325
326
327
328
  } else {
    ++isZero1Count;
  }
}

//============================================================================
// AttributeEncoder Members

329
void
330
AttributeEncoder::encode(
331
  const SequenceParameterSet& sps,
332
  const AttributeDescription& desc,
333
  const AttributeParameterSet& attr_aps,
334
  const AttributeBrickHeader& abh,
335
  PCCPointSet3& pointCloud,
336
  PayloadBuffer* payload)
337
{
338
  std::vector<Quantizers> quantLayers = deriveQuantizerLayers(attr_aps, abh);
339

340
  PCCResidualsEncoder encoder;
341
  encoder.start(sps, int(pointCloud.getPointCount()));
342

343
  if (desc.attr_num_dimensions == 1) {
344
345
    switch (attr_aps.attr_encoding) {
    case AttributeEncoding::kRAHTransform:
346
      encodeReflectancesTransformRaht(
347
        desc, attr_aps, quantLayers, pointCloud, encoder);
348
      break;
349

350
    case AttributeEncoding::kPredictingTransform:
351
      encodeReflectancesPred(desc, attr_aps, quantLayers, pointCloud, encoder);
352
      break;
353

354
    case AttributeEncoding::kLiftingTransform:
355
      encodeReflectancesLift(desc, attr_aps, quantLayers, pointCloud, encoder);
356
357
      break;
    }
358
  } else if (desc.attr_num_dimensions == 3) {
359
360
    switch (attr_aps.attr_encoding) {
    case AttributeEncoding::kRAHTransform:
361
362
      encodeColorsTransformRaht(
        desc, attr_aps, quantLayers, pointCloud, encoder);
363
364
365
      break;

    case AttributeEncoding::kPredictingTransform:
366
      encodeColorsPred(desc, attr_aps, quantLayers, pointCloud, encoder);
367
368
369
      break;

    case AttributeEncoding::kLiftingTransform:
370
      encodeColorsLift(desc, attr_aps, quantLayers, pointCloud, encoder);
371
372
373
      break;
    }
  } else {
374
    assert(desc.attr_num_dimensions == 1 || desc.attr_num_dimensions == 3);
375
  }
376

377
378
379
380
  uint32_t acDataLen = encoder.stop();
  std::copy_n(
    encoder.arithmeticEncoder.buffer(), acDataLen,
    std::back_inserter(*payload));
381
382
383
384
}

//----------------------------------------------------------------------------

385
386
387
388
int64_t
AttributeEncoder::computeReflectanceResidual(
  const uint64_t reflectance,
  const uint64_t predictedReflectance,
389
  const Quantizer& quant)
390
391
392
{
  const int64_t quantAttValue = reflectance;
  const int64_t quantPredAttValue = predictedReflectance;
393
394
  const int64_t delta = quant.quantize(
    (quantAttValue - quantPredAttValue) << kFixedPointAttributeShift);
395

396
  return IntToUInt(delta);
397
398
399
400
401
402
}

//----------------------------------------------------------------------------

void
AttributeEncoder::computeReflectancePredictionWeights(
403
  const AttributeParameterSet& aps,
404
  const PCCPointSet3& pointCloud,
405
406
  const std::vector<uint32_t>& indexesLOD,
  const uint32_t predictorIndex,
407
408
  PCCPredictor& predictor,
  PCCResidualsEncoder& encoder,
409
  PCCResidualsEntropyEstimator& context,
410
  const Quantizer& quant)
411
412
{
  predictor.computeWeights();
413
  predictor.maxDiff = 0;
414
415
416
  if (predictor.neighborCount > 1) {
    int64_t minValue = 0;
    int64_t maxValue = 0;
417
    for (size_t i = 0; i < predictor.neighborCount; ++i) {
418
      const uint64_t reflectanceNeighbor = pointCloud.getReflectance(
419
        indexesLOD[predictor.neighbors[i].predictorIndex]);
420
421
422
423
424
425
426
427
      if (i == 0 || reflectanceNeighbor < minValue) {
        minValue = reflectanceNeighbor;
      }
      if (i == 0 || reflectanceNeighbor > maxValue) {
        maxValue = reflectanceNeighbor;
      }
    }
    const int64_t maxDiff = maxValue - minValue;
428
    predictor.maxDiff = maxDiff;
429
    if (maxDiff >= aps.adaptive_prediction_threshold) {
430
      uint64_t attrValue =
431
        pointCloud.getReflectance(indexesLOD[predictorIndex]);
432
433
434

      // base case: weighted average of n neighbours
      predictor.predMode = 0;
435
      uint64_t attrPred = predictor.predictReflectance(pointCloud, indexesLOD);
436
      int64_t attrResidualQuant =
437
        computeReflectanceResidual(attrValue, attrPred, quant);
438

439
      double best_score = attrResidualQuant
440
        + kAttrPredLambdaR * (quant.stepSize() >> kFixedPointAttributeShift);
441
442
443
444
445

      for (int i = 0; i < predictor.neighborCount; i++) {
        if (i == aps.max_num_direct_predictors)
          break;

446
447
        attrPred = pointCloud.getReflectance(
          indexesLOD[predictor.neighbors[i].predictorIndex]);
448
        attrResidualQuant =
449
          computeReflectanceResidual(attrValue, attrPred, quant);
450
451

        double idxBits = i + (i == aps.max_num_direct_predictors - 1 ? 1 : 2);
452
        double score = attrResidualQuant
453
454
          + idxBits * kAttrPredLambdaR
            * (quant.stepSize() >> kFixedPointAttributeShift);
455
456
457
458
459
460
461

        if (score < best_score) {
          best_score = score;
          predictor.predMode = i + 1;
          // NB: setting predictor.neighborCount = 1 will cause issues
          // with reconstruction.
        }
462
463
464
465
466
467
468
      }
    }
  }
}

//----------------------------------------------------------------------------

469
void
470
AttributeEncoder::encodeReflectancesPred(
471
  const AttributeDescription& desc,
472
  const AttributeParameterSet& aps,
473
  const std::vector<Quantizers>& quantLayers,
474
475
476
  PCCPointSet3& pointCloud,
  PCCResidualsEncoder& encoder)
{
477
478
  const uint32_t pointCount = pointCloud.getPointCount();
  std::vector<PCCPredictor> predictors;
479
480
  std::vector<uint32_t> numberOfPointsPerLOD;
  std::vector<uint32_t> indexesLOD;
481

482
  buildPredictorsFast(
483
    aps, pointCloud, 0, predictors, numberOfPointsPerLOD, indexesLOD);
484

485
  const int64_t clipMax = (1ll << desc.attr_bitdepth) - 1;
486
  PCCResidualsEntropyEstimator context;
487
488
489
490
491
492
  int zero_cnt = 0;
  std::vector<int> zerorun;
  zerorun.reserve(pointCount);
  std::vector<uint32_t> residual;
  residual.resize(pointCount);

493
  int quantLayer = 0;
494
495
  for (size_t predictorIndex = 0; predictorIndex < pointCount;
       ++predictorIndex) {
496
497
498
499
    if (predictorIndex == numberOfPointsPerLOD[quantLayer]) {
      quantLayer = std::min(int(quantLayers.size()) - 1, quantLayer + 1);
    }
    auto& quant = quantLayers[quantLayer];
500
    auto& predictor = predictors[predictorIndex];
501

502
    computeReflectancePredictionWeights(
503
      aps, pointCloud, indexesLOD, predictorIndex, predictor, encoder, context,
504
      quant[0]);
505

506
    const uint32_t pointIndex = indexesLOD[predictorIndex];
507
    const uint64_t reflectance = pointCloud.getReflectance(pointIndex);
508
    const uint16_t predictedReflectance =
509
      predictor.predictReflectance(pointCloud, indexesLOD);
510
511
    const int64_t quantAttValue = reflectance;
    const int64_t quantPredAttValue = predictedReflectance;
512
513
    const int64_t delta = quant[0].quantize(
      (quantAttValue - quantPredAttValue) << kFixedPointAttributeShift);
514
    const uint32_t attValue0 = uint32_t(IntToUInt(long(delta)));
515
516
    const int64_t reconstructedDelta =
      divExp2RoundHalfUp(quant[0].scale(delta), kFixedPointAttributeShift);
517
518
    const int64_t reconstructedQuantAttValue =
      quantPredAttValue + reconstructedDelta;
519
520
521
    const uint16_t reconstructedReflectance =
      uint16_t(PCCClip(reconstructedQuantAttValue, int64_t(0), clipMax));

522
523
524
525
526
527
528
    if (!attValue0)
      ++zero_cnt;
    else {
      zerorun.push_back(zero_cnt);
      zero_cnt = 0;
    }
    residual[predictorIndex] = attValue0;
529
    pointCloud.setReflectance(pointIndex, reconstructedReflectance);
530
  }
531
532
533
534
535
536
537
538
539

  zerorun.push_back(zero_cnt);
  int run_index = 0;
  encoder.encodeZeroCnt(zerorun[run_index], pointCount);
  zero_cnt = zerorun[run_index++];

  for (size_t predictorIndex = 0; predictorIndex < pointCount;
       ++predictorIndex) {
    auto& predictor = predictors[predictorIndex];
540
    if (predictor.maxDiff >= aps.adaptive_prediction_threshold) {
541
542
543
544
545
546
547
548
549
550
551
552
      encoder.encodePredMode(
        predictor.predMode, aps.max_num_direct_predictors);
    }
    if (zero_cnt > 0)
      zero_cnt--;
    else {
      encoder.encode(residual[predictorIndex]);
      if (predictorIndex != pointCount - 1)
        encoder.encodeZeroCnt(zerorun[run_index], pointCount);
      zero_cnt = zerorun[run_index++];
    }
  }
553
554
555
556
}

//----------------------------------------------------------------------------

557
Vec3<int64_t>
558
AttributeEncoder::computeColorResiduals(
559
560
  const Vec3<uint8_t> color,
  const Vec3<uint8_t> predictedColor,
561
  const Quantizers& quant)
562
{
563
  Vec3<int64_t> residuals;
564
565
  const int64_t quantAttValue = color[0];
  const int64_t quantPredAttValue = predictedColor[0];
566
567
  const int64_t delta = quant[0].quantize(
    (quantAttValue - quantPredAttValue) << kFixedPointAttributeShift);
568
  residuals[0] = IntToUInt(delta);
569
570
571
  for (size_t k = 1; k < 3; ++k) {
    const int64_t quantAttValue = color[k];
    const int64_t quantPredAttValue = predictedColor[k];
572
573
    const int64_t delta = quant[1].quantize(
      (quantAttValue - quantPredAttValue) << kFixedPointAttributeShift);
574
    residuals[k] = IntToUInt(delta);
575
576
577
578
579
580
581
582
  }
  return residuals;
}

//----------------------------------------------------------------------------

void
AttributeEncoder::computeColorPredictionWeights(
583
  const AttributeParameterSet& aps,
584
  const PCCPointSet3& pointCloud,
585
586
  const std::vector<uint32_t>& indexesLOD,
  const uint32_t predictorIndex,
587
588
  PCCPredictor& predictor,
  PCCResidualsEncoder& encoder,
589
  PCCResidualsEntropyEstimator& context,
590
  const Quantizers& quant)
591
592
{
  predictor.computeWeights();
593
  predictor.maxDiff = 0;
594
595
596
  if (predictor.neighborCount > 1) {
    int64_t minValue[3] = {0, 0, 0};
    int64_t maxValue[3] = {0, 0, 0};
597
    for (int i = 0; i < predictor.neighborCount; ++i) {
598
      const Vec3<uint8_t> colorNeighbor =
599
        pointCloud.getColor(indexesLOD[predictor.neighbors[i].predictorIndex]);
600
601
602
603
604
605
606
607
608
609
610
611
      for (size_t k = 0; k < 3; ++k) {
        if (i == 0 || colorNeighbor[k] < minValue[k]) {
          minValue[k] = colorNeighbor[k];
        }
        if (i == 0 || colorNeighbor[k] > maxValue[k]) {
          maxValue[k] = colorNeighbor[k];
        }
      }
    }
    const int64_t maxDiff = (std::max)(
      maxValue[2] - minValue[2],
      (std::max)(maxValue[0] - minValue[0], maxValue[1] - minValue[1]));
612
613
    predictor.maxDiff = maxDiff;

614
    if (maxDiff >= aps.adaptive_prediction_threshold) {
615
616
      Vec3<uint8_t> attrValue =
        pointCloud.getColor(indexesLOD[predictorIndex]);
617
618
619

      // base case: weighted average of n neighbours
      predictor.predMode = 0;
620
      Vec3<uint8_t> attrPred = predictor.predictColor(pointCloud, indexesLOD);
621
      Vec3<int64_t> attrResidualQuant =
622
        computeColorResiduals(attrValue, attrPred, quant);
623
624

      double best_score = attrResidualQuant[0] + attrResidualQuant[1]
625
        + attrResidualQuant[2]
626
627
        + kAttrPredLambdaC
          * (double)(quant[0].stepSize() >> kFixedPointAttributeShift);
628
629
630
631
632

      for (int i = 0; i < predictor.neighborCount; i++) {
        if (i == aps.max_num_direct_predictors)
          break;

633
634
        attrPred = pointCloud.getColor(
          indexesLOD[predictor.neighbors[i].predictorIndex]);
635
        attrResidualQuant = computeColorResiduals(attrValue, attrPred, quant);
636
637
638

        double idxBits = i + (i == aps.max_num_direct_predictors - 1 ? 1 : 2);
        double score = attrResidualQuant[0] + attrResidualQuant[1]
639
          + attrResidualQuant[2]
640
641
          + idxBits * kAttrPredLambdaC
            * (quant[0].stepSize() >> kFixedPointAttributeShift);
642
643
644
645
646
647
648

        if (score < best_score) {
          best_score = score;
          predictor.predMode = i + 1;
          // NB: setting predictor.neighborCount = 1 will cause issues
          // with reconstruction.
        }
649
650
651
652
653
654
655
      }
    }
  }
}

//----------------------------------------------------------------------------

656
void
657
AttributeEncoder::encodeColorsPred(
658
  const AttributeDescription& desc,
659
  const AttributeParameterSet& aps,
660
  const std::vector<Quantizers>& quantLayers,
661
662
663
  PCCPointSet3& pointCloud,
  PCCResidualsEncoder& encoder)
{
664
  const size_t pointCount = pointCloud.getPointCount();
665
  std::vector<PCCPredictor> predictors;
666
667
  std::vector<uint32_t> numberOfPointsPerLOD;
  std::vector<uint32_t> indexesLOD;
668

669
  buildPredictorsFast(
670
    aps, pointCloud, 0, predictors, numberOfPointsPerLOD, indexesLOD);
671

672
  const int64_t clipMax = (1ll << desc.attr_bitdepth) - 1;
673
674
  uint32_t values[3];
  PCCResidualsEntropyEstimator context;
675
676
677
678
679
680
681
  int zero_cnt = 0;
  std::vector<int> zerorun;
  std::vector<uint32_t> residual[3];
  for (int i = 0; i < 3; i++) {
    residual[i].resize(pointCount);
  }

682
  int quantLayer = 0;
683
684
  for (size_t predictorIndex = 0; predictorIndex < pointCount;
       ++predictorIndex) {
685
686
687
688
    if (predictorIndex == numberOfPointsPerLOD[quantLayer]) {
      quantLayer = std::min(int(quantLayers.size()) - 1, quantLayer + 1);
    }
    auto& quant = quantLayers[quantLayer];
689
    auto& predictor = predictors[predictorIndex];
690

691
    computeColorPredictionWeights(
692
      aps, pointCloud, indexesLOD, predictorIndex, predictor, encoder, context,
693
      quant);
694
    const auto pointIndex = indexesLOD[predictorIndex];
695
696
    const Vec3<uint8_t> color = pointCloud.getColor(pointIndex);
    const Vec3<uint8_t> predictedColor =
697
      predictor.predictColor(pointCloud, indexesLOD);
698

699
    Vec3<uint8_t> reconstructedColor;
700
701
702
703
704
705
706
707
708
709
710
711
    for (int k = 0; k < 3; ++k) {
      const auto& q = quant[std::min(k, 1)];
      int64_t residual = color[k] - predictedColor[k];

      int64_t residualQ = q.quantize(residual << kFixedPointAttributeShift);
      int64_t residualR =
        divExp2RoundHalfUp(q.scale(residualQ), kFixedPointAttributeShift);

      values[k] = uint32_t(IntToUInt(long(residualQ)));

      int64_t recon = predictedColor[k] + residualR;
      reconstructedColor[k] = uint8_t(PCCClip(recon, int64_t(0), clipMax));
712
    }
713
    pointCloud.setColor(pointIndex, reconstructedColor);
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733

    if (!values[0] && !values[1] && !values[2]) {
      ++zero_cnt;
    } else {
      zerorun.push_back(zero_cnt);
      zero_cnt = 0;
    }

    for (int i = 0; i < 3; i++) {
      residual[i][predictorIndex] = values[i];
    }
  }

  zerorun.push_back(zero_cnt);
  int run_index = 0;
  encoder.encodeZeroCnt(zerorun[run_index], pointCount);
  zero_cnt = zerorun[run_index++];
  for (size_t predictorIndex = 0; predictorIndex < pointCount;
       ++predictorIndex) {
    auto& predictor = predictors[predictorIndex];
734
    if (predictor.maxDiff >= aps.adaptive_prediction_threshold) {
735
736
737
738
739
740
741
742
743
744
745
746
747
748
      encoder.encodePredMode(
        predictor.predMode, aps.max_num_direct_predictors);
    }
    if (zero_cnt > 0)
      zero_cnt--;
    else {
      for (size_t k = 0; k < 3; k++)
        values[k] = residual[k][predictorIndex];

      encoder.encode(values[0], values[1], values[2]);
      if (predictorIndex != pointCount - 1)
        encoder.encodeZeroCnt(zerorun[run_index], pointCount);
      zero_cnt = zerorun[run_index++];
    }
749
750
751
  }
}

752
753
//----------------------------------------------------------------------------

754
755
void
AttributeEncoder::encodeReflectancesTransformRaht(
756
  const AttributeDescription& desc,
757
  const AttributeParameterSet& aps,
758
  const std::vector<Quantizers>& quantLayers,
759
760
761
  PCCPointSet3& pointCloud,
  PCCResidualsEncoder& encoder)
{
762
763
764
  const int voxelCount = int(pointCloud.getPointCount());
  std::vector<MortonCodeWithIndex> packedVoxel(voxelCount);
  for (int n = 0; n < voxelCount; n++) {
765
    packedVoxel[n].mortonCode = mortonAddr(pointCloud[n], 0);
766
767
768
769
770
    packedVoxel[n].index = n;
  }
  sort(packedVoxel.begin(), packedVoxel.end());

  // Allocate arrays.
771
  int64_t* mortonCode = new int64_t[voxelCount];
772
  const int attribCount = 1;
David Flynn's avatar
David Flynn committed
773
774
  int* attributes = new int[attribCount * voxelCount];
  int* coefficients = new int[attribCount * voxelCount];
775
776
777
778

  // Populate input arrays.
  for (int n = 0; n < voxelCount; n++) {
    mortonCode[n] = packedVoxel[n].mortonCode;
779
780
    const auto reflectance = pointCloud.getReflectance(packedVoxel[n].index);
    attributes[attribCount * n] = reflectance;
781
782
783
784
  }

  // Transform.
  regionAdaptiveHierarchicalTransform(
785
    aps.raht_prediction_enabled_flag, quantLayers, mortonCode, attributes,
786
    attribCount, voxelCount, coefficients);
787

788
  // Entropy encode.
789
  int zero_cnt = 0;
790
  uint32_t value;
791
  for (int n = 0; n < voxelCount; ++n) {
David Flynn's avatar
David Flynn committed
792
    const int64_t detail = IntToUInt(coefficients[n]);
793
    assert(detail < std::numeric_limits<uint32_t>::max());
794
    value = uint32_t(detail);
795
796
797
798
799
800
801
    if (!value)
      ++zero_cnt;
    else {
      encoder.encodeZeroCnt(zero_cnt, voxelCount);
      encoder.encode(value);
      zero_cnt = 0;
    }
802
  }
803
  encoder.encodeZeroCnt(zero_cnt, voxelCount);
804

805
806
  const int64_t maxReflectance = (1 << desc.attr_bitdepth) - 1;
  const int64_t minReflectance = 0;
807
  for (int n = 0; n < voxelCount; n++) {
David Flynn's avatar
David Flynn committed
808
    int64_t val = attributes[attribCount * n];
809
810
811
    const uint16_t reflectance =
      (uint16_t)PCCClip(val, minReflectance, maxReflectance);
    pointCloud.setReflectance(packedVoxel[n].index, reflectance);
812
813
814
815
816
  }

  // De-allocate arrays.
  delete[] mortonCode;
  delete[] attributes;
David Flynn's avatar
David Flynn committed
817
  delete[] coefficients;
818
819
820
821
}

//----------------------------------------------------------------------------

822
823
void
AttributeEncoder::encodeColorsTransformRaht(
824
  const AttributeDescription& desc,
825
  const AttributeParameterSet& aps,
826
  const std::vector<Quantizers>& quantLayers,
827
828
829
  PCCPointSet3& pointCloud,
  PCCResidualsEncoder& encoder)
{
830
831
832
  const int voxelCount = int(pointCloud.getPointCount());
  std::vector<MortonCodeWithIndex> packedVoxel(voxelCount);
  for (int n = 0; n < voxelCount; n++) {
833
    packedVoxel[n].mortonCode = mortonAddr(pointCloud[n], 0);
834
835
836
837
838
    packedVoxel[n].index = n;
  }
  sort(packedVoxel.begin(), packedVoxel.end());

  // Allocate arrays.
839
  int64_t* mortonCode = new int64_t[voxelCount];
840
  const int attribCount = 3;
David Flynn's avatar
David Flynn committed
841
842
  int* attributes = new int[attribCount * voxelCount];
  int* coefficients = new int[attribCount * voxelCount];
843
844
845
846
847
848
849
850
851
852
853
854

  // Populate input arrays.
  for (int n = 0; n < voxelCount; n++) {
    mortonCode[n] = packedVoxel[n].mortonCode;
    const auto color = pointCloud.getColor(packedVoxel[n].index);
    attributes[attribCount * n] = color[0];
    attributes[attribCount * n + 1] = color[1];
    attributes[attribCount * n + 2] = color[2];
  }

  // Transform.
  regionAdaptiveHierarchicalTransform(
855
    aps.raht_prediction_enabled_flag, quantLayers, mortonCode, attributes,
856
    attribCount, voxelCount, coefficients);
857
858

  // Entropy encode.
859
  uint32_t values[attribCount];
860
  int zero_cnt = 0;
861
  for (int n = 0; n < voxelCount; ++n) {
862
    for (int d = 0; d < attribCount; ++d) {
David Flynn's avatar
David Flynn committed
863
      const int64_t detail = IntToUInt(coefficients[voxelCount * d + n]);
864
865
      assert(detail < std::numeric_limits<uint32_t>::max());
      values[d] = uint32_t(detail);
866
    }
867
868
869
870
871
872
873
    if (!values[0] && !values[1] && !values[2])
      ++zero_cnt;
    else {
      encoder.encodeZeroCnt(zero_cnt, voxelCount);
      encoder.encode(values[0], values[1], values[2]);
      zero_cnt = 0;
    }
874
  }
875
876
  encoder.encodeZeroCnt(zero_cnt, voxelCount);

877
  const int clipMax = (1 << desc.attr_bitdepth) - 1;
878
  for (int n = 0; n < voxelCount; n++) {
David Flynn's avatar
David Flynn committed
879
880
881
    const int r = attributes[attribCount * n];
    const int g = attributes[attribCount * n + 1];
    const int b = attributes[attribCount * n + 2];
882
    Vec3<uint8_t> color;
883
884
885
    color[0] = uint8_t(PCCClip(r, 0, clipMax));
    color[1] = uint8_t(PCCClip(g, 0, clipMax));
    color[2] = uint8_t(PCCClip(b, 0, clipMax));
886
887
888
889
890
891
    pointCloud.setColor(packedVoxel[n].index, color);
  }

  // De-allocate arrays.
  delete[] mortonCode;
  delete[] attributes;
David Flynn's avatar
David Flynn committed
892
  delete[] coefficients;
893
}
894
895
896
897
//----------------------------------------------------------------------------

void
AttributeEncoder::encodeColorsLift(
898
  const AttributeDescription& desc,
899
  const AttributeParameterSet& aps,
900
  const std::vector<Quantizers>& quantLayers,
901
902
903
904
  PCCPointSet3& pointCloud,
  PCCResidualsEncoder& encoder)
{
  const size_t pointCount = pointCloud.getPointCount();
905
  std::vector<PCCPredictor> predictors;
906
907
  std::vector<uint32_t> numberOfPointsPerLOD;
  std::vector<uint32_t> indexesLOD;
908

909
  buildPredictorsFast(
910
    aps, pointCloud, 0, predictors, numberOfPointsPerLOD, indexesLOD);
911

912
913
914
915
  for (size_t predictorIndex = 0; predictorIndex < pointCount;
       ++predictorIndex) {
    predictors[predictorIndex].computeWeights();
  }
916
  std::vector<uint64_t> weights;
917
918
919
920
921
922
  if (!aps.scalable_lifting_enabled_flag) {
    PCCComputeQuantizationWeights(predictors, weights);
  } else {
    computeQuantizationWeightsScalable(
      predictors, numberOfPointsPerLOD, pointCount, 0, weights);
  }
923
  const size_t lodCount = numberOfPointsPerLOD.size();
924
  std::vector<Vec3<int64_t>> colors;
925
926
927
928
929
  colors.resize(pointCount);

  for (size_t index = 0; index < pointCount; ++index) {
    const auto& color = pointCloud.getColor(indexesLOD[index]);
    for (size_t d = 0; d < 3; ++d) {
930
      colors[index][d] = int32_t(color[d]) << kFixedPointAttributeShift;
931
932
933
934
935
936
937
938
939
940
941
942
    }
  }

  for (size_t i = 0; (i + 1) < lodCount; ++i) {
    const size_t lodIndex = lodCount - i - 1;
    const size_t startIndex = numberOfPointsPerLOD[lodIndex - 1];
    const size_t endIndex = numberOfPointsPerLOD[lodIndex];
    PCCLiftPredict(predictors, startIndex, endIndex, true, colors);
    PCCLiftUpdate(predictors, weights, startIndex, endIndex, true, colors);
  }

  // compress
943
  int zero_cnt = 0;
944
  int quantLayer = 0;
945
946
  for (size_t predictorIndex = 0; predictorIndex < pointCount;
       ++predictorIndex) {
947
948
949
950
951
    if (predictorIndex == numberOfPointsPerLOD[quantLayer]) {
      quantLayer = std::min(int(quantLayers.size()) - 1, quantLayer + 1);
    }
    auto& quant = quantLayers[quantLayer];

952
    const int64_t quantWeight = weights[predictorIndex];
953
    auto& color = colors[predictorIndex];
954
    const int64_t delta = quant[0].quantize(color[0] * quantWeight);
955
    const int64_t detail = IntToUInt(delta);
956
    assert(detail < std::numeric_limits<uint32_t>::max());
957
    const int64_t reconstructedDelta = quant[0].scale(delta);
958
959
960
961
    color[0] = reconstructedDelta / quantWeight;
    uint32_t values[3];
    values[0] = uint32_t(detail);
    for (size_t d = 1; d < 3; ++d) {
962
      const int64_t delta = quant[1].quantize(color[d] * quantWeight);
963
      const int64_t detail = IntToUInt(delta);
964
      assert(detail < std::numeric_limits<uint32_t>::max());
965
      const int64_t reconstructedDelta = quant[1].scale(delta);
966
967
968
      color[d] = reconstructedDelta / quantWeight;
      values[d] = uint32_t(detail);
    }
969
970
971
972
973
974
975
    if (!values[0] && !values[1] && !values[2])
      ++zero_cnt;
    else {
      encoder.encodeZeroCnt(zero_cnt, pointCount);
      encoder.encode(values[0], values[1], values[2]);
      zero_cnt = 0;
    }
976
  }
977
  encoder.encodeZeroCnt(zero_cnt, pointCount);
978
979
980
981
982
983
984
985

  // reconstruct
  for (size_t lodIndex = 1; lodIndex < lodCount; ++lodIndex) {
    const size_t startIndex = numberOfPointsPerLOD[lodIndex - 1];
    const size_t endIndex = numberOfPointsPerLOD[lodIndex];
    PCCLiftUpdate(predictors, weights, startIndex, endIndex, false, colors);
    PCCLiftPredict(predictors, startIndex, endIndex, false, colors);
  }
986

987
  const int64_t clipMax = (1 << desc.attr_bitdepth) - 1;
988
  for (size_t f = 0; f < pointCount; ++f) {
989
990
    const auto color0 =
      divExp2RoundHalfInf(colors[f], kFixedPointAttributeShift);
991
    Vec3<uint8_t> color;
992
    for (size_t d = 0; d < 3; ++d) {
993
      color[d] = uint8_t(PCCClip(color0[d], 0.0, clipMax));
994
    }
995
    pointCloud.setColor(indexesLOD[f], color);
996
997
998
999
1000
1001
1002
  }
}

//----------------------------------------------------------------------------

void
AttributeEncoder::encodeReflectancesLift(
1003
  const AttributeDescription& desc,
1004
  const AttributeParameterSet& aps,
1005
  const std::vector<Quantizers>& quantLayers,