Particles! (2/?)

Español


Hello!

Another article in the particles series! The first one presented all the results, so now I want to talk about different parts of the final product. In this case I want to talk about the colorization, or how the particles end up being red, white, blue or some other mix of those.

Red, White and Density

Like I said in the first entry, I wanted to show the properties of each type of random position generation using a density visualization, in essence that means that the color is an indicator of how populated is the space around each particle: redder particles have less neighbors, white particles have more. The detail is “less” and “more” are defined relatively to the maximum value in the entire amount of particles. The value of density then is very simple (conceptually) to calculate but a little costly.

Calculating Density

The most basic version of the density calculation is as simple as:

void update_density(std::vector<Particle>& particles) {
    for (auto& p1 : particles) {
        for (auto& p2 : particles) {
            if (p1.dist2(p2) < THRESHOLD_DIST) {
                p1.density++;
            }
        }
    }
}

The two most important things to notice here are,

  • There is a lot of wasted computing repeating operations.
  • This function order of execution is O(N^2) with N being the number of particles

Don’t repeat work

The first approach to improving this 2 problems is very simple and involves avoid doing the same calculation twice. On each loop the same particle will take the place of p1 and then p2 so we’ll eventually do p1.dist2(p2) & then p2.dist2(p1). Let’s avoid this:

void update_density(std::vector<Particle>& particles) {
    for (auto ix = 0u; ix < particles.size(); ix++) {
        auto& p1 = particles[ix];
        for (auto jx = ix + 1; jx < particles.size(); jx++) {
            auto& p2 = particles[jx];
            if (p1.dist2(p2) < THRESHOLD_DIST) {
                p1.density++;
                p2.density++;
            }
        }
    }
}

With this small change, we avoid doing lots of calculations twice, because each iteration updates the density of both particles (and we prevent from calculating again both values by only comparing each particles with the following ones).

The order of execution did not changed that much though, we now do comparisons in this series:


The order still is quadratic. We need a way to drastically reduce the amount of calculations we do.

Locality of the problem

We’re trying to visualize the particle density value at each particle which means we need a threshold value small enough to have some sense in the final image. Since we’re plotting particles in the surface of the unit sphere (1 unit radius), the maximum distance between 2 particles is 2 units (diameter). I picked sqrt(0.004) ~ 0.063 (the tests are done with all distances squared so that we don’t need to use the square root) as my threshold value, it could’ve been any number but this one resulted we from the beginning.

If you compare both values you know that many, MANY particles are far away from each other, we can use a data structure suitable to just discard as many comparisons as possible before even doing them. Structures for this kind of range based distance calculations are within the family of the spatial tree structures like quadtrees, octrees, BSP trees and kd-trees to name a few.
In fact, in ray tracing using the Photon Mapping technique, it’s very common to use kd-trees as the Photon Map base for the quality of the distance range queries to retrieve photons closer to a given position in space.

The only caveat it’s that it might get costly to maintain a kd-tree if the points are dynamically created an removed (the tree will eventually have to be re balanced many times), since I wanted to eventually animate the particles, that was an inconvenience.

What I ended up implementing is a hacky-dirty version of a fixed size leaf octree (anyone can argue that this code could be licensed under CRAPL instead of MIT, but hey, it worked fine and was easy ;P). Instead of dynamically creating a tree structure, I defined 2 parameters, a bounding box for the model (given by the maximum and minimum values any coordinate can take) and a number of intervals each axis will be subdivided.

The base structure is a hashmap of vectors (sets might’ve worked too, even better…) where each key is the indentifier generated based on a position and the value is a vector of all the elements that are within one leaf of this thing.

Now, instead of comparing with all the particles, a great portion of them is discarded with little trouble and only the best candidates are used in the comparison.

SPP or “Space Partitioned Positions”

This is the accelerator structure I’ve used:

class spp {
public:
    spp(uint8_t intervals_per_axis, float min_val, float max_val);
    ~spp() = default;

    uint32_t add(const glm::vec3& pos, size_t external_idx);
    void remove(const glm::vec3& pos, size_t external_idx);
    void remove(uint32_t bucket_id, size_t external_idx);
    std::vector<uint32_t> get_buckets_area(const glm::vec3& pos) const;
    std::vector<uint32_t> get_buckets_area(uint32_t bucket_id) const;
    const std::vector<size_t>& get_bucket(uint32_t bucket_id) const;
private:
    std::unordered_map<uint32_t, std::vector<size_t>> m_buckets;
    uint8_t m_intervals_per_axis;
    glm::vec3 m_min_vec;
    float m_normalize_value;

    uint32_t get_bucket(const glm::vec3& pos) const;
};

Step by step:

spp::spp(const uint8_t intervals_per_axis, const float min_val, const float max_val)
    : m_intervals_per_axis(intervals_per_axis)
    , m_min_vec{ min_val, min_val, min_val }
    , m_normalize_value(max_val - min_val) {
    SPL_ASSERT(m_normalize_value != 0., "SPP interval can't be zero.");
    SPL_ASSERT(max_val > min_val, "Consider swapping min and max!");
}

This is the constructor, initializes the auxiliary members and makes some simple controls. The members are a vector (to simplify sintax, it could’ve been just the value) for the minimum value, another for the interval count (with a maximum number of intervals equal to 0xFF or 255) and a third one for the length of the bounding box (per axis). (By how it works, using either 0 or 1 as the number of intervals will have the same effect).

uint32_t spp::add(const glm::vec3& pos, const size_t external_idx) {
    const auto bid = get_bucket(pos);
    m_buckets[bid].push_back(external_idx);
    return bid;
}

This method adds an element to the structure, the element is referenced with the external_idx (in my case I was using indices to reference each particle outside) and the pos variable is the position of the element.

As you can see, it generates the key associated with the storage bucket using the position, stores the id in the bucket and finally returns the bucket identifier. Since calculating the id from the position is a costly operation, in the rest of the code I optimized a bit the usage of the structure by storing the id of the bucket each particle belongs to. This is not necessary since the position alone is enough to map a particle to a bucket.

void spp::remove(const glm::vec3& pos, const size_t external_idx) {
    remove(get_bucket(pos), external_idx);
}
void spp::remove(const uint32_t bucket_id, const size_t external_idx) {
    auto& b = m_buckets[bucket_id];
    b.erase(std::find(b.begin(), b.end(), external_idx));
}

Now we have the remove methods, one overload using the position and the other using directly the bucket id, this is the operation that could be used to argue in favor of using a hashmap of sets instead of one of vectors, but I didn’t do any performance comparison to make the change (others parts of the code might work worse in a set than in a vector, because of cache misses and data locality).

Now let’s jump to the interesting method, get_bucket:

uint32_t spp::get_bucket(const glm::vec3& pos) const {
    const auto norm = (pos - m_min_vec) / m_normalize_value;
    return pack(
        static_cast<uint8_t>(std::floor(m_intervals_per_axis * norm.x)),
        static_cast<uint8_t>(std::floor(m_intervals_per_axis * norm.y)),
        static_cast<uint8_t>(std::floor(m_intervals_per_axis * norm.z)));
}

You can see here what’s being done, each coordinate is mapped to the [0, 1] interval using the constructor parameters and then that value is used to obtain a uint8_t value based on the number of intervals. Finally the 3 values are packed on a single uint32_t value using the pack method.

How does this works? Let’s think about it in a single dimension:

Using 8 as the number of intervals (intervals_per_axis), suppose we get pos = -0.3:


So our bucket id is number 2.
Doing the same for 0.5 we get: bucket id = 6. And finally, using a value very close to 0.5 we should get either the same bucket or the contiguous ones, for 0.55 -> 6, for 0.45 -> 5, for 0.75 -> 7. We are effectively dividing the interval in smaller ones with the lenght equal to the size of the bounding box (per axis) divided by the number of intervals and assigning the ids based on within which interval the position falls.

If you’ve ever read about scalar quantization or float packing (like color channels into a single integer) this is the same.

The general result is:

  • It is certain that the maximum distance within elements on the same bucket is equal to:
  • Unfortunately that doesn’t mean a single bucket is enough to calculate the density around each particle, since 2 very close particles could end up in neighboring buckets but…
  • We can obtain all the particles that are within X value by using the bucket the particle falls and all the neighboring ones provided that X is smaller than m_normalize_value/m_intervals_per_axis.

In our one dimensional case, the maximum distance between 2 values in the same bucket is 0.25, and taking any position P within bucket B, we can calculate all the positions with a distance smaller than 0.25 to P by analyzing only the particles in buckets B, B – 1 (if B != 0) and B + 1 (if B != 7). So by dividing our original set in 8, we’ve reduced the number of positions we need to test to those only in 3 of the 8 buckets (in the worst case).

Moving on:

std::vector<uint32_t> spp::get_buckets_area(const uint32_t bucket_id) const {
    std::vector<uint32_t> buckets;
    buckets.reserve(3 * 3 * 3); // Max adjacent buckets (think of a rubik's cube)
    const auto unpacked = unpack(bucket_id);
    for (int8_t ix = -1; ix < 2; ix++) {
        const int8_t xval = unpacked[0] + ix;
        if (xval >= 0 && xval < m_intervals_per_axis) {
            for (int8_t iy = -1; iy < 2; iy++) {
                const int8_t yval = unpacked[1] + iy;
                if (yval >= 0 && yval < m_intervals_per_axis) {
                    for (int8_t iz = -1; iz < 2; iz++) {
                        const int8_t zval = unpacked[2] + iz;
                        if (zval >= 0 && zval < m_intervals_per_axis) {
                            const auto potential_id = pack(xval, yval, zval);
                            if (m_buckets.find(potential_id) != m_buckets.end()) {
                                buckets.push_back(potential_id);
                            }
                        }
                    }
                }
            }
        }
    }
    return buckets;
}

This code allows us detect this set of bucket+neighbors by calculating the valid values for each part of the packed id (imagine a rubik’s cube center, and the positions of all the other pieces). Given a single bucket id we get that one and the 26 neighboring ones (at worst). Finally, pack and unpack are methods to pack 3 (they could be 4) 8 bit integers into a single 32 bit one.

static uint32_t pack(const uint8_t x, const uint8_t y, const uint8_t z) {
    return  x << 16 | y << 8 | z;
}
static std::array<uint8_t, 3> unpack(uint32_t value) {
    return{
        (value >> 16) & 0xFF,
        (value >> 8) & 0xFF,
        value & 0xFF
    };
}

Using the SPP

Now that we have our SPP structure up and running, instead of computing the density by calculating the distance to all particles, we:

  • Add all particles to the SPP and then,
  • For each particle updated we obtain the set of buckets neighboring it and,
  • Calculate the distance to all particles within those buckets and those only:
void load_particles(std::vector<Particle>& particles) {
    for (auto ix = 0u; ix < particles.size(); ix++) {
        // Note we store the bucket id outside to reduce the amount of times we need to calculate it later.
        p.bucket_id = spp.add(p.position, ix);
    }
}

void update_density(std::vector<Particle>& particles) {
    for (auto ix = 0u; ix < particles.size(); ix++) {
        auto& p1 = particles[ix];
        auto buckets = spp.get_buckets_area(p1.bucket_id);
        for (auto jx = 0u; jx < buckets.size(); jx++) {
            auto& bucket = spp.get_bucket(buckets[jx]);
            for (auto kx = 0u; kx < bucket.size(); kx++) {
                auto& p2 = particles[bucket[kx]];
                if (p1.dist2(p2) < THRESHOLD_DIST) {
                    p1.density++;
                }
            }
        }
    }
}

Note that now, I’m again repeating a bit of work, since I’m not updating both density values in the same loop, the thing is I can no longer ensure how the particles will be iterated, so trying to avoid counting twice gets more complicated. The savings of using the spp are more than enough to compensate for this (and the possible complexity of an algorithm that would correctly skip already made comparisons).

After this, calculating the color for each particle it’s just a matter of finding the maximum density and getting an RGBA value based on the density of each particle in comparison with the maximum one. I’ve used someting along this lines:

Color color_for_particle(Particle& p, const uint32_t max_density) {
    float alive_value = p.isAlive() ? 1.f : 0.f; //Just cast it if you like 🙂
    return Color {
        alive_value,
        (p.density / max_density) * alive_value,
        (p.density / max_density) * alive_value,
        alive_value
    };
}

You can see in the history of the repo, that at first this calculation was performed on CPU and then forwarding the color directly on the shader. Obviously this is not the best thing, in the current version, only the maximum value is calculated in CPU then passed as a uniform, and the density por particle is passed as another attribute. That way, the calculation can be performed on the GPU.

Blue, Red and distance to center

In this demo, the calculations are much simpler, nothing is done on the CPU and the shader outputs Red if the particle distance to the center is less than 1 and Blue otherwise.

Color color_for_particle_blue(Particle& p) {
    float alive_value = p.isAlive() ? 1.f : 0.f; //Just cast it if you like 🙂
    float r = (p.length <= 1.) ? 1.f : 0.f; //Just cast it if you like 🙂
    return Color {
        r,
        0.f,
        1.f - r,
        alive_value
    };
}

Footnote on the SPP number of intervals

The number of buckets directly influences the amount of comparisons we need to perform but also the maximum distance 2 particles can have. We need to instantiate the SPP with the correct value so to skip as many particles as possible and to avoid missing particles that should be part of the density. To calculate it we need to ‘enclose’ the volume of closeness around the particle in the cube created by the bucket the particle belongs and all the neighbors:

  • Having chosen a density threshold equal to srqr(0.004) ~ 0.0632 we need an bucket size (per axis) that is at least this big.
  • We know the maximum value a coordinate can have is 1, and the minimum is -1, so our bbox is 2 units wide on each axis





This number of intervals per axis, give us a bucket size (per axis) of ~0.645 units, the smallest value not smaller than our threshold. This ensures that even the worst case scenario is covered, and all particles that will contribute to the density of a given particle are tested.

The worst case scenario is given by a particle laying on any of the bordering planes of the bucket. Thinking again in our 1 dimensional sample, if our threshold is 0.26, having 8 intervals is too much: a sample with the position 0.5 will be in bucket 6, and a sample with position 0.245 will end up in bucket 4. Those 2 won’t be compared since are not in neighboring buckets even if the distance between them is less than 0.26.

OK! The end! You’ve made it through a long article about details of an unimportant demo, congrats! At least I hope you get something useful out of it. You can refer to the previous article for the original story and the source code.
Back to Top




¡Buenas!

¡Otro artículo en la serie de las partículas! En el primero presenté los resultados, así que ahora quiero hablar de las diferentes partes del producto final. En este caso, sobre cómo colorear las partículas para que terminen siendo rojas, blancas, azules o alguna combinación.

Rojo, Blanco y la Densidad

Como comenté en la primera entrada, quería mostrar las propiedades de cada función para generar posiciones aleatorias usando una visualización de densidad, en esencia significa que el color es un indicador de qué tan poblado es el espacio alrededor de cada partícula: mientras más rojo, menos vecinos, a más blanco más. El detalle es que “más” y “menos” están definidos de manera relativa al valor máximo de densidad entre todas las partículas. Este valor de densidad es bastante simple de calcular (conceptualmente) pero un poco costoso.

Calculando la densidad

La forma más básica de calcular esta densidad es tan simple como:

void update_density(std::vector<Particle>& particles) {
    for (auto& p1 : particles) {
        for (auto& p2 : particles) {
            if (p1.dist2(p2) < THRESHOLD_DIST) {
                p1.density++;
            }
        }
    }
}

Las dos cosas más importantes a notar acá son,

  • Hay un una gran cantidad de poder de cómputo perdido repitiendo operaciones.
  • El órden de ejecución de esta función es O(N^2) con N siendo la cantidad de partículas

No repetir trabajo

La primera aproximación a mejorar estos 2 problemas es muy simple e implica evitar hacer 2 veces los mismos cálculos. En cada iteración de la función, la misma partícula va a tomar la posición de p1 y luego p2 así que vamos a ejecutar, eventualmente p1.dist2(p2) y luego p2.dist2(p1). Evitemos esto:

void update_density(std::vector<Particle>& particles) {
    for (auto ix = 0u; ix < particles.size(); ix++) {
        auto& p1 = particles[ix];
        for (auto jx = ix + 1; jx < particles.size(); jx++) {
            auto& p2 = particles[jx];
            if (p1.dist2(p2) < THRESHOLD_DIST) {
                p1.density++;
                p2.density++;
            }
        }
    }
}

Con este pequeño cambio nos ahorramos el realizar muchos cálculos 2 veces, porque cada iteración actualiza la densidad de las dos partículas (y evitamos calcular de nuevo ambos valores comparando cada partícula solo con las siguientes en la lista).

Sin embargo, el orden de ejecución no cambia demasiado, ahora realizamos comparaciones en base a la siguiente serie:


El orden sigue siendo cuadrático. Necesitamos una forma de reducir la cantidad de cálculos que realizamos drásticamente.

Localidad del problema

Tratamos de visualizar la densidad de partículas en cada posición de cada partícula, lo que implica que necesitamos una valor de umbral lo sificientemente pequeño para que tenga algún sentido en la imagen final. Como tratamos de dibujar las partículas en la superficie de la esfera unitaria (una unidad de radio), la distancia máxima entre 2 partículas va a ser de 2 unidades (diámetro). Elegí como valor de umbral sqrt(0.004) ~ 0.063 (las pruebas se hacen con todas las distancias cuadradas para ahorrar el costo de realizar la raíz cuadrada), podría haber sido cualquier número, pero este resultó bien desde el principio.

Si uno compara ambos valores, se nota que muchas, MUCHAS partículas están muy lejos entre ellas, por lo que se puede utilizar una estructura de datos que sea útil para descartar la mayor cantidad posible de comparaciones aún antes de realizarlas. Estructuras de este tipo entran en la familia de árboles de división espacial como quadtrees, octrees, árboles BSP y kd-trees por nombrar algunos.
De hecho, en ray tracing usando la técnica de Photon Mapping, es muy común utilizar kd-trees como estructura básica del Photon Map por la calidad de las consultas de puntos por rangos de distancia (al consultar por los fotones en el área circundante a algún punto en el espacio).

El tema, es que se puede tornar costos mantener el kd-tree si las posiciones se agregan y quitan dinámicamete (el árbol termina necesitando re-balanceos muchas veces), como yo quería animar las partículas eventualmente, esto se tornó una inconveniencia.

Lo que terminé implementando es una versión hackeada-trucha, de un octree con tamaño de hojas fijo (cualquiera me puede argumentar que esto se podría haber licenciado bajo CRAPL en lugar de MIT, pero che, funcionó y fue bastante fácil ;P). En lugar de crear dinámicamente la estructura del árbol, definí 2 parámetros, la bounding box del modelo (solo a través de los valores máximo y mínimo que cualquier coordenada pueda tomar) y el número de intervalos en los que cada eje se va a subdividir.

La estructura básica es un hashmap de vectores (sets podrían haber funcionado también, capaz mejor…) donde cada key es el identificador generado a partir de la posición y el valor es un vector con todos los elementos que pertenecen a una hoja de esta cosa.

En lugar de comparar todas las partículas, una gran porción de las mismas es descartada sin problema y solo los mejores candidatos se utilizan en el proceso de comparación.

SPP o “Space Partitioned Positions”

Esta es la estructura aceleradora que utilicé:

class spp {
public:
    spp(uint8_t intervals_per_axis, float min_val, float max_val);
    ~spp() = default;

    uint32_t add(const glm::vec3& pos, size_t external_idx);
    void remove(const glm::vec3& pos, size_t external_idx);
    void remove(uint32_t bucket_id, size_t external_idx);
    std::vector<uint32_t> get_buckets_area(const glm::vec3& pos) const;
    std::vector<uint32_t> get_buckets_area(uint32_t bucket_id) const;
    const std::vector<size_t>& get_bucket(uint32_t bucket_id) const;
private:
    std::unordered_map<uint32_t, std::vector<size_t>> m_buckets;
    uint8_t m_intervals_per_axis;
    glm::vec3 m_min_vec;
    float m_normalize_value;

    uint32_t get_bucket(const glm::vec3& pos) const;
};

Paso a paso:

spp::spp(const uint8_t intervals_per_axis, const float min_val, const float max_val)
    : m_intervals_per_axis(intervals_per_axis)
    , m_min_vec{ min_val, min_val, min_val }
    , m_normalize_value(max_val - min_val) {
    SPL_ASSERT(m_normalize_value != 0., "SPP interval can't be zero.");
    SPL_ASSERT(max_val > min_val, "Consider swapping min and max!");
}

Este es el cosntructor, inicializa algunos miembros auxiliares y hace algunos controles simples. Estos miembros son un vector (solo para simplificar sintaxis, bien podría haber sido un valor solo) para el valor mínimo, otro para la cantidad de intervalos (con un máximo de 0xFF o 255 intervalos por eje) y un tercer valor para el largo de la boundinb box (por eje). (Por como funciona, tanto 0 como 1 en el número de intervalos resultan en el mismo efecto).

uint32_t spp::add(const glm::vec3& pos, const size_t external_idx) {
    const auto bid = get_bucket(pos);
    m_buckets[bid].push_back(external_idx);
    return bid;
}

Este método agrega un elemento a la estructura, este elemento es referencioado con el external_idx (en mi caso estaba utilizando índices para referenciar cada partícula afuera) y la variable pos refiere a la posición de dicho elemento.

Como se ve, genera una key asociada al vector donde se almacenan los elementos usando la posición, guarda el id en el vector y finalmente retorna el identificador de este vector. Dado que calcular el id a partir de la posición es algo costoso, en el resto del código lo optimicé un poco almacenando este identificador en cada partícula. De todos modos esto no es necesario ya que siempre se puede obtener el identificador a partir de la posición.

void spp::remove(const glm::vec3& pos, const size_t external_idx) {
    remove(get_bucket(pos), external_idx);
}
void spp::remove(const uint32_t bucket_id, const size_t external_idx) {
    auto& b = m_buckets[bucket_id];
    b.erase(std::find(b.begin(), b.end(), external_idx));
}

Luego vienen los métodos de quitar elementos, uno utilizando la posición y el otro directamente el identificador de la hoja/vector, Esta es la operación en la que uno podría argumentar que un set sería mejor ‘hoja’ que un vector, pero no he realizado ninguna comparación de performance como para justificar el cambio (otras partes del código podrían funcionar peor con un set, a causa de cache misses y localidad de datos).

Bueno, finalmente el método intersante, get_bucket:

uint32_t spp::get_bucket(const glm::vec3& pos) const {
    const auto norm = (pos - m_min_vec) / m_normalize_value;
    return pack(
        static_cast<uint8_t>(std::floor(m_intervals_per_axis * norm.x)),
        static_cast<uint8_t>(std::floor(m_intervals_per_axis * norm.y)),
        static_cast<uint8_t>(std::floor(m_intervals_per_axis * norm.z)));
}

Pueden ver lo que se hace, cada coordenada es mapeada al intervalo [0, 1] usando los parámetros del constructor y luego el valor obtenido se utiliza para generar un uint8_t por cada coordenada basado en la cantidad de intervalos. Finalmente los 3 valores se empacan en un solo uint32_t usando el método pack.

¿Cómo funciona? Pensémoslo en una sola dimensión:

Usando 8 como el número de intervalos (intervals_per_axis), supongamos que recibimos pos = -0.3:


Así que el identificador de la hoja es 2.
Haciendo lo mismo para 0.5 tenemos: id = 6. Y finalmente, usando un valor muy cercano a 0.5 deberíamos obtener o bien el mismo identificador o uno contiguo, para 0.55 -> 6, para 0.45 -> 5, para 0.75 -> 7. Estamos dividiendo el intervalo original en otros más pequeños cuyo largo es igual al largo de la bounding box (por eje) dividido por el número de intervalos y asignando identificadores basados en dentro de cuál subintervalo cae.

Si alguna vez leyeron sobre scalar quantization o empaquetado de floats (como para los canales de colores) esto es lo mismo.

El resultado general es que:

  • Es un hecho que la distancia máxima entre elementos dentro de la misma hoja es igual a:
  • Desafortunadamente eso no significa que una sola hoja es suficiente para calcular la densidad alrededor de cada partícula, dado que 2 partículas muy cercanas podrían terminar en hojas aledañas pero…
  • Se pueden obtener todas las partículas que están dentro de un radio X utilizando la hoja en la que cayó la partícula y todas las hojas vecinas, siempre y cuando X sea menor que m_normalize_value/m_intervals_per_axis.

En nuestro caso unidimensional, la distancia máxima entre 2 valores dentro de la misma hoja es 0.25, y tomando cualquier posición P dentro de la hoja H, se pueden calcular todas las posiciónes con una distancia a P menor a 0.25 analizando solo las partículas en las hojas B, B-1 (si B != 0) y B+1 (si B != 7). Al dividir nuestro set original en 8, se reducen la cantidad de posiciones a analizar a aquellas que estén dentro de 3 de las 8 hojas (en el peor caso).

Prosiguiendo:

std::vector<uint32_t> spp::get_buckets_area(const uint32_t bucket_id) const {
    std::vector<uint32_t> buckets;
    buckets.reserve(3 * 3 * 3); // Max adjacent buckets (think of a rubik's cube)
    const auto unpacked = unpack(bucket_id);
    for (int8_t ix = -1; ix < 2; ix++) {
        const int8_t xval = unpacked[0] + ix;
        if (xval >= 0 && xval < m_intervals_per_axis) {
            for (int8_t iy = -1; iy < 2; iy++) {
                const int8_t yval = unpacked[1] + iy;
                if (yval >= 0 && yval < m_intervals_per_axis) {
                    for (int8_t iz = -1; iz < 2; iz++) {
                        const int8_t zval = unpacked[2] + iz;
                        if (zval >= 0 && zval < m_intervals_per_axis) {
                            const auto potential_id = pack(xval, yval, zval);
                            if (m_buckets.find(potential_id) != m_buckets.end()) {
                                buckets.push_back(potential_id);
                            }
                        }
                    }
                }
            }
        }
    }
    return buckets;
}

Este código nos permite detectar el set de hoja+vecinas calculando los valores válidos para cada parte del identificador empaquetado (imaginen el centro de un cubo rubik, y las posiciones de todas las restantes piezas). Dado un solo identificador obtenemos ese mismo ID y todos los de los 26 vecinos (en el peor caso). Finalmente, pack y unpack son métodos que empacan 3 (podrían ser 4) enteros de 8 bits en uno de 32.

static uint32_t pack(const uint8_t x, const uint8_t y, const uint8_t z) {
    return  x << 16 | y << 8 | z;
}
static std::array<uint8_t, 3> unpack(uint32_t value) {
    return{
        (value >> 16) & 0xFF,
        (value >> 8) & 0xFF,
        value & 0xFF
    };
}

Usando el SPP

Una vez que el SPP está pronto, en lugar de computar la densidad calculando las distancias entre todas las partículas, podemos:

  • Agregar todas las partículas al SPP y luego,
  • Para cada partícula a actualizar obtenemos el set de hojas que son vecinas y,
  • Calcular la distancia a las partículas dentro de esas hojas únicamente:
void load_particles(std::vector<Particle>& particles) {
    for (auto ix = 0u; ix < particles.size(); ix++) {
        // Note we store the bucket id outside to reduce the amount of times we need to calculate it later.
        p.bucket_id = spp.add(p.position, ix);
    }
}

void update_density(std::vector<Particle>& particles) {
    for (auto ix = 0u; ix < particles.size(); ix++) {
        auto& p1 = particles[ix];
        auto buckets = spp.get_buckets_area(p1.bucket_id);
        for (auto jx = 0u; jx < buckets.size(); jx++) {
            auto& bucket = spp.get_bucket(buckets[jx]);
            for (auto kx = 0u; kx < bucket.size(); kx++) {
                auto& p2 = particles[bucket[kx]];
                if (p1.dist2(p2) < THRESHOLD_DIST) {
                    p1.density++;
                }
            }
        }
    }
}

Noten que nuevamente estoy repitiendo algo de trabajo, porque ya no actualizo la densidad de ambas partículas en la misma iteración, esto se debe a que ya no se puede asegurar el orden en el que las partículas van a ser visitadas, así que tratar de evitar este doble cálculo se vuelve más complicado. El ahorro de utilizar el spp justifica más que suficientemente por esto (y por la posible complejidad del algoritmo para saltearse correctamente las comparaciones ya realizadas).

Luego de esto, calcular el color de cada partícula es solo cosa de encontrar la densidad máxima y obtener un valor de RGBA basado en la proporción de densidad de cada partícula sobre el máximo. Yo utilicé algo parecido a esto:

Color color_for_particle(Particle& p, const uint32_t max_density) {
    float alive_value = p.isAlive() ? 1.f : 0.f; //Se podría castear simplemete 🙂
    return Color {
        alive_value,
        (p.density / max_density) * alive_value,
        (p.density / max_density) * alive_value,
        alive_value
    };
}

Si analizan la historia del repositorio, pueden ver que en un principio calculaba el color en el CPU y luego se forwardeaba directamente al shader. Obviamente esto no es la mejor opción, en la versión actual, solo el valor máximo se calcula en el CPU y luego es pasado como un uniform, y la densidad por partícula se pasa como otro atributo. De esta manera, el cálculo se puede realizar en el GPU.

Azul, Rojo y la distancia al centro

En esta demo, los cálculos son más simples, nada se realiza en el CPU y el shader genera rojo si la distancia de la partícula al centro es menor a 1 y azul si no.

Color color_for_particle_blue(Particle& p) {
    float alive_value = p.isAlive() ? 1.f : 0.f; //Se podría castear simplemete 🙂
    float r = (p.length <= 1.) ? 1.f : 0.f; //Se podría castear simplemete 🙂
    return Color {
        r,
        0.f,
        1.f - r,
        alive_value
    };
}

Nota al pie sobre el número de intervalos en el SPP

El número de hojas influencia directamente la cantidad de comparaciones que debemos realizar pero además la distancia máxima que dos partículas pueden tener. Necesitamos instanciar el SPP con el valor correcto para saltear el mayor número de comparaciones entre partículas sin perdernos ninguna de las partículas que deben ser parte del cálculo de densidad. Para caulcular esta densidad hay que ‘encerrar’ el volumen de cercanía alrededor de la partícula en el cubo creado por la hoja a la que la partícula pertenece y todas las hojas vecinas:

  • Habiendo elegido un umbral igual a srqr(0.004) ~ 0.0632 es necesario un tamaño por hoja (por eje) por lo menos de ese valor.
  • Sabemos que el valor máximo de coordenadas es uno y el mínimo es -1, así que la bounding box tiene 2 unidades de largo por eje





Este número de intervalos por eje nos da un tamaño de hoja (por eje) de ~0.645 unidades, el valor más pequeño no menor al umbral. Esto asegura que aún el peor caso queda cubierto, y todas las particulas que pueden contribuir al valor de densidad de una partícula dada son evaluadas.

El peor caso, está dado por una partícula que cae en el plano de borde de una hoja. Pensando nuevamente en el caso unidimensional, si el umbral fuese 0.26, tener 8 intervalos es demasiado: una muestra con la posición 0.5 termina en la hoja 6, y una muestra con la posición 0.245 termina en la hoja 4. Estas dos muestras no se van a comparar porque no están en hojas vecinas aún cuando la distancia entre ellas es menor a 0.26.

¡Fin! ¡Bien! Terminaste un artícula largo sobre los detalles de una demo sin importancia, ¡felicidades! Por lo menos espero que le hayas sacado algún provecho. Podés referirte al artículo previo para ver la historia original y el código fuente.
Inicio

Advertisements

One thought on “Particles! (2/?)

Leave a Reply (Deja una respuesta)

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s