python algorithm markov-chains procedural-generation

Problemas al implementar el algoritmo "Función de colapso de onda" en Python



algorithm markov-chains (2)

En una palabra:

Mi implementación del algoritmo de la función de colapso de onda en Python 2.7 es defectuosa, pero no puedo identificar dónde se encuentra el problema. Necesitaría ayuda para descubrir lo que posiblemente me estoy perdiendo o haciendo mal.

¿Qué es el algoritmo de la función de colapso de onda ?

Es un algoritmo escrito en 2016 por Maxim Gumin que puede generar patrones de procedimiento a partir de una imagen de muestra. Puede verlo en acción here (modelo de superposición 2D) y here (modelo de mosaico 3D).

Objetivo de esta implementación:

Reducir el algoritmo (modelo superpuesto 2D) a su esencia y evitar las redundancias y la torpeza del guión original de C # (sorprendentemente largo y difícil de leer). Este es un intento de hacer una versión más corta, más clara y pitónica de este algoritmo.

Características de esta implementación:

Estoy usando Processing (modo Python), un software para diseño visual que facilita la manipulación de imágenes (sin PIL, sin Matplotlib, ...). Los principales inconvenientes son que estoy limitado a Python 2.7 y NO puedo importar numpy.

A diferencia de la versión original, esta implementación:

  • no está orientado a objetos (en su estado actual), lo que hace que sea más fácil de entender / más cercano al pseudocódigo
  • está utilizando matrices 1D en lugar de matrices 2D
  • está utilizando el corte de matriz para la manipulación de matrices

El algoritmo (según lo entiendo)

1 / Lea el mapa de bits de entrada, almacene todos los patrones NxN y cuente sus ocurrencias. ( opcional: aumentar los datos del patrón con rotaciones y reflexiones).

Por ejemplo, cuando N = 3:

2 / Precalcule y almacene todas las relaciones de adyacencia posibles entre patrones. En el siguiente ejemplo, los patrones 207, 242, 182 y 125 pueden superponerse al lado derecho del patrón 246

3 / Cree una matriz con las dimensiones de la salida (llamada W para onda). Cada elemento de esta matriz es una matriz que contiene el estado ( True False ) de cada patrón.

Por ejemplo, supongamos que contamos 326 patrones únicos en la entrada y queremos que nuestra salida tenga dimensiones de 20 por 20 (400 celdas). Luego, la matriz "Wave" contendrá 400 (20x20) matrices, cada una de ellas con 326 valores boolan.

Al inicio, todos los valores booleanos se establecen en True porque cada patrón está permitido en cualquier posición de la Onda.

W = [[True for pattern in xrange(len(patterns))] for cell in xrange(20*20)]

4 / Cree otra matriz con las dimensiones de la salida (llamada H ). Cada elemento de esta matriz es un flotante que contiene el valor de "entropía" de su celda correspondiente en la salida.

La entropía aquí se refiere a la entropía de Shannon y se calcula en función del número de patrones válidos en una ubicación específica en la onda. Cuanto más una celda tenga patrones válidos (establecida en True en la onda), mayor será su entropía.

Por ejemplo, para calcular la entropía de la celda 22, observamos su índice correspondiente en la onda ( W[22] ) y contamos el número de booleanos establecido en True . Con ese recuento ahora podemos calcular la entropía con la fórmula de Shannon. El resultado de este cálculo se almacenará en H en el mismo índice H[22]

Al inicio, todas las celdas tienen el mismo valor de entropía (mismo flotador en cada posición en H ) ya que todos los patrones están establecidos en True , para cada celda.

H = [entropyValue for cell in xrange(20*20)]

Estos 4 pasos son pasos introductorios, son necesarios para inicializar el algoritmo. Ahora comienza el núcleo del algoritmo:

5 / Observación:

Encuentre el índice de la celda con la entropía mínima distinta de cero (tenga en cuenta que en la primera iteración todas las entropías son iguales, por lo que debemos elegir el índice de una celda al azar).

Luego, observe los patrones aún válidos en el índice correspondiente en la onda y seleccione uno de ellos al azar, ponderado por la frecuencia con la que aparece el patrón en la imagen de entrada (opción ponderada).

Por ejemplo, si el valor más bajo en H está en el índice 22 ( H[22] ), observamos todos los patrones establecidos en True en W[22] y elegimos uno aleatoriamente en función del número de veces que aparece en la entrada. (Recuerde que en el paso 1 hemos contado el número de ocurrencias para cada patrón). Esto asegura que los patrones aparezcan con una distribución similar en la salida como se encuentra en la entrada.

6 / colapso:

Ahora asignamos el índice del patrón seleccionado a la celda con la entropía mínima. Lo que significa que todos los patrones en la ubicación correspondiente en la Ola están configurados en False excepto el que se ha elegido.

Por ejemplo, si el patrón 246 en W[22] se configuró en True y se seleccionó, entonces todos los demás patrones se configurarán en False . La celda 22 tiene asignado el patrón 246 . En la celda de salida 22 se rellenará con el primer color (esquina superior izquierda) del patrón 246. (azul en este ejemplo)

7 / Propagación:

Debido a las restricciones de adyacencia, esa selección de patrones tiene consecuencias en las celdas vecinas de la Ola. Las matrices de booleanos correspondientes a las celdas a la izquierda y a la derecha, encima y encima de la celda recientemente colapsada deben actualizarse en consecuencia.

Por ejemplo, si la celda 22 se ha colapsado y asignado con el patrón 246 , entonces W[21] (izquierda), W[23] (derecha), W[2] (arriba) y W[42] (abajo) tienen que modificarse así que solo mantienen True los patrones adyacentes al patrón 246 .

Por ejemplo, mirando hacia atrás en la imagen del paso 2, podemos ver que solo los patrones 207, 242, 182 y 125 se pueden colocar a la derecha del patrón 246. Eso significa que W[23] (a la derecha de la celda 22 ) necesita mantenga los patrones 207, 242, 182 y 125 como True y establezca todos los demás patrones en la matriz como False . Si estos patrones ya no son válidos (ya están configurados en False debido a una restricción previa), entonces el algoritmo enfrenta una contradicción .

8 / Actualización de entropías

Debido a que una celda se ha colapsado (un patrón seleccionado, establecido en True ) y sus celdas circundantes se actualizaron en consecuencia (estableciendo patrones no adyacentes en False ), la entropía de todas estas celdas ha cambiado y debe calcularse nuevamente. (Recuerde que la entropía de una celda está correlacionada con el número de patrones válidos que tiene en la onda).

En el ejemplo, la entropía de la celda 22 ahora es 0, ( H[22] = 0 , porque solo el patrón 246 se establece en True en W[22] ) y la entropía de sus celdas vecinas ha disminuido (patrones que no eran adyacentes al patrón 246 se han establecido en False ).

Por ahora, el algoritmo llega al final de la primera iteración y recorrerá los pasos 5 (buscar celda con entropía mínima distinta de cero) a 8 (actualizar entropías) hasta que todas las celdas se colapsen.

Mi guion

Necesitará Processing con modo Python instalado para ejecutar este script. Contiene alrededor de 80 líneas de código (cortas en comparación con las ~ 1000 líneas del script original) que están completamente anotadas para que pueda entenderse rápidamente. También deberá descargar la imagen de entrada y cambiar la ruta en la línea 16 en consecuencia.

from collections import Counter from itertools import chain, izip import math d = 20 # dimensions of output (array of dxd cells) N = 3 # dimensions of a pattern (NxN matrix) Output = [120 for i in xrange(d*d)] # array holding the color value for each cell in the output (at start each cell is grey = 120) def setup(): size(800, 800, P2D) textSize(11) global W, H, A, freqs, patterns, directions, xs, ys, npat img = loadImage(''Flowers.png'') # path to the input image iw, ih = img.width, img.height # dimensions of input image xs, ys = width//d, height//d # dimensions of cells (squares) in output kernel = [[i + n*iw for i in xrange(N)] for n in xrange(N)] # NxN matrix to read every patterns contained in input image directions = [(-1, 0), (1, 0), (0, -1), (0, 1)] # (x, y) tuples to access the 4 neighboring cells of a collapsed cell all = [] # array list to store all the patterns found in input # Stores the different patterns found in input for y in xrange(ih): for x in xrange(iw): '''''' The one-liner below (cmat) creates a NxN matrix with (x, y) being its top left corner. This matrix will wrap around the edges of the input image. The whole snippet reads every NxN part of the input image and store the associated colors. Each NxN part is called a ''pattern'' (of colors). Each pattern can be rotated or flipped (not mandatory). '''''' cmat = [[img.pixels[((x+n)%iw)+(((a[0]+iw*y)/iw)%ih)*iw] for n in a] for a in kernel] # Storing rotated patterns (90°, 180°, 270°, 360°) for r in xrange(4): cmat = zip(*cmat[::-1]) # +90° rotation all.append(cmat) # Storing reflected patterns (vertical/horizontal flip) all.append(cmat[::-1]) all.append([a[::-1] for a in cmat]) # Flatten pattern matrices + count occurences '''''' Once every pattern has been stored, - we flatten them (convert to 1D) for convenience - count the number of occurences for each one of them (one pattern can be found multiple times in input) - select unique patterns only - store them from less common to most common (needed for weighted choice)'''''' all = [tuple(chain.from_iterable(p)) for p in all] # flattern pattern matrices (NxN --> []) c = Counter(all) freqs = sorted(c.values()) # number of occurences for each unique pattern, in sorted order npat = len(freqs) # number of unique patterns total = sum(freqs) # sum of frequencies of unique patterns patterns = [p[0] for p in c.most_common()[:-npat-1:-1]] # list of unique patterns sorted from less common to most common # Computes entropy '''''' The entropy of a cell is correlated to the number of possible patterns that cell holds. The more a cell has valid patterns (set to ''True''), the higher its entropy is. At start, every pattern is set to ''True'' for each cell. So each cell holds the same high entropy value'''''' ent = math.log(total) - sum(map(lambda x: x * math.log(x), freqs)) / total # Initializes the ''wave'' (W), entropy (H) and adjacencies (A) array lists W = [[True for _ in xrange(npat)] for i in xrange(d*d)] # every pattern is set to ''True'' at start, for each cell H = [ent for i in xrange(d*d)] # same entropy for each cell at start (every pattern is valid) A = [[set() for dir in xrange(len(directions))] for i in xrange(npat)] #see below for explanation # Compute patterns compatibilities (check if some patterns are adjacent, if so -> store them based on their location) '''''' EXAMPLE: If pattern index 42 can placed to the right of pattern index 120, we will store this adjacency rule as follow: A[120][1].add(42) Here ''1'' stands for ''right'' or ''East''/''E'' 0 = left or West/W 1 = right or East/E 2 = up or North/N 3 = down or South/S '''''' # Comparing patterns to each other for i1 in xrange(npat): for i2 in xrange(npat): for dir in (0, 2): if compatible(patterns[i1], patterns[i2], dir): A[i1][dir].add(i2) A[i2][dir+1].add(i1) def compatible(p1, p2, dir): ''''''NOTE: what is refered as ''columns'' and ''rows'' here below is not really columns and rows since we are dealing with 1D patterns. Remember here N = 3'''''' # If the first two columns of pattern 1 == the last two columns of pattern 2 # --> pattern 2 can be placed to the left (0) of pattern 1 if dir == 0: return [n for i, n in enumerate(p1) if i%N!=2] == [n for i, n in enumerate(p2) if i%N!=0] # If the first two rows of pattern 1 == the last two rows of pattern 2 # --> pattern 2 can be placed on top (2) of pattern 1 if dir == 2: return p1[:6] == p2[-6:] def draw(): # Equivalent of a ''while'' loop in Processing (all the code below will be looped over and over until all cells are collapsed) global H, W, grid ### OBSERVATION # Find cell with minimum non-zero entropy (not collapsed yet) ''''''Randomly select 1 cell at the first iteration (when all entropies are equal), otherwise select cell with minimum non-zero entropy'''''' emin = int(random(d*d)) if frameCount <= 1 else H.index(min(H)) # Stoping mechanism '''''' When ''H'' array is full of ''collapsed'' cells --> stop iteration '''''' if H[emin] == ''CONT'' or H[emin] == ''collapsed'': print ''stopped'' noLoop() return ### COLLAPSE # Weighted choice of a pattern '''''' Among the patterns available in the selected cell (the one with min entropy), select one pattern randomly, weighted by the frequency that pattern appears in the input image. With Python 2.7 no possibility to use random.choice(x, weight) so we have to hard code the weighted choice '''''' lfreqs = [b * freqs[i] for i, b in enumerate(W[emin])] # frequencies of the patterns available in the selected cell weights = [float(f) / sum(lfreqs) for f in lfreqs] # normalizing these frequencies cumsum = [sum(weights[:i]) for i in xrange(1, len(weights)+1)] # cumulative sums of normalized frequencies r = random(1) idP = sum([cs < r for cs in cumsum]) # index of selected pattern # Set all patterns to False except for the one that has been chosen W[emin] = [0 if i != idP else 1 for i, b in enumerate(W[emin])] # Marking selected cell as ''collapsed'' in H (array of entropies) H[emin] = ''collapsed'' # Storing first color (top left corner) of the selected pattern at the location of the collapsed cell Output[emin] = patterns[idP][0] ### PROPAGATION # For each neighbor (left, right, up, down) of the recently collapsed cell for dir, t in enumerate(directions): x = (emin%d + t[0])%d y = (emin/d + t[1])%d idN = x + y * d #index of neighbor # If that neighbor hasn''t been collapsed yet if H[idN] != ''collapsed'': # Check indices of all available patterns in that neighboring cell available = [i for i, b in enumerate(W[idN]) if b] # Among these indices, select indices of patterns that can be adjacent to the collapsed cell at this location intersection = A[idP][dir] & set(available) # If the neighboring cell contains indices of patterns that can be adjacent to the collapsed cell if intersection: # Remove indices of all other patterns that cannot be adjacent to the collapsed cell W[idN] = [True if i in list(intersection) else False for i in xrange(npat)] ### Update entropy of that neighboring cell accordingly (less patterns = lower entropy) # If only 1 pattern available left, no need to compute entropy because entropy is necessarily 0 if len(intersection) == 1: H[idN] = ''0'' # Putting a str at this location in ''H'' (array of entropies) so that it doesn''t return 0 (float) when looking for minimum entropy (min(H)) at next iteration # If more than 1 pattern available left --> compute/update entropy + add noise (to prevent cells to share the same minimum entropy value) else: lfreqs = [b * f for b, f in izip(W[idN], freqs) if b] ent = math.log(sum(lfreqs)) - sum(map(lambda x: x * math.log(x), lfreqs)) / sum(lfreqs) H[idN] = ent + random(.001) # If no index of adjacent pattern in the list of pattern indices of the neighboring cell # --> mark cell as a ''contradiction'' else: H[idN] = ''CONT'' # Draw output '''''' dxd grid of cells (squares) filled with their corresponding color. That color is the first (top-left) color of the pattern assigned to that cell '''''' for i, c in enumerate(Output): x, y = i%d, i/d fill(c) rect(x * xs, y * ys, xs, ys) # Displaying corresponding entropy value fill(0) text(H[i], x * xs + xs/2 - 12, y * ys + ys/2)

Problema

A pesar de todos mis esfuerzos por poner cuidadosamente en código todos los pasos descritos anteriormente, esta implementación arroja resultados muy extraños y decepcionantes:

Ejemplo de una salida de 20x20

Parece que se respetan tanto la distribución del patrón como las restricciones de adyacencia (la misma cantidad de colores azul, verde, amarillo y marrón que en la entrada y el mismo tipo de patrones: suelo horizontal, tallos verdes).

Sin embargo, estos patrones:

  • a menudo están desconectados
  • a menudo están incompletos (falta de "cabezas" compuestas de pétalos de 4 amarillos)
  • encuentra demasiados estados contradictorios (celdas grises marcadas como "CONT")

En ese último punto, debo aclarar que los estados contradictorios son normales pero que deben ocurrir muy raramente (como se indica en el medio de la página 6 de this documento y en this artículo)

Horas de depuración me convencieron de que los pasos introductorios (1 a 5) son correctos (contar y almacenar patrones, cálculos de adyacencia y entropía, inicialización de matrices). Esto me ha llevado a pensar que algo debe estar mal con la parte central del algoritmo (pasos 6 a 8) . O estoy implementando uno de estos pasos incorrectamente o me falta un elemento clave de la lógica.

¡Cualquier ayuda con respecto a este asunto sería inmensamente apreciada!

Además, cualquier respuesta que se base en la secuencia de comandos proporcionada (usando Procesamiento o no) es bienvenida .

Recursos adicionales útiles:

Este this detallado de Stephen Sherratt y este this explicativo de Karth & Smith. Además, para comparar, sugeriría verificar esta otra implementación de Python (contiene un mecanismo de retroceso que no es obligatorio).

Nota: hice todo lo posible para que esta pregunta fuera lo más clara posible (explicación completa con GIF e ilustraciones, código completamente anotado con enlaces útiles y recursos), pero si por alguna razón decide rechazarla, deje un breve comentario para explicar por qué lo estás haciendo


Mientras mira la demostración en vivo vinculada en uno de sus ejemplos, y basada en una revisión rápida del código del algoritmo original, creo que su error radica en el paso "Propagación".

La propagación no es solo actualizar las 4 celdas vecinas a la celda contraída. También debe actualizar todos esos vecinos de celdas, y luego los vecinos de esas celdas, etc., de forma recursiva. Bueno, para ser específicos, tan pronto como actualice una sola celda vecina, entonces actualice su vecina (antes de llegar a los otros vecinos de la primera celda), es decir, actualizaciones de profundidad, no amplitud. Al menos, eso es lo que deduzco de la demostración en vivo.

La implementación real del código C # del algoritmo original es bastante complicada y no lo entiendo completamente, pero los puntos clave parecen ser la creación del objeto "propagador" here , así como la función Propagate en sí, here .


La hipótesis sugerida por @mbrig y @Leon de que el paso de propagación itera sobre una pila completa de celdas (en lugar de limitarse a un conjunto de 4 vecinos directos) era correcta. El siguiente es un intento de proporcionar más detalles al responder mis propias preguntas.

El problema ocurrió en el paso 7, mientras se propagaba. El algoritmo original actualiza los 4 vecinos directos de una celda específica PERO:

  • el índice de esa celda específica se reemplaza por turnos por los índices de los vecinos actualizados previamente .
  • Este proceso en cascada se activa cada vez que una celda se colapsa
  • y duran mientras los patrones adyacentes de una celda específica estén disponibles en 1 de su celda vecina

En otras palabras, y como se menciona en los comentarios, este es un tipo de propagación recursiva que actualiza no solo a los vecinos de la celda colapsada, sino también a los vecinos de los vecinos ... y así sucesivamente, siempre que sean posibles las adyacencias.

Algoritmo Detallado

Una vez que una celda se contrae, su índice se coloca en una pila. Esa pila está destinada más tarde a almacenar temporalmente índices de celdas vecinas

stack = set([emin]) #emin = index of cell with minimum entropy that has been collapsed

La propagación durará mientras esa pila esté llena de índices:

while stack:

Lo primero que hacemos es pop() el último índice contenido en la pila (el único por ahora) y obtener los índices de sus 4 celdas vecinas (E, W, N, S). Tenemos que mantenerlos dentro de los límites y asegurarnos de que se adapten.

while stack: idC = stack.pop() # index of current cell for dir, t in enumerate(mat): x = (idC%w + t[0])%w y = (idC/w + t[1])%h idN = x + y * w # index of neighboring cell

Antes de continuar, nos aseguramos de que la celda vecina aún no esté contraída (no queremos actualizar una celda que solo tenga 1 patrón disponible):

if H[idN] != ''c'':

Luego verificamos todos los patrones que podrían colocarse en esa ubicación. Ej: si la celda vecina está a la izquierda de la celda actual (lado este), observamos todos los patrones que se pueden colocar a la izquierda de cada patrón contenido en la celda actual.

possible = set([n for idP in W[idC] for n in A[idP][dir]])

También observamos los patrones que están disponibles en la celda vecina:

available = W[idN]

Ahora nos aseguramos de que la celda vecina realmente tenga que actualizarse. Si todos sus patrones disponibles ya están en la lista de todos los patrones posibles -> no hay necesidad de actualizarlo (el algoritmo omite este vecino y pasa al siguiente):

if not available.issubset(possible):

Sin embargo, si no es un subconjunto de la lista possible -> miramos la intersección de los dos conjuntos (todos los patrones que se pueden colocar en esa ubicación y que, "por suerte", están disponibles en esa misma ubicación):

intersection = possible & available

Si no se cruzan (patrones que podrían haberse colocado allí pero no están disponibles) significa que nos encontramos con una "contradicción". Tenemos que detener todo el algoritmo WFC.

if not intersection: print ''contradiction'' noLoop()

Si, por el contrario, se cruzan -> actualizamos la celda vecina con esa lista refinada de índices de patrones:

W[idN] = intersection

Debido a que esa celda vecina se ha actualizado, su entropía también debe actualizarse:

lfreqs = [freqs[i] for i in W[idN]] H[idN] = (log(sum(lfreqs)) - sum(map(lambda x: x * log(x), lfreqs)) / sum(lfreqs)) - random(.001)

Finalmente, y lo más importante, agregamos el índice de esa celda vecina a la pila para que se convierta en la siguiente celda actual por turnos (aquella cuyos vecinos se actualizarán durante el siguiente ciclo while):

stack.add(idN)

Script actualizado completo

from collections import Counter from itertools import chain from random import choice w, h = 40, 25 N = 3 def setup(): size(w*20, h*20, P2D) background(''#FFFFFF'') frameRate(1000) noStroke() global W, A, H, patterns, freqs, npat, mat, xs, ys img = loadImage(''Flowers.png'') iw, ih = img.width, img.height xs, ys = width//w, height//h kernel = [[i + n*iw for i in xrange(N)] for n in xrange(N)] mat = ((-1, 0), (1, 0), (0, -1), (0, 1)) all = [] for y in xrange(ih): for x in xrange(iw): cmat = [[img.pixels[((x+n)%iw)+(((a[0]+iw*y)/iw)%ih)*iw] for n in a] for a in kernel] for r in xrange(4): cmat = zip(*cmat[::-1]) all.append(cmat) all.append(cmat[::-1]) all.append([a[::-1] for a in cmat]) all = [tuple(chain.from_iterable(p)) for p in all] c = Counter(all) patterns = c.keys() freqs = c.values() npat = len(freqs) W = [set(range(npat)) for i in xrange(w*h)] A = [[set() for dir in xrange(len(mat))] for i in xrange(npat)] H = [100 for i in xrange(w*h)] for i1 in xrange(npat): for i2 in xrange(npat): if [n for i, n in enumerate(patterns[i1]) if i%N!=(N-1)] == [n for i, n in enumerate(patterns[i2]) if i%N!=0]: A[i1][0].add(i2) A[i2][1].add(i1) if patterns[i1][:(N*N)-N] == patterns[i2][N:]: A[i1][2].add(i2) A[i2][3].add(i1) def draw(): global H, W emin = int(random(w*h)) if frameCount <= 1 else H.index(min(H)) if H[emin] == ''c'': print ''finished'' noLoop() id = choice([idP for idP in W[emin] for i in xrange(freqs[idP])]) W[emin] = [id] H[emin] = ''c'' stack = set([emin]) while stack: idC = stack.pop() for dir, t in enumerate(mat): x = (idC%w + t[0])%w y = (idC/w + t[1])%h idN = x + y * w if H[idN] != ''c'': possible = set([n for idP in W[idC] for n in A[idP][dir]]) if not W[idN].issubset(possible): intersection = possible & W[idN] if not intersection: print ''contradiction'' noLoop() return W[idN] = intersection lfreqs = [freqs[i] for i in W[idN]] H[idN] = (log(sum(lfreqs)) - sum(map(lambda x: x * log(x), lfreqs)) / sum(lfreqs)) - random(.001) stack.add(idN) fill(patterns[id][0]) rect((emin%w) * xs, (emin/w) * ys, xs, ys)

Mejoras generales

Además de estas correcciones, también realicé una pequeña optimización de código para acelerar los pasos de observación y propagación, y acortar el cálculo de la opción ponderada.

  • La "Ola" ahora está compuesta de conjuntos de índices de Python cuyo tamaño disminuye a medida que las celdas se "colapsan" (reemplazando grandes listas de tamaños fijos de booleanos).

  • Las entropías se almacenan en un defaultdict cuyas claves se eliminan progresivamente.

  • El valor de entropía inicial se reemplaza por un entero aleatorio (el primer cálculo de entropía no es necesario ya que el alto nivel de incertidumbre equiparable al inicio)

  • Las celdas se muestran una vez (evitando almacenarlas en una matriz y redibujando en cada cuadro)

  • La opción ponderada ahora es de una sola línea (evitando varias líneas prescindibles de comprensión de la lista)