Better quaternion quantization #36
Conversation
|
This looks great, thank you for contributing this! Please give me until early next week to look through this a bit more carefully :) |
|
My colleague (@keyboardspecialist) and I took a closer look at this and ran some of our own 3DGS PLY files through this to test things out, given our interest in improved quaternion precision in SPZ. We were interested in two angles:
Precision testingTo test the impact on precision @keyboardspecialist put together a histogram measurement of the angular error across several of our datasets with and without the fix applied. The results between the datasets are identical so we are only sharing one set of data here. The results I'm sharing are taken from running the source data for our substation example through the SPZ library with and without @f-dy's fixes. Overall, the results look good to us. The error is substantially reduced, with the maximum angular error being reduced from 11-degrees to 7-degrees. More importantly, the error is skewed much further towards 0, with the 0.5-degree to 1.0-degree bucket seeing a rise of about 22%. Subsequent buckets have a much sharper drop off compared to the uncorrected version. Charts provided below.
Performance & file size testingI performed performance & file size testing as I had concerns about how this would impact the creation and read processes. Read in particular is extremely important since that is a runtime factor for us. My test cases used mechanisms within the SPZ library to eliminate outside variables. To test creation, I loaded my ply file and created an SPZ. To test read, I loaded an SPZ and created a PLY. I ran each for 10 iterations since we're above the second range for creation even on NOTE: Something on the original test machine I used caused me to get weird results with my testing. A different machine had way better results, resulting in only a 25% uplift: 2.23s to 2.98s with the fix.
Results with
|
|
Thank you for testing! I'm a bit surprised by the 4x time factor for the compression, is it with compile-time optimization turned on? Does it include the zip compression? My timings were for the whole quantization-and-compression process, and I observed a much smaller overhead (10%), but that may be scene-dependent |
|
@f-dy, so am I, and it didn't sit well with me all night. I spent some significant time tweaking things to eliminate variables, including disabling cpu scaling, but kept getting the same result. This morning, I decided to take my test harness and run it on a different machine and got much better results. On this other machine, I'm seeing about a 25% uplift: 2.23s to 2.98s with the fix. Only theory I have right now is that something on my work machine is interfering with this somehow. I've updated my original post. |
|
Thanks for that extended testing! 25% seems very reasonable. One parameter to tune for speed is the maximum radius. I set it to 8 because I got one case (out of millions) where this yielded a better result, but you see from the two examples above (look at "Best radius found") that max radius=3 is probably fine. If targeting real-time compression, I would even limit to radius=2 and hardcode the point offsets in the radius=1 and radius=2 shells. I noticed that, even if a higher radius yields a better result, it is usually only marginally better, i.e. it's almost never bringing an error of 5° down to 0.5°, but it will maybe get something like 4.5°. I didn't figure out what metric could demonstrate that. Do you have an idea, @weegeekps or @keyboardspecialist? |
|
Hi. Thanks for adressing this topic. At Esri, we are interested in guaranteeing a survey-grade accuracy to gaussian splats when using 3D Tiles with a GLTF extension based on SPZ. In this context, a maximum error of 12 degrees or 7 degrees does not work for us. We are proposing to introduce a new version of the SPZ format using the classical "smallest three" quantization of quaternions, on 4 bytes (9 bits is enough for each of the 3 smallest components + 2 bits to remind which component is largest). The PR is here and your feedback is very welcome: #41 The current PR may still be useful to improve the accuracy of SPZ version 2. |
|
Hi @jeanphilippepons ! Hope you're doing well. Smallest three sure makes sense to mitigate that issue, but even with that solution, rounding may not be the best solution (to be tested). |
|
Are you the Frédéric Devernay I've known at INRIA 20 years ago? |
|
Does this become obsolete with #41 ? |
Yes
https://marc-b-reynolds.github.io/quaternions/2017/05/02/QuatQuantPart1.html |
|
Is it not better to encode the orientation as spherical coordinates + a roll angle? This will give an upper-bound on the rotation error. Here’s the Monte Carlo histogram for the (φ, θ, ψ) = (11, 10, 11) bit allocation for a total budget of 32bit. About 500 000 random orientations were quantized and compared: 
import numpy as np
import matplotlib.pyplot as plt
def quat_from_axis_angle(axis, angle):
axis = axis / np.linalg.norm(axis, axis=-1, keepdims=True)
half = angle / 2
s = np.sin(half)
return np.stack([np.cos(half), *(axis.T * s)], axis=-1)
def quat_mult(q1, q2):
w1, x1, y1, z1 = q1.T
w2, x2, y2, z2 = q2.T
return np.stack([
w1*w2 - x1*x2 - y1*y2 - z1*z2,
w1*x2 + x1*w2 + y1*z2 - z1*y2,
w1*y2 - x1*z2 + y1*w2 + z1*x2,
w1*z2 + x1*y2 - y1*x2 + z1*w2
], axis=-1)
def quat_conj(q):
q = np.array(q)
q[...,1:] *= -1
return q
def quat_angle_error(q1, q2):
dq = quat_mult(q1, quat_conj(q2))
dq /= np.linalg.norm(dq, axis=-1, keepdims=True)
ang = 2 * np.arccos(np.clip(np.abs(dq[...,0]), -1, 1))
return np.degrees(ang)
# Bit allocation
b_phi, b_theta, b_psi = 11, 10, 11
# Quantization steps
step_phi = 2*np.pi / (2**b_phi)
step_theta = np.pi / (2**b_theta)
step_psi = 2*np.pi / (2**b_psi)
# Monte Carlo sampling
N = 500_000
phi = np.random.rand(N) * 2*np.pi
theta = np.arccos(2*np.random.rand(N) - 1)
psi = np.random.rand(N) * 2*np.pi
# Quantize angles
phi_q = (np.floor(phi / step_phi + 0.5) * step_phi) % (2*np.pi)
theta_q = np.clip(np.floor(theta / step_theta + 0.5) * step_theta, 0, np.pi)
psi_q = (np.floor(psi / step_psi + 0.5) * step_psi) % (2*np.pi)
# Convert to quaternions
axes = np.stack([np.sin(theta)*np.cos(phi),
np.sin(theta)*np.sin(phi),
np.cos(theta)], axis=-1)
axes_q = np.stack([np.sin(theta_q)*np.cos(phi_q),
np.sin(theta_q)*np.sin(phi_q),
np.cos(theta_q)], axis=-1)
# Build quaternions
q_true = quat_from_axis_angle(axes, psi)
q_quant = quat_from_axis_angle(axes_q, psi_q)
# Compute angular errors
errors_deg = quat_angle_error(q_true, q_quant)
# Statistics
percentiles = np.percentile(errors_deg, [50, 90, 99, 99.9])
worst = np.max(errors_deg)
# Plot histogram
plt.figure(figsize=(6,4))
plt.hist(errors_deg, bins=300, range=(0,0.5), density=True)
plt.xlabel("Rotation error (degrees)")
plt.ylabel("Probability density")
plt.title(
f"Monte Carlo orientation quantization errors (11,10,11 bits)\n"
f"median={percentiles[0]:.3f}°, 90%={percentiles[1]:.3f}°, "
f"99%={percentiles[2]:.3f}°, max≈{worst:.3f}°"
)
plt.grid(True)
plt.show() |
|
This PR could/should be closed. Alternative suggestions for encoding could/should be discussed in dedicated issues. |
Fixes #31
As explained in #31, the method used by spz for quaternion quantization may cause very large errors. An example is given in that issue that has a 8° rotation error, and while the effect of such an error may be barely visible for small Gaussians, it may clearly stand out as artifacts if the model has large elongated Gaussians.
This PR fixes it by keeping backward compatibility with the SPZ format, but adopting a slightly better quantization strategy, which works by gradually exploring "shells" of increasing radius in the xyz 8-bit quantized quaternion space, centered around the default value (which is just the result of rounding xyz).
The code outputs statistics about rotation quantization. For example, in this case:
29.93% of the points had a reconstruction error less than 0.5° with the current (rounding) strategy.
51.13% of the points had a reconstruction error more than 0.5° with the current strategy, but less than 1.0° when exploring the first shell as well (3x3x3 neighborhood)
etc. (those thresholds were set experimentally and they are in the code)
for 63.76% of the points, the best reconstruction error was obtained using the current strategy
for 32.56% of the points, the best reconstruction error was obtained from a value in the first shell
etc.
In terms of execution speed, the resulting code is ~10% slower than the original version
It has been found experimentally that going beyond r=8 never brougth a better solution (in the example above, only 1 point out of 116380 has the best solution for r=8).
Although one could argue that the reduction of the average error is small, but most of the error is made of a few large outliers that really stand out visually if they correspond to large Gaussians.
Here are the results on the lego ply produced by the original code:
As you can see, it's highly scene dependent (the lego scene is perfectly axis-aligned, and my scene above is slightly misaligned, which is the worst case for the original quaternion quantization method), but the good news is that:
This PR also includes a small standalone utility to do the compression/decompression (inspired by another fork), and that utility can also run the original compression method, as well as set the antialiasing/mip-splatting flag in the SPZ.
It also improves PLY parsing but skipping comment lines (as those produced by nerfstudio's exporter).