Visualizando un wormhole como Interstellar

El equipo de DNEG creó los VFX de Interstellar y también publicaron un paper en el cual describen como imaginaron un wormhole para la cinta. Fueron tan geniales que incluyeron pasos para replicar un visual similar. ¡Agradezco mucho al equipo de DNEG por su trabajo en el filme y por compartir su proceso!

Este post es para organizar un poco los pasos prácticos del paper para ayudar a aquellos que quieren implementar la visualización tambien.

La herramienta que usaremos es Python. Específicamente:

1. PIL para abrir y combinar las imágenes
2. SciPy para calcular las ecuaciones
3. Polars para acomodar los resultados

Puedes instalarlos uno a uno, o si usas Anaconda puedes clonar el código del Github y usar:

conda env create -f environment.yml
conda activate visualizing-wormhole

Luego abre el archivo main.py para seguir los pasos a continuacion.

Imagenes iniciales

Necesitamos dos imagenes para crear la visualización. Una para el cielo como visto a cada entrada del wormhole. Afortunadamente, DNEG proveyó las siguientes en su paper:

Puedes usar estas o las que quieras. Segun te guste como vaya quedando al final!

Ecuaciones

Son unas cuantas que iremos organizando. Además, leerlas en LateX me confundio unas veces (ej. estan multiplicando o dividiendo?) so optaré por presentarlas en el código en vez.

Constantes

Solo hay una, que es el ratio entre W y M,

W_M = -np.log( 1/np.cos(np.pi/(2*np.sqrt(2))) ) + (np.pi/(2*np.sqrt(2)))*np.tan(np.pi/(2*np.sqrt(2)))
# W_M = 1.42953...

Aja, trato de no mirar eso mucho pero ya está set.

Posición de la cámara

Hay que elegir desde donde estamos mirando el wormhole. Todo lo que sigue creará una mirada de 360 alrededor de este punto. Los números para theta y phi son “estables” segun el paper, y el l_c es como el “zoom” al wormhole.

# Elegir posición de la cámara. 'c' significa 'camera'
l_c = -50 * a
theta_c = np.pi / 2
phi_c = np.pi/2

Parametros del wormhole

Estas son las variables que definen el tamaño y comportamiento del wormhole. Usaremos los parametros mencionados en el paper que fueron usados para la película.

# Parametros del wormhole
throat_radius = 1 # km
lensing = 0.05 * throat_radius # km
throat_length = 0.01 * throat_radius # km

# Son los mismos parametros pero con los nombres que usan las ecuaciones
rho = throat_radius # km
W = lensing
a = throat_length / 2 # km

# Constantes del wormhole
M = W / W_M
dr_dl = lambda l: (2 / np.pi) * np.arctan( 2*l / (np.pi * M))

Para cada rayo de luz…

…hay que calcular las ecuaciones. Por esto puede tomar mucho tiempo correr el código. La primera parte cuenta los rayos a calcular y la segunda organiza los loops con coordenadas para cada uno.

# Contabilidad de cuantos rayos se van a calcular
n_phi_samples = out_size[0]
n_theta_samples = out_size[1]
n_total_samples = n_phi_samples * n_theta_samples

i = 0 # Contador de rayos
results = [] # Contenedor de todos los resultados
for (yidx_cs, theta_cs) in enumerate(np.linspace(0, np.pi, n_theta_samples, endpoint=False)):
    for (xidx_cs, phi_cs) in enumerate(np.linspace(0, np.pi*2, n_phi_samples, endpoint=False)):
        print('i / total:', i, '/', n_total_samples)

        # Diccionario de resultados para este rayo de luz
        res = dict(i=i, theta_cs = theta_cs, phi_cs = phi_cs, xidx_cs = xidx_cs, yidx_cs = yidx_cs)

        ... # Todo lo que sigue:

Nota: Muchos de los siguientes pasos ocurren dentro de este loop doble. Para cualquier duda sobre que va donde, consulte el código que está en el orden correcto en el Github.

Condiciones iniciales del rayo

Estas son el punto de partida para las ecuaciones. Es decirle al código donde comienza el rayo de luz para que pueda marchar hacia delante a través del wormhole. Todas las variables cuyo nombre culmina en “0” son valores iniciales que pueden cambiar luego.

# Unit vector pointing in one direction to sky
# 'cs' significa 'celesital sphere', NO camera
N_x = np.sin(theta_cs) * np.cos(phi_cs)
N_y = np.sin(theta_cs) * np.sin(phi_cs)
N_z = np.cos(theta_cs)

# Componentes del rayo en global spherical polar basis
n_l = -N_x
n_phi = -N_y
n_theta = +N_z

# Condiciones iniciales del rayo
l_0 = l_c
theta_0 = theta_c
phi_0 = phi_c

r_0 = get_r(l_0) # basado en l_0

# Canonical momenta del rayo
p_l0 = n_l
p_theta0 = r_0 * n_theta
p_phi0 = r_0 * np.sin(theta_0) * n_phi

# Constantes de movimiento del rayo
b = p_phi0
B = np.sqrt( p_theta0**2 + (p_phi0 / np.sin(theta_0))**2 )

Ecuación del rayo

Esta es LA ecuación que se usa para dar marcha atras en el tiempo. Le pasamos el estado (“state”) del rayo y calculamos el cambio (“dstate”) que va a ocurrir.

# Ecuacion diferencial del rayo
def ray_geodesic(t, state):
    l, theta, phi, p_l, p_theta = state
    r = get_r(l)

    # Derivadas del estado ('state') del rayo en este punto
    dl = p_l
    dtheta = p_theta / r**2
    dphi = b / (r * np.sin(theta))**2
    dp_l = B**2 * (dr_dl(l)) / r**3 # FILL IN dr_dl args
    dp_theta = (b / r)**2 * np.cos(theta) / np.sin(theta)**3

    dstate = [dl, dtheta, dphi, dp_l, dp_theta]

    return dstate

Resolvemos la ecuación pasando los valores iniciales y la ecuacion del rayo al solve_ivp de SciPy:

# Resolver la ecuacion diferencial
y0 = [l_0, theta_0, phi_0, p_l0, p_theta0]
t_span = [0, -10]
result_solve_ivp = solve_ivp(ray_geodesic, t_span, y0,
                            # method='LSODA',
                            # max_step=0.01
                            )
# print(result_solve_ivp.message)
# print(result_solve_ivp)

Y guardamos resultados y cerramos el doble loop!

# Extraer los tiempos marcados y los estados del rayo
ts = result_solve_ivp.t
ys = result_solve_ivp.y

tlast = ts[-1]

# Guardar resultados
res['tlast'] = tlast
res['l_last'] = ys[0, -1]
res['theta_last'] = ys[1, -1]
res['phi_last'] = ys[2, -1]

results.append(res)
i += 1

Resultados

Usamos Polars para manejar como DataFrame los distintos resultados de todos los rayos. Hay ajustes que realizar para que sean mejor compatibles para crear imagenes:

# Procesar conjunto de resultados
pl.Config.set_tbl_cols(20)
results = (
    pl.DataFrame(results)
      .with_columns([

        # Arreglar angulos para caer entre 0 y 2pi
        (pl.col('theta_last').apply(lambda x: x % (np.pi*2))).alias('theta_last_mod'),
        (pl.col('phi_last').apply(lambda x: x % (np.pi*2))).alias('phi_last_mod'),

      ])
      .with_columns([

        # Asignar hemisferio
        pl.when(pl.col('l_last') > 0)
            .then( pl.lit('upper') )
            .otherwise( pl.lit('lower') )
            .alias('hemisphere_last'),

        # Arreglar angulos theta para caer entre 0 y pi
        pl.when(pl.col('theta_last_mod') > np.pi)
            .then( 2*np.pi - pl.col('theta_last_mod') )
            .otherwise( pl.col('theta_last_mod') )
            .alias('theta_last_fixed'),

        # Arreglar angulos phi para esos puntos el cual se arreglo theta (lo anterior)
        pl.when(pl.col('theta_last_mod') > np.pi)
            .then( (pl.col('phi_last_mod') + np.pi).apply(lambda x: x % (np.pi*2)) )
            .otherwise( pl.col('phi_last_mod') )
            .alias('phi_last_fixed'),

      ])
      .with_columns([

        # Normalizar para poder compararlo con las imagenes
        (pl.col('theta_last_fixed') / np.pi).alias('theta_pct'),
        (pl.col('phi_last_fixed') / (np.pi*2)).alias('phi_pct'),


      ])
      .with_columns([

        # Buscar coordenadas correspondiente en la imagen
        # para phi
        pl.when(pl.col('hemisphere_last') == 'lower')
            .then( pl.col('phi_pct') * lower_im.size[0] )
            .otherwise( pl.col('phi_pct') * upper_im.size[0] )
            .floor()
            .cast(int)
            .alias('phi_im_x'),

        # para theta
        pl.when(pl.col('hemisphere_last') == 'lower')
            .then( pl.col('theta_pct') * lower_im.size[1] )
            .otherwise( pl.col('theta_pct') * upper_im.size[1] )
            .floor()
            .cast(int)
            .alias('theta_im_y'),
      ])
)

La idea es asignar si el rayo de luz culmino en el hemisferio superior o inferior. Con esa información sabemos cual de las dos imágenes hay que mirar.

Luego traducimos esas coordenadas finales (phi, theta) a coordenadas de la imagen (x, y).

Construir la imagen

Cada resultado tiene un pixel en una imagen que vamos copy-pasting a la imagen final

# Pintar resultados sobre la imagen nueva
for row in results.iter_rows(named=True):
    i = row['xidx_cs']
    j = row['yidx_cs']
    if row['hemisphere_last'] == 'lower':
        im = lower_im
    else:
        im = upper_im
    i_im = min(row['phi_im_x'], im.size[0]-1)
    j_im = min(row['theta_im_y'], im.size[1]-1)
    im_pixel = im.getpixel((i_im, j_im))
    # print('placing im_pixel', im_pixel, 'at', (i,j))
    out_im.putpixel((i, j), im_pixel)

Y conviene rodar la imagen final para que el wormhole quede en el centro:

roll_delta = out_im.size[0]//2
rolled_out = roll_image(out_im, roll_type='horizontal', delta=roll_delta)

Ya está!

rolled_out.show() # Enseñar la imagen final!

Recuerda guardar la imagen a un archivo antes de cerrar la ventanilla en que aparece!

James Webb wormhole

So recuerda que puedes usar las imagenes que quieras para los cielos. Por ejemplo, aquí tomamos dos imagenes capturadas por el James Webb Space Telescope:

1. Stephan’s Quintet
2. Southern Ring Nebula

Las editamos y rodamos un poco y asi queda:

Final

El código está disponible en Github. Para experimentar, puedes cambiar las imágenes, croppearlas y rodarlas, aumentar la resolución de la imagen final entre otras. Si estan tardando demasiado los renders, reduce la resolución de la imagen final o de las imagenes iniciales.