Automatas Celulares para simulación de explosiones en juegos

Ayer me puse a implementar un Autómata Celular. El objetivo final es simular una explosión en tiempo real. Voy a mostrar unos resultados preliminares que obtuve con Python, pero primero explico la idea.

En general, simular fenómenos físicos como explosiones, fluidos y fuego, en forma realista es demasiado costoso numéricamente como para que se pueda hacer en tiempo real en un videojuego. Sin embargo, hay una técnica que permite obtener aproximaciones visualmente plausibles utilizando los llamados Autómatas Celulares (AC).

Autómatas Celulares

Un Autómata Celular se puede pensar como una grilla cuyas celdas tienen asociado uno o varios valores. Cada celda tiene vecinos, generalmente 4. Al principio, se asignan valores iniciales, y luego se va haciendo evolucionar el AC siguiendo una regla matemática predefinida. Con este esquema general se pueden simular muchos fenómenos biológicos y físicos interesantes. Ver por ejemplo una implementación en Flash

Este concepto se puede generalizar mediante grafos dinámicos. Un grafo es en conjunto de Vértices, unidos mediante conexiones llamadas Aristas. Tanto los Vértices como las Aristas pueden tener propiedades asignadas.

 Ejemplo de grafo: mis contactos de Facebook. Cada Vértice (círculos) es un amigo. Cada Arista (lineas) indica que mis contactos se tienen en sus respectivas listas de contactos. También hay propiedades asociadas: cada Arista tiene un tamaño dado por la actividad del contacto, una foto y un color que indica que el contacto pertenece a un "cluster". Los lados tienen un numero relacionado que indica la cantidad de fotos en la que ambos aparecen juntos.

Usando ese concepto general, se puede generalizar un Autómata Celular. El ejemplo básico de una grilla 2D, se traslada a un grafo en que cada Vértice tiene una posición (lo que originalmente era una coordenada en la grilla) y cuatro Aristas que conectan al Vértice con sus vecinos.

Esquema básico de un grafo que reemplaza a una grilla. Cada vértice de las paredes está conectado con un vértice de la pared opuesta, para crear un "mundo infinito". Todo lo que "sale" por la derecha "entra" por la izquierda.

Este esquema es mas general en varios sentidos. Primero, permite desacoplar fácilmente el espacio "real" de nuestro Autómata Celular, de esa forma ahora cada nodo del grafo puede estar en posiciones arbitrarias del espacio y, en particular, ubicados en forma no uniforme. Se puede entonces simular con mas detalle la física en algunas partes que en otras. Solamente es cuestión de ser cuidadoso a la hora de crear nuestro algoritmo. Otra virtud que tiene este esquema es que se pueden generar obstrucciones en el grafo para simular obstrucciones reales. Por ejemplo

Simular obstrucciones es solo cuestión de romper Aristas o Vértices. La imagen muestra cómo se puede relacionar un grafo con un mapa de nuestro juego. En esquemas mas complicados se pueden mantener todos los Vértices y Aristas y asignar pesos a las Aristas. De esta forma se pueden simular distintos tipos de "rigidez" o incluso destrucción de paredes.

Simulando explosiones

En este post explico como usar este esquema para simular explosiones plausibles en un juego, con poco costo numérico. El primer código que uso es python puro, con fines demostrativos. Incluso con la lentitud de implementar un algoritmo que recorre nodos en python, es posible ejecutar en tiempo real una grilla de 30x30 que ya presenta resultados interesantes.
Para simular una explosión se necesita un grafo que tenga definida una presión en cada Vértice. En cada paso de tiempo, se mira Vértice por Vértice y se compara la presión del Vértice con la presión de sus vecinos. Luego se calcula un flujo de "aire" entre el Vértice y el vecino y se lo utiliza para modificar sus presiones.
El esquema es básico e intuitivo: El nodo que tiene mas presión le entrega un poco a cada uno de los vecinos que tienen menos. Y así, paso por paso, se genera una frente de presión.

Voy a mostrar los resultados primero, y dejar el código comentado al final para quien quiera leerlo.

Explosión al lado de un obstáculo pequeño. Los tiempos son t=1, 4, 10, 20, 40, 80 y 200.

Explosión en una "puerta". La linea marrón es una pared. La simulación es para tiempos t=1, 4, 8, 20, 100 y 200.

Trasladando los efectos de la explosión al juego

Supongamos que el objeto que nos interesa se encuentra cerca de la explosión. Podemos calcular el vector de flujo sumando los vectores de flujo entre un Vértice y sus vecinos. En el ejemplo de la imagen siguiente, el vector de flujo sobre el nodo en que se encuentra el personaje tiene dirección sur oeste. Ese flujo, multiplicado por algún valor que represente la sección de área del personaje, nos indica una fuerza que actúa sobre el personaje.

Es interesante notar que en una explosión, todo sucede muy rápido en los primeros instantes de tiempo. En un juego uno podría simular decenas de pasos de un CA por cada frame del juego y utilizar un promedio en el tiempo del vector de flujo. Este promedio es el impulso a transferir a los objetos del juego en cada frame. De esta forma se puede simular el efecto medio de la onda expansiva en los jugadores, objetos ubicados en el mapa y en el mapa mismo.

También se puede sumar las magnitudes de flujo, para simular efectos como sordera momentanea, ruptura de vidrios, o daño general.

 Cómo trasladar el efecto de la onda expansiva al personaje: Las flechas azules son los vectores de flujo entre el nodo del personaje y sus vecinos. La flecha roja es la suma de los vectores azules, que multiplicadas por el cross-section del personaje dan la fuerza que actúa sobre el personaje debido a la explosión. Promediando esta fuerza durante varios pasos del AC se obtiene el impulso sobre el personaje.

Código Python

Este es el código de los dos ejemplos mostrados. El próximo paso es convertir el loop principal a código C y exponer funciones para calcular los impulsos transmitidos a objetos del mapa.

def iterate_fl(N=30, steps=10):
    '''Automata Celular simple para calcular una explosión.
    Implementación de algoritmo de Tom Forsyth.
    Mas detalles:
    http://home.comcast.net/~tom_forsyth/papers/cellular_automata_for_physical_modelling.html'''

    from itertools import product, combinations
    def clamp(value, max, min):
        '''
        Función clamp, devuelve value si value está entre min y max, devuelve min o max si no.
        '''
        if value<min:
            return min
        if value>max:
            return max
        return value
   
    class Node():
        '''
        Decidí usar una implementación python pura para las pruebas preliminares.
        Este modelo permite CA en 3D y simular paredes u objetos inmoviles mediante desconexiones de nodos.
        '''
        pressure=0.0
        new_pressure=0.0
        t=0
       
        def __init__(self, position):
            self.position=position
            '''Los grafos tienen un conjunto de vecinos.'''
            self.neighbours=set()
           
        def __repr__(self):
            return 'node pressure: %s'%self.pressure
   
   
    '''Mantengo un diccionario posicion->nodo para poder acceder los nodos facilmente'''
    nodes=dict([(position,Node(position)) for position in product(range(N),range(N))])
   
    '''Cada nodo representa una celda en un espacio 2D, cuyo vecino es el vecino físico'''
    for node in nodes.values():
        '''En este ejemplo muestro dos tipos de condiciones de contorno en un grafo:
        paredes rígidas y condiciones periodicas'''
        if node.position[0]+1<N-1:
            '''Hay paredes rigidas a la derecha del mapa.'''
            nodes[node.position].neighbours.add(nodes[(node.position[0]+1, node.position[1])])
        if node.position[0]-1>0:
            '''Hay paredes rígidas a la izquierda del mapa'''
            nodes[node.position].neighbours.add(nodes[(node.position[0]-1, node.position[1])])

        '''La funcion mod (%) me asegura que la parte superior e inferior del mapa están
        conectadas. Lo que sale por arriba, aparece por abajo.
        '''
        nodes[node.position].neighbours.add(nodes[(node.position[0],(node.position[1]+1)%N)])
        nodes[node.position].neighbours.add(nodes[(node.position[0],(node.position[1]-1)%N)])
   
    '''Crea pared, removiendo nodos del grafo. Se pueden hacer paredes removiendo solamente conexiones
    Tipos de pared:'''
    wall=product((10,),range(0,15)+range(17,N)) #pared que cruza el mapa en x=10, con un agujero
    #wall=product((10,),range(14,17)) #pared de tamaño limitado en x=10
    for edg_pos in wall:
        node=nodes[edg_pos]
        for neigh in node.neighbours:
            neigh.neighbours.remove(node)
        del nodes[edg_pos]
   
    '''La explosión propiamente dicha.'''
    #explosion=product(range(4,10),(5,15)) # Un rectangulo
    explosion=((9,15),) #Solo un punto
    for node in explosion:
        nodes[node].pressure=300.0 #La explosion son puntos de alta presion
        nodes[node].new_pressure=300.0 #Este es un detalle del algoritmo
   
       
    '''Ejecutar CA'''   
    i=0
    while i<steps:
        i=i+1
        '''Algoritmo CA para explosiones'''
        for node in nodes.values():
            '''
            Recorremos todo el grafo
            '''
            if node.t!=i:
                '''Detalle de implementación: iteramos por todo el grafo aplicando las reglas del CA,
                los calculos se realizan con el estado actual del grafo y se almacenan los resultados en
                new_pressure, en el siguiente paso de tiempo, se actualizan los valores de la presion usando
                los valores calculados en el paso anterior.'''
                node.t=i
                node.pressure=node.new_pressure
               
           
            for neight in node.neighbours:
                '''Para cada nodo, se miran los vecinos. Este ejemplo es 2D, pero se puede hacer
                3D. En 2D tenemos 4 vecinos por cada posicion del espacio, pero puede haber menos
                vecinos si hay una obstrucción.'''
                if neight.t!= i:
                    '''Nos aseguramos de estar usando los valores de presion
                    calculados en el paso anterior'''
                    neight.t=i
                    neight.pressure=neight.new_pressure
                   
                '''Se calcula la dif. de presion entre el nodo actual y su vecino'''
                dpres=node.pressure-neight.pressure
               
                '''Calculamos el flujo como una fracción de la diferencia de presion.
                0.12 es un numero magico que indica la fluidez del medio. Un valor muy chico
                hace que la presion se propague mas lentamente.
                Usamos la funcion clamp para evitar valores negativos en el nodo.
                Se hace que el flujo no produsca una presión negativa'''
                flow=clamp(0.12*dpres, node.pressure/4.0, -neight.pressure/4.0)
               
                #Actualizamos la presion (usando new_pressure en lugar de pressure) substrayendo
                #presion del nodo actual y agregandola a su vecino. Un flujo negativo implica
                #que la presion del nodo actual aumenta y la del vecino disminuye.
                node.new_pressure -= flow
               
                neight.new_pressure += flow
       
    return nodes