【点云压缩】Haar小波变换与RAHT自适应区域层级变换

Haar小波

小波变换由一堆小波基和其系数组成，小波基又分为母小波（低频的）和父小波（高频的）。

常用于二维图形处理的小波变换是Haar小波变换，Haar小波变换具有压缩比、抗干扰、速度快的特点，经过小波变换后的系数数据会变得具有规律性，方便后续处理算法进行压缩，同时一些值较小的分量置0不影响图片整体观感。

截取了PCL-AVS-PCC一段小波变换点云压缩的代码

void WaveletCoreTransform(FXPoint* attributes, const int attribCount, const int voxelCount,
                          int* integerizedAttributes) {//attribCount=1
  size_t M, N;
  size_t S, d, i, j;

  N = voxelCount;
  
  FXPoint* attributesTransformed = new FXPoint[voxelCount * attribCount];
  d = 0;  // 层数

  while (N > 1) {
    d++;
    i = 0;  //attributes 当前节点
    M = 0;  //attributesTransformed 当前节点
    S = N;  //当前层节点数

    while (i < S) {
      FXPoint attributes222[3][2] = {};
      for (int k = 0; k < attribCount; k++) {
        attributes222[k][0] = attributes[i * attribCount + k];
        attributes222[k][1] = attributes[(i + 1) * attribCount + k];
      }
      N--;
      i += 2;

      fwdTransform2Node(attribCount, attributes222);

      for (int kk = 0; kk < attribCount; kk++) {
        attributesTransformed[M * attribCount + kk] = attributes222[kk][0];
        attributesTransformed[N * attribCount + kk] = attributes222[kk][1];
      }
      M++;

      if (i == S - 1) {//奇数个节点
        for (size_t k = 0; k < attribCount; k++) {
          attributesTransformed[M * attribCount + k] = attributes[i * attribCount + k];
          attributesTransformed[M * attribCount + k] *= FXPoint(1.41421356);
        }
        i++;
        M++;
      }
    }

    i = S * attribCount;
    j = N * attribCount;
    while (i-- > j)
      attributes[i] = attributesTransformed[i];

    std::swap(attributes, attributesTransformed);
  }

  if (d % 2) {
    attributesTransformed[0] = attributes[0];
    std::swap(attributes, attributesTransformed);
  }
  delete[] attributesTransformed;

  for (i = 0; i < voxelCount; i++)
    for (size_t k = 0; k < attribCount; k++) {
      integerizedAttributes[i + voxelCount * k] = attributes[i * attribCount + k].round();
    }
}

代码逻辑讲解：

attributes数组存放所有属性，attributesTransformed数组与attributes等长，用于存放变换结果。最外层循环进行的条件是当前处理的节点数N>1，内层循环进行的条件是当前处理的节点索引i未大于当前层节点数S，取出待处理的i和i+1处的attributes，进行haar小波变换得到均值和差值的均值，N--，i+=2，haar变换后得到的值，均值放入attributesTransformed[M]，差值的均值放入attributesTransformed[N]，M++（N和M可以看做用来存放结果的游标，i则是读取属性的游标）。如果当前层的节点数是奇数个，直接将最后一个节点的值当做均值存储（假设有五个节点，则存储数组3+2,3个存均值,2个存差值）。

对应的逆变换：

void WaveletCoreInverseTransform(FXPoint* attributes, const int attribCount, const int voxelCount,
                                 int* integerizedAttributes) {
  size_t M, N;
  size_t S, d, i;

  M = N = voxelCount;
  FXPoint* attributesTransformed = new FXPoint[voxelCount * attribCount];
  for (i = 0; i < voxelCount; i++)
    for (size_t k = 0; k < attribCount; k++) {
      attributesTransformed[i * attribCount + k] =
        FXPoint(integerizedAttributes[i + k * voxelCount]);
    }

  size_t S_buffer[32] = {};

  // compute number of nodes for each level
  d = 0;
  while (N > 1) {
    S_buffer[d] = M;
    d++;
    if (M % 2)
      M = M / 2 + 1;
    else
      M = M / 2;
    N = M;
  }

  // Inverse transform
  while (d--) {
    S = S_buffer[d];  //当前层节点数
    M = 0;
    N = S;
    i = 0;

    while (i < S) {
      int flag = i;
      FXPoint attributes222[3][2] = {0};
      N--;
      for (int k = 0; k < attribCount; k++) {
        attributes222[k][0] = attributesTransformed[M * attribCount + k];
        attributes222[k][1] = attributesTransformed[N * attribCount + k];
      }

      invTransform2Node(attribCount, attributes222);

      for (int kk = 0; kk < attribCount; kk++) {
        attributes[i * attribCount + kk] = attributes222[kk][0];
        attributes[(i + 1) * attribCount + kk] = attributes222[kk][1];
      }
      M++;
      i += 2;

      if (i == S - 1) {
        for (size_t k = 0; k < attribCount; k++) {
          attributes[i * attribCount + k] = attributesTransformed[M * attribCount + k];
          attributes[i * attribCount + k] /= FXPoint(1.41421356);
        }
        i++;
        M++;
      }
    }
    i = S * attribCount;
    while (i--)
      attributesTransformed[i] = attributes[i];
  }
  delete[] attributesTransformed;
}

RAHT

Region -adaptive Haar(or hierarchical ) transform 一种区域自适应的哈尔（或层级）变换，目前仍旧是时间复杂度、压缩信噪比的SOTA，于2016年发表于Trans on Image Processing，目前已有三百多次引用。

核心算法区域自适应变换的motivation：

对于三维点云的属性压缩与二维图片的属性压缩不同点在于，三维点云并不紧凑，不能保证每个点都有数据，因此无法直接使用DCT变换进行压缩。

但是本文使用八叉树对每个点进行划分，每个点划定得到的八叉树网格具有最小的Voxel，通过这个最小的Voxel可以把空间点云划分成均匀的voxels，每个voxels存在occupied 和unoccupied两种不同属性。在二维图片的小波变换中，不允许出现空洞和不存在的值的pixels，但是三维空间中存在没有值的Voxels。如何解决这个问题？作者举了一个二维的例子，

对于一个2D平面(a)，并不是所有的区域都存在a，a可以认为是一个点云的属性,所有的a的权重w一开始都是1。图(a)进行水平方向上的变换，两两聚合，对于都存在a的区域，使用变换公式进行计算，并在对应的H（高频信息集合）中填上高频信息；对于只存在一个a的区域，a保持不变；对于都不存在a的区域，保持空。在上图(b)中可以看到，聚合后的b0以及H对应位置所填写的系数c0，权重w等于原来的两个权重w相加。再进行竖直方向聚合，聚合结果如(d)(e)所示。最后，所有区域将会被聚合成样本均值h0，在保存和传输时，只需要记录下h0,i0,g0,g1,e0,e1,e2,c0即可复原出所有的颜色！

只是进行了变换，就能压缩颜色了吗？

由于颜色提取到了变换域，对于较小的变换系数的改变对整体颜色的影响较小，同时变换系数数值间更为接近，提高了压缩效率。压缩是对变换的系数进行量化和编码。

为什么只需要记录这几个系数，不需要记录权重w就能进行逆变换呢，变换矩阵不是w组成的吗？

其实w来自于几何信息，几何信息复原了，生成了八叉树记录下占据网格后，不需要颜色属性也能计算出w来。

上述的例子只是二维平面上的，怎么对应到三维空间中呢？

小波变换具有线性性，可以先沿着x轴这个维度做上述变换，再y，z轴做同样的一维变换即可，或者先做2D变换，再做1D变换也行，顺序改变不会影响变换结果。

参考资料：

小波变换(三): 从实例代码看Haar小波分解和重构 · GitBook

Wavelet Transform

【图像处理笔记】小波变换 - 湾仔码农 - 博客园