65.9K
CodeProject 正在变化。 阅读更多。
Home

使用球面坐标的 16 位普通向量压缩

johnwhile
starIconstarIconstarIconstarIconstarIcon

5.00/5 (13投票s)

2023年8月22日

CPOL

3分钟阅读

viewsIcon

16866

downloadIcon

265

一种使用 C# 将 3D 单位向量编码和解码为 16 位值的简单方法

引言

在使用 Directx 库编写 C# 渲染工具代码时，我需要减小法向量的大小。我意识到，随着顶点数量的增加，发送到 GPU 的数据和/或简单网格几何的文件存储的大小会过度增加。

法线通常用于 3D 图形中以确定表面的方向，它由一个具有单位长度的 3xfloat 向量组成，因此 x y z 的值在 [-1.0, 1.0] 范围内。

通过 Google 搜索，我找到了几种方法，但对我来说最有效的方法是本文中描述的向量的球形分布

我会把文章留给你们阅读，以获取更多信息。

Directx 坐标系中的极坐标

由于目的是使用 Directx 库，因此有必要为左手坐标系定义一些基本的数学运算。

public struct Vector3f
{
   public float x, y, z;
   public float Normalize()
   {
      var length = (float)Math.Sqrt(x * x + y * y + z * z);
      var l = length > float.Epsilon ? 1f / length : 0;
      x *= l;
      y *= l;
      z *= l;
      return length;
   }
}

将极坐标向量转换为单位向量，反之亦然。在笛卡尔坐标系到球坐标系的转换中，我们必须考虑到我们只使用法向量，并且必须确保没有被零除的情况。

public static Vector3f SphericalToCartesian(float r, float theta, float phi)
{
   return new Vector3f(
     r * Math.Cos(theta) * Math.Sin(phi),
     r * Math.Cos(phi),
     r * Math.Sin(theta) * Math.Sin(phi));
}

public static (float r, float theta, float phi) CartesianToSpherical(Vector3f cartesian)
{
   //now x, y, z are divided by r
   float r = cartesian.Normalize();
            
   //work for vector -z
   float theta = (float)Math.Atan2(cartesian.z, cartesian.x); 
   //float theta = (float)Math.Atan(cartesian.z /cartesian.x);

   float phi = (float)Math.Acos(cartesian.y);
   return (r, theta, phi);
}

单位向量量化

为了降低复杂性，我决定将问题简化为仅考虑正笛卡尔坐标，并将法线的符号存储在前 3 位（我们在小端模式下工作）。

ushort value = 0;
if (normal.x < 0) { value |= 1 << 15; normal.x *= -1; }
if (normal.y < 0) { value |= 1 << 14; normal.y *= -1; }
if (normal.z < 0) { value |= 1 << 13; normal.z *= -1; }

这样，phi 和 theta 角位于 [0, π/2] 范围内。要量化角度，您可以简单地选择您喜欢的细分，例如，对于 N 个细分，我们将得到

d_phi = π/2 / N
d_theta = π/2 / N

因此，角度可以用索引 i 和 j 表示，其中

phi = d_phi * i 
theta = d_theta * j

其中 i 和 j 在 [0,N] 范围内

但是，如文章中所写，我们将获得朝向 Y 极的高密度。为了提高量化向量分布的均匀性，我选择根据 phi 改变 theta 角。考虑到对于 phi=0，我们有被零除的情况，但任何 theta 角都是允许的，我们将写成

float phi = i * d_phi;
float theta = i > 0 ? (j/i * π/2) : 0;

现在的问题是将这两个索引 i 和 j 存储在一个 13 位 数字中，因为 3 位 保留给符号。

表 [i,j]

如果我们构造一个表格，其中行是 i, 列是 j，并且该序列的编号是渐进的

int n = 0;
for (byte i = 0; i <= N; i++)
    for (int j = 0; j <= i; j++)
        i_tab[n++] = i;

我们有一个简单的序列，另请参阅 https://en.wikipedia.org/wiki/Triangular_number。

使用公式 n = (i+3)*i/2 (或 n=(i+1)*i/2+i )，我们可以计算用于量化法向量的四分之一球体上的最大点数。

对于 i 和 j，使用 N=126 细分，该表生成 8128 个点，幸运的是，使用 13 位数字，我们可以对从 0 到 8191 的数字进行编码。

编码

static int sum(int i) => (i + 1) * i / 2;

public static ushort Encode(Vector3f normal)
{
   ushort value = 0;
   if (normal.x < 0) { value |= 1 << 15; normal.x *= -1; }
   if (normal.y < 0) { value |= 1 << 14; normal.y *= -1; }
   if (normal.z < 0) { value |= 1 << 13; normal.z *= -1; }

   Mathelp.CartesianToSpherical(normal, out _, out float theta, out float phi);

   //need round due to C# division approximation, 
   //example: it can return 7.99999 instead 8
   int i = (int)Math.Round(phi / d_phi); 
   int j = (int)Math.Round(theta * i * 2 / PI);

   int n = sum(i) + j;
#if DEBUG
   if (n > MAX || n < 0) throw new Exception("Fail");
#endif       
   value |= (ushort)n;

   return value;
}

解码

要反向计算，我们首先需要从 13 位值中获取 i 和 j 索引。我们可以使用两种方法

  1. 使用预先计算的表，但需要存储一个 byte[8128] 表
    int i = i_tab[value];
int j = value - sum(i);
  2. 使用 inverse 函数。例如，这些公式可以在 HLSL Shader 代码中使用。目前，sqrt 函数根本不影响性能（与较旧的图形架构相比）。但是，必须检查每个值，因为 float 的精度可能会给出不同的值。
    static void inverse(int n, out int i, out int j)
{
   //check for n[0,8127]
   i = (int)(Math.Sqrt(1 + 8 * n) - 1) / 2;
   j = n - sum(i);
}
/// <summary>
/// approximation of <see cref="inverse"/> (with sqrt)
/// </summary>
static void inverse_aprox1(int n, out int i, out int j)
{
   //check for n[0,8127]
   i = (int)(Math.Sqrt(n) * rad2);
   j = n - sum(i);
   if (j < 0) { j += i; i--; }
}
/// <summary>
/// approximation of <see cref="inverse"/> (without sqrt)
/// </summary>
static void inverse_aprox2(int n, out int i, out int j)
{
   //check for n[0,8127]
   i = (int)(Math.Exp(0.5 * Math.Log(n)) * rad2);
   j = n - sum(i);
   if (j < 0) { j += i; i--; }
}

HLSL shader 代码

已使用 directx11 测试，但我现在不知道如何衡量性能

float3 DecodeUnitVector16(min16uint encode)
{
    int n = encode & 0x1FFF;
    int i = (sqrt(1 + 8 * n) - 1) / 2;
    int j = n - (i + 1) * i / 2;
    
    float phi = i * 1.5707963267 / 126;
    float theta = i > 0 ? j * 1.5707963267 / i : 0;

    float3 normal = float3(cos(theta) * sin(phi),cos(phi),sin(theta) * sin(phi));
    
    if ((encode & 0x8000) != 0) normal.x *= -1;
    if ((encode & 0x4000) != 0) normal.y *= -1;
    if ((encode & 0x2000) != 0) normal.z *= -1;
    
    return normal;
}

结果

如果我们使用此代码创建所有可能的值，我们将获得密集且均匀的单位向量分布，总共有 65.024 个可能的点，如图所示。

for (int sign = 0; sign < 8; sign++)
    for (int n = 0; n <= 8127; n++)
    {
        int code = n | sign << 13;
        Vertices.Add(Decode((ushort)code));
    }

欢迎提出任何意见或建议以改进代码。

扩展到 24 位

在源代码中，我添加了 24 位扩展，其中 N=2046 为四分之一球体生成 2.096.128 个点，因此总共有 16.769.024 个可能的法线。

扩展到 32 位没有意义，也因为我遇到内存溢出错误。

性能

瓶颈是由于三角函数造成的。通过用近似值替换标准函数可以获得两倍的速度（正在测试中）

/// <summary>
/// Bhāskara I's sin
/// </summary>
static float Sin(float r)
{
    float n = (Mathelp.PI - r)*r * 4;
    return n * 4 / (5 * Mathelp.PI * Mathelp.PI - n);
}
/// <summary>
/// Bhāskara I's cos
/// </summary>
static float Cos(float r)
{
    return (Mathelp.PI2 - 4 * r * r) / (Mathelp.PI2 + r * r);
}
/// <summary>
/// Trey Reynolds' ArcCos
/// </summary>
static float AcosTR(float x)
{
    return (float)(8 / 3.0 * Math.Sqrt(2 - Math.Sqrt(2 + 2 * x)) - 1 / 3.0 * 
                             Math.Sqrt(2 - 2 * x));
}
/// <summary>
/// Sebastien Lagarde's ArcCos
/// </summary>
static float AcosSL(float x)
{
    var z = Math.Abs(x);
    z = (float)((-0.168577f * z + 1.56723f) * Math.Sqrt(1 - z));
    return x < 0 ? z + Mathelp.PI : -z;
}

历史

  • 2023 年 8 月 22 日：初始版本
© . All rights reserved.