particle filling

[Mr-Ma-yikun]

I encountered excessively long runtime (no result after two days) while using particle filling.
May I ask if there are any parameters that need adjustment or improvement?

  particle_filling: {
        "n_grid": 128,
        "density_threshold": 40.0,
        "search_threshold": 0.5,
        "search_exclude_direction": 5,
        "ray_cast_direction": 0,
        "max_particles_num": 2000000,
        "max_particles_per_cell": 4,
        "boundary": [
            0.2,
            0.8,
            0.2,
            0.8,
            0.1,
            1.0
        ],
        "smooth": 
    }

[r13944003]

Hi! I also encountered an issue where fill_particles() appeared to hang indefinitely. It seemed like the problem was inside densify_grids() kernel, where an unstable covariance matrix can cause the radius r to explode.
Inside densify_grids(), the radius of influence r is computed from the eigenvalues of the covariance matrix. If any eigenvalue is extremely small (or numerically unstable), the computed r can become huge, like this:

[WARNING] Large r at particle 12119 -> 2147483647.000000

The fix is to add a safety upper bound to r. Here's the fixed version of densify_grids():

@ti.kernel
def densify_grids(
    init_particles: ti.template(),
    opacity: ti.template(),
    cov_upper: ti.template(),
    grid: ti.template(),
    grid_density: ti.template(),
    grid_dx: float,
):
    for pi in range(init_particles.shape[0]):
        pos = init_particles[pi]
        x, y, z = pos[0], pos[1], pos[2]
        i = ti.floor(x / grid_dx, dtype=int)
        j = ti.floor(y / grid_dx, dtype=int)
        k = ti.floor(z / grid_dx, dtype=int)

        if 0 <= i < grid.shape[0] and 0 <= j < grid.shape[1] and 0 <= k < grid.shape[2]:
            ti.atomic_add(grid[i, j, k], 1)

        # Construct full covariance
        cov = ti.Matrix([
            [cov_upper[pi][0], cov_upper[pi][1], cov_upper[pi][2]],
            [cov_upper[pi][1], cov_upper[pi][3], cov_upper[pi][4]],
            [cov_upper[pi][2], cov_upper[pi][4], cov_upper[pi][5]],
        ])

        # Check for NaN or inf
        valid = True
        for p in ti.static(range(3)):
            for q in ti.static(range(3)):
                if ti.math.isnan(cov[p, q]) or ti.math.isinf(cov[p, q]):
                    valid = False
        if not valid:
            print(f'[ERROR] Particle {pi} has invalid value in cov')
            continue  # skip this particle

        # Eigen decomposition
        sig, Q = ti.sym_eig(cov)
        for idx in ti.static(range(3)):
            sig[idx] = ti.max(sig[idx], 1e-8)

        sig_mat = ti.Matrix(
            [[1.0 / sig[0], 0, 0], [0, 1.0 / sig[1], 0], [0, 0, 1.0 / sig[2]]]
        )
        cov = Q @ sig_mat @ Q.transpose()

        # Compute r
        r = 0.0
        for idx in ti.static(range(3)):
            s = ti.sqrt(sig[idx])
            sig[idx] = s
            r = ti.max(r, s)
        r = ti.ceil(r / grid_dx, dtype=int)
        
        max_r = 6  # Safety Upper Bound For r
        if r > max_r:
            print("[WARNING] Large r at particle", pi, "->", r)
            r = max_r

        for dx in range(-r, r + 1):
            for dy in range(-r, r + 1):
                for dz in range(-r, r + 1):
                    gi, gj, gk = i + dx, j + dy, k + dz
                    if (
                        0 <= gi < grid_density.shape[0] and
                        0 <= gj < grid_density.shape[1] and
                        0 <= gk < grid_density.shape[2]
                    ):
                        density = compute_density(
                            ti.Vector([gi, gj, gk]), pos, opacity[pi], cov, grid_dx
                        )
                        ti.atomic_add(grid_density[gi, gj, gk], density)

Hope it helps!