Ruido Perlin
Vamos con una explicación (seguida por la implementación de rigor) del algoritmo de ruido más usado en todo el mundo; el ruido Perlin.
El modelo conceptual del ruido Perlin fue descubierto en 1982 por Ken Perlin como resultado de su trabajo de generación de texturas para la película Tron. Desde entonces, el ruido Perlin, y pequeñas variantes posteriores (como el ruido Simplex) son la base para generar efectos como fuego, humo, nubes… texturas como madera, mármol, granito… mapas de altura, etc, de forma totalmente automática. En esencia, nos permite generar patrones de forma muy similar a como lo suele hacer la naturaleza.
En esta entrada vamos a explicar todo el proceso, la naturaleza del ruido, a implementarlo, y finalmente a usarlo para generar terrenos en 3D como el de la siguiente imagen (y un visor interactivo para probarlos, al final de la entrada):
RUIDO BLANCO
Primero, entendamos las características del ruido aleatorio para comprender en qué se diferencia el ruido Perlin de él. Si generamos un mapa 2D de ruido blanco clásico (valores (pseudo)aleatorios), obtenemos algo similar a la siguiente imagen:
Exacto, es muy parecido a lo que veíamos en nuestras viejas televisiones cuando no habíamos seteado ningún canal; sólo que entonces era el resultado del ruido blanco electromagnético del ambiente, completamente caótico, y en este caso sólo son bits al azar.
¿Qué características importantes vemos en esa imagen?
1) Desde el punto de vista de la teoría de la información, la entropía es máxima. Así que esta imagen no puede comprimirse al codificarse sólo con los símbolos usados en su lenguaje (blanco/negro).
2) Desde un punto de vista físico, su densidad media es prácticamente constante (tiende fuertemente a serlo) para todos los fragmentos de dimensiones superiores a 1×1. Es decir, si hacemos un desenfoque gaussiano la imagen se volverá esencialmente plana, incluso con radios de tan sólo 2px.
RUIDO PERLIN (1 dimensión)
En el ruido Perlin las características son las opuestas, que es lo deseable cuando trabajamos en generación procedural o automática. Si entendemos cómo se crea el ruido, es fácil y lógico entender el por qué de esta naturaleza.
Así que empecemos:
Vamos a calcular el ruido perlin como la generación de onda n-dimensional compuesta por k octavas. Si pensamos en la función seno para representar una onda:
Vemos como su amplitud va desde -1 a 1 (2 en total) y su longitud de onda es 6. Esta bien podría ser la primera componente de nuestra onda final, en el siguiente paso, queremos una onda de la misma naturaleza, pero con mayor frecuencia (menor longitud de onda) y menor amplitud. Por ejemplo:
La tercera y cualta deberían obedecer al mismo principio (suponiendo que sólo usáramos 4):
De esta forma, si sumamos las componentes anteriores, y en general cualesquiera en las que cada una aumenta la frecuencia y reduce la amplitud de forma progresiva, obtendremos algo como lo siguiente:
Esta onda unidimensional, que se asemeja al perfil de un paisaje montañoso cumple el siguiente principio, que tanto observamos en la naturaleza: A mayor escala, menor número de variaciones pero con más influencia, a menor escala, mayor número de variaciones pero con menor influencia. O dicho de otra forma, a altas frecuencias bajas amplitudes y viceversa.
Pero en la naturaleza los procesos no están formados a base de ondas periódicas como el seno que hemos usado hasta ahora. Debemos crear una onda pseudo-aleatoria simplemente eligiendo valores aleatorios en los “nodos”. Los nodos son los puntos en los que una onda periódica alcanzaría un máximo o mínimo, es decir, debemos introducir valores cada λ unidades de espacio (longitud de onda).
Ejemplo:
puntos generadores aleatoriamente para generar una onda
Una vez generados los puntos, debemos interpolar entre los valores de los nodos para rellenar todo el dominio en el que se define la onda. Podemos usar aquí cualquier tipo de interpolación: linear, trigonométrica, Lagrange, B-splines… pero usualmente con una simple interpolación linear es suficiente. Los resultados finales no mejoran demasiado con interpolaciones más complejas (como los B-Splines, que sería la ideal) y que tienen un coste computacional muy alto. No obstante sois libres de usar la que mejor os parezca.
Nota: Tengo pendiente una entrada sobre interpolación, usando B-Splines (splines cúbicos) y Lagrange, que sí son muy útiles en otros cientos de aplicaciones.
interpolación lineal
interpolación por cosenos
interpolación cúbica
Para utilizar una nomenclatura adecuada, a cada componente vamos a llamarla octava, y usaremos un caso típico en el que en cada octava la frecuencia se multiplica por dos y la amplitud si divide por dos (de hecho, es así cómo se definen las octavas en música).
Con todo esto, ya podemos empezar a generar nuestra primera onda de ruido Perlin real. Empezaremos con una frecuencia de 2 (dos cambios en todo el dominio) y una amplitud de 128 (para representarla en una imagen de 256*256 siendo el centro del eje de ordenadas el centro de la imagen.
frecuencia: 2
amplitud: 128
frecuencia: 4
amplitud: 64
frecuencia: 8
amplitud: 32
frecuencia: 16
amplitud: 16
frecuencia: 32
amplitud: 6
frecuencia: 64
amplitud: 4
frecuencia: 128
amplitud: 2
El resultado de sumar estas 7 octavas es el siguiente:
Ruido Perlin final con 7 octavas
Y el código que he escrito para generarlas es el siguiente:
//=============
//Interpolation
//=============
var interpolate={}
interpolate.linear=function(x1,y1,x2,y2,x)
{
return ((x-x1)*(y2-y1)/(x2-x1))+y1;
}
//=============
//Perlin Noise
//=============
var perlin={}
perlin.octave(frequency, amplitude, length)
{
var wave=new Array(length);
var step_length=Math.floor(length/frequency);
for(var i=0;i<=length;i+=step_length)
{
//Create the nodes of the wave
wave[i]=(amplitude/2)-(Math.random()*amplitude);
//Interpolation between nodes
if(i>=step_length)
{
for(var j=i-(step_length-1);j<i;j++)
{
var x1=i-(step_length);
var x2=i;
wave[j]=interpolate.linear(x1,wave[x1],x2,wave[x2],j);
}
}
}
return wave;
}
Como veis el código es muy sencillo. El objeto “interpolate” está preparado para expandirse añadiendo nuevos métodos de interpolación si fuese necesario. Para generar todas las octavas y sumarlas, usando una progresión de doble frecuencia/mitad amplitud, podríamos simplemente hacer algo ad hoc como lo siguiente:
var lenght=256;
var perlin_wave=new Array(lenght);
var n_octaves=7;
for(var i=0;i<perlin_wave.length;i++) //wave initialization
perlin_wave[i]=0;
for(var i=1;i<=n_octaves;i++)
{
var current_pow=Math.pow(2,i);
var octave=perlin.octave(current_pow,lenght/current_pow,lenght);
for(var j=0;j<perlin_wave.length;j++)
perlin_wave[j]+=octave[j];
}
O bien, podemos crear un método dentro del objeto perlin con el código anterior si lo que queremos es tenerlo encapsulado.
Con esto ya tenemos nuestro generador de ruido Perlín en una dimensión. Aunque lo interesante viene para el caso en 2D.
Prueba a generar nuevas ondas con nuestra implementación haciendo clic en el botón
Generar una nueva onda
RUIDO PERLIN (2 dimensiones)
Para el caso de 2 dimensiones debemos pensar en la onda y sus octavas como mapas o texturas. Una octava 2D con frecuencia 2 creará dos valores aleatorios por dimensión, es decir, un total de 2×2=4 valores que estarán situados en las 4 esquinas del mapa. Si la amplitud es de 256 con valores entre -128 a 128 podemos normalizarla a valores entre 0 y 256 y dibujar en escala de grises (un gris para cada valor) la onda de forma muy gráfica y representativa.
Ejemplo:
En este caso, hemos rellenado el dominio vacío de la onda, con el color del punto definido que manda sobre el cuadrante; dado que como hemos dicho, sólo los vértices del mapa tienen un valor establecido. Sin embargo, lo ideal es de nuevo usar algún método de interpolación para rellenar todo el dominio. De nuevo el más usado para este caso es la interpolación lineal aplicada a dos dimensiones (bilinear), que está perfectamente explicada en la wikipedia. De esta forma, obtenemos algo como lo siguiente:
Usando estos principios, podemos empezar a generar octavas y crear nuestro ruido Perlin en 2D. Por ejemplo:
frecuencia: 2
amplitud: 128
frecuencia: 4
amplitud: 64
frecuencia: 8
amplitud: 32
frecuencia: 16
amplitud: 16
frecuencia: 32
amplitud: 8
frecuencia: 64
amplitud: 4
frecuencia: 128
amplitud: 2
La suma de todas estas octavas da como resultado:
Ruido Perlin 2D final con 7 octavas
Esta imagen ya tiene semejanza con algo que parecen ser nubes o humo sobre un fondo, de hecho, si dejamos el canal B a 255 y el canal G como el máximo de (185,valor_de_la_onda), habremos creado un sencillo filtro que produce imágenes como la siguiente:
nubes fluorescentes lanshóricas
El código para el ruido Perlin en 2D que he escrito para generar los ejemplos anteriores es el siguiente:
//=============
//Interpolation
//=============
var interpolate={}
interpolate.linear=function(x1,y1,x2,y2,x)
{
return ((x-x1)*(y2-y1)/(x2-x1))+y1;
}
interpolate.bilinear=function(x1,y1,x2,y2,v1,v2,v3,v4,tx,ty)
{
//P1:{x1,y1,v1} - P2:{x2,y1,v2} - P3:{x1,y2,v3} - P4:{x2,y2,v4}
//Target:{tx,ty}
//NOTE: Bind each area with the oposite value
var area_v1=Math.abs((tx-x1)*(ty-y1))*v4;
var area_v2=Math.abs((tx-x2)*(ty-y1))*v3;
var area_v3=Math.abs((tx-x1)*(ty-y2))*v2;
var area_v4=Math.abs((tx-x2)*(ty-y2))*v1;
var area_total=(x2-x1)*(y2-y1);
return (area_v1+area_v2+area_v3+area_v4)/area_total;
}
//=============
var perlin={}
perlin.octave_2d=function(frequency, amplitude, length)
{
var step_length=Math.floor(length/frequency);
var wave=new Array(length+1);
for(var i=0;i<=length;i++)
wave[i]=new Array(length);
for(var i=0;i<=length;i+=step_length)
{
for(var j=0;j<=length;j+=step_length)
{
//Create the nodes of the wave
wave[i][j]=(amplitude/2)-(Math.random()*amplitude);
//Interpolation between nodes
if(i>=step_length && j>=step_length)
{
for(var a=i-(step_length);a<=i;a++)
{
for(var b=j-(step_length);b<=j;b++)
{
var x1=i-(step_length);
var x2=i;
var y1=j-(step_length);
var y2=j;
wave[a][b]=interpolate.bilinear(x1,y1,x2,y2,wave[x1][y1],wave[x2][y1],wave[x1][y2],wave[x2][y2],a,b);
}
}
}
}
}
return wave;
}
//=============
function wave2D(data)
{
this.map=data;
this.add=function(data)
{
for(var i=0;i<data.length;i++)
{
for(var j=0;j<data.length;j++)
this.map[i][j]+=data[i][j];
}
}
this.render=function(ctx)
{
for(var i=0;i<this.map.length;i++)
{
for(var j=0;j<this.map.length;j++)
{
//Normalize from 0 to 255
var depth=Math.floor(128+this.map[i][j]);
//LaNsHoRic clouds: ("+depth+","+Math.max(185,depth)+",255)
ctx.fillStyle="rgb("+depth+","+depth+","+depth+")";
ctx.fillRect(i,j,1,1);
}
}
}
this.toURL=function()
{
var canvas=new Canvas();
canvas.width=canvas.height=this.map.length;
var ctx=canvas.getContext("2d");
this.render(ctx);
return canvas.toDataURL();
}
}
Con esto, podemos adaptar el código de generación anterior a lo siguiente:
var final_wave=new wave2D(perlin.octave_2d(1,256,256));
for(var i=1;i<8;i++) //7 octaves
{
var current_pow=Math.pow(2,i);
var octave=perlin.octave_2d(current_pow,256/current_pow,256);
final_wave.add(octave);
}
Prueba a generar nuevos mapas con nuestra implementación haciendo clic en el botón
Generar un nuevo mapa
Podemos cambiar la forma en que aumenta la frecuencia y disminuye la amplitud para obtener diferentes resultados. Jugando con la progresión y los valores podemos obtener texturas como granito, nubes dispersas, humo denso, etc. O usar el ruido como perturbación de un mapa base, con esto podemos generar lava, madera, mármol, relámpagos, etc. Sólo es cuestión de ir jugando con valores y experimentar un poco.
El siguiente uso popular de los mapas de ruido Perlin como el que acabamos de generar, es como mapas de altura. De esta forma podemos generar terrenos de forma automática, que en base a nuestros parámetros de generación podrán ser más o menos montañosos, más escarpados, más suaves…
En el siguiente cuadro puedes ver como usando nuestro generador de ruido Perlin creamos terrenos en 3D. Puedes hacer clic con el ratón y arrastrar para rotar la cámara. La implementación del modelado en 3D y del visor del mapa de altura están fuera del contenido de la entrada, pero prometo preparar una hablando de ello si os interesa (notificádmelo).
Como resumen rápido, dividimos un plano en tantas secciones de como frecuencia f tiene la mayor octaba. El plano queda con f*f vértices que asociamos con los píxeles del mapa 2D generado, para elevarlos según el valor del mapa en ese punto. Así el plano queda deformado como en los siguientes ejemplos:
Si os sirve de ayuda no es necesario mención pero sí sería cortés un agradecimiento aunque fuese en un comentario 🙂
Leer másAtractores, caos, sistemas complejos y mariposas
Siempre nos sentimos más cómodos con aquello que creemos que conocemos; con aquello que ya sabemos y no puede cambiar. Pero en en realidad, siendo un poco puristas, la estabilidad no existe. Ningún sistema en la naturaleza (sistemas reales, no artefactos matemáticos imaginarios) es estable; en el sentido de que no va a permanecer inalterado para siempre. No existe nada en este universo, que haya detenido su evolución nunca.
En la práctica, cuando hacemos ciencia y estudiamos sistemas “aislados”, rebajamos considerablemente nuestros requisitos de estabilidad; y decimos que un sistema ha alcanzado la estabilidad cuando su evolución cesa durante “un tiempo lo suficientemente largo”.
Sí, es tan poco riguroso como suena…
Pero lo divertido está en la incertidumbre, y una vez perdemos el miedo a lo desconocido podemos empezar a disfrutar de ella. Veamos algunos tipos de sistemas con ejemplos sencillos y clasifiquémoslos:
Sistemas Simples
Un ejemplo de sistema simple es una taza de café caliente suspendida en una atmósfera de aire. Poco a poco la taza irá enfriándose hasta alcanzar la temperatura ambiente del aire (el sistema alcanzará la estabilidad). En sus inicios, el sistema no es estable, pero nos sentimos casi tan cómodos con estos sistemas como con los estables porque es fácil predecir su evolución.
Sistemas Complejos
Un sistema complejo es un sistema cuya evolución es muy difícil de determinar; ya que existen factores no evidentes que afectan de forma dramática a cómo evoluciona el sistema. Un ejemplo es el campo gravitatorio, cuando hay varios cuerpos en juego y en movimiento; es difícil calcular la evolución de las trayectorias.
Para verlo mejor; he escrito un pequeño script para hacer simulaciones (del que hablaremos después). En la primera simulación tendréis 6 puntos negros que representan cuerpos “fijos” que atraen a otro “móvil” cuya trayectoria queda marcada por una estela azul. La simulación es infinita y no se estabiliza nunca (aunque durante un rato el objeto móvil se salga de la pantalla, tarde o temprano volverá).
Para empezar la simulación haced clic en “empezar simulación” (podéis reiniciarla cuando queráis). ¿Podéis calcular mentalmente su trayectoria (la evolución del sistema)? Entonces, es un sistema complejo. Otro ejemplo de sistema complejo es el clima de la tierra.
Sistemas Caóticos
Un sistema caótico es un sistema, que además de complejo, cumple la siguiente propiedad fundamental: un mínimo y casi insignificante cambio en las condiciones iniciales hace que el sistema evolucione de forma radicalmente distinta.
Para ilustrarlo vamos a hacer una simulación como la anterior; pero esta vez, habrá tres objetos móviles (uno rojo, uno verde y uno azul). Los tres objetos estarán en la misma posición, pero con una diferencia de 0.001 píxeles. Esta diferencia, provocará que con el tiempo, los tres objetos sigan trayectorias completamente distintas.
comenzar simulación
La vida es un propio ejemplo de sistema caótico. A veces pensamos que ya hemos vivido situaciones similares a las que estamos viviendo en un momento dado; y creemos que la experiencia pasada nos puede servir, pensando que ya sabemos qué va a ocurrir. Pero en realidad, cualquier mínimo detalle puede ser determinante para que la evolución de sucesos sea completamente distinta.
De los sistemas caóticos como este viene “la teoría del caos”; y la famosa metáfora de la mariposa que bate sus alas y cambia radicalmente los acontecimientos futuros en una suerte de efecto dominó cuyas consecuencias se van amplificando. Es una forma de decir; que el más mínimo detalle es crucial para la evolución del sistema.
NOTA: Los sistemas caóticos, a pesar de su nombre, no son azarosos en absoluto; son completamente deterministas.
Atractores
Los atractores, como mención importante (fuera de los propósitos de esta entrada), son sistemas complejos y o caóticos, pero, en los que las trayectorias próximas permanecen siempre próximas, siento permisivos con ciertas perturbaciones en sus condiciones iniciales y evolución.
A todos os sonará el más famoso de todos: el atractor de Lorenz (prometo dedicar una entrada para hablar exclusivamente de atractores ;))
El código
Propongo como ejercicio escribir un sistema de simulación como el que hemos usado en esta entrada. Lo único un poco más complicado (sigue siendo sencillo) es hacer las estelas que se van difumando.
Si queréis la solución directa, mi código, con una inicialización de ejemplo que ejecuta una simulación con 25 atractores y 150 objetos móviles (todos aleatorios) es:
function Attractor(canvas)
{
//-[utils]---------------- *********************************************
function randChannel() {return Math.floor(Math.random()*255)}
function randPosition() {return [Math.random()*canvas.width,Math.random()*canvas.height]}
//-[balls]---------------- *********************************************
function Ball(x, y, mass, rgb)
{
//-{init}
this.x=x || randPosition()[0];
this.y=y || randPosition()[1];
this.mass=mass || 4;
this.vx=0;
this.vy=0;
this.radius=this.mass/2;
if(rgb)
this.color="rgb("+rgb[0]+","+rgb[1]+","+rgb[2]+")";
else
this.color="rgb("+randChannel()+","+randChannel()+","+randChannel()+")";
balls.push(this);
//-{methods}
this.move=function()
{
var tmp_vx=0;
var tmp_vy=0;
for(var i=0;i<acelerators.length;i++)
{
var dx=acelerators[i].x-this.x;
var dy=acelerators[i].y-this.y;
var theta=Math.atan2(dy,dx);
var mod=Math.sqrt(Math.pow(dx,2)+Math.pow(dy,2));
var force=Math.min(acelerators[i].force/Math.pow(mod,2),3);
tmp_vx+=force*Math.cos(theta);
tmp_vy+=force*Math.sin(theta);
}
this.vx+=tmp_vx/this.mass;
this.vy+=tmp_vy/this.mass;
this.x+=this.vx;
this.y+=this.vy;
}
this.render=function()
{
current_ctx.save();
current_ctx.fillStyle=this.color;
current_ctx.beginPath();
current_ctx.arc(this.x, this.y, this.radius, 0, 2*Math.PI, false);
current_ctx.closePath();
current_ctx.fill();
current_ctx.restore();
}
}
//-[acelerator]---------------- *********************************************
function Accelerator(x, y, force)
{
//-{init}
this.x=x || randPosition()[0];
this.y=y || randPosition()[1];
this.force=force || 200;
this.radius=this.force/50;
acelerators.push(this);
//-{methods}
this.render=function()
{
ctx.save();
ctx.fillStyle="#000";
ctx.beginPath();
ctx.arc(this.x, this.y, this.radius, 0, 2*Math.PI, false);
ctx.closePath();
ctx.fill();
ctx.restore();
}
}
//-[attractor]----------- *********************************************
//-[setup]----------------
var ctx=canvas.getContext("2d"); //main context
var balls=[]; //balls array
var acelerators=[]; //acelerators array
//-[wake]-----------------
var wake_resolution=20; //number of canvas for the wake
var ctx_buffer=[]; //context array for the wake
var current_ctx=0; //current context object
var arrow=0; //pointer to the current context
var wake_time=2000/wake_resolution; //5 seconds of wake
var last_timestamp=0; //last time of buffer change
var frames=0; //frame count
//-[init]-----------------
for(var i=0;i<wake_resolution;i++)
{
var tmp_canvas=document.createElement("canvas");
tmp_canvas.width=canvas.width;
tmp_canvas.height=canvas.height;
ctx_buffer[i]=tmp_canvas.getContext("2d");
}
current_ctx=ctx_buffer[0];
//-[run]------------------
this.run=function()
{
frames++;
ctx.clearRect(0,0,canvas.width,canvas.height);
//chance of change current context
var now=Date.now();
if(now-last_timestamp>wake_time)
{
last_timestamp=now;
current_ctx=ctx_buffer[(++arrow%ctx_buffer.length)];
current_ctx.clearRect(0,0,canvas.width,canvas.height);
}
//render balls at current context
for(var i=0;i<balls.length;i++)
{
balls[i].move();
balls[i].render();
}
//composition of the wake at main context
var step_alpha=1/ctx_buffer.length;
for(var i=1;i<=ctx_buffer.length;i++)
{
var tmp_ctx=ctx_buffer[(arrow+i)%ctx_buffer.length];
ctx.save();
ctx.globalAlpha=i*step_alpha;
ctx.drawImage(tmp_ctx.canvas,0,0);
ctx.restore();
}
//render acelerators at main context
for(var i=0;i<acelerators.length;i++)
acelerators[i].render();
}
//-[stop]-----------------
this.stop=function()
{
clearInterval(this.id_int); //external setter
}
//-[classes]-----------------
this.Ball=Ball;
this.Accelerator=Accelerator;
}
//gogogo!
(function gogogo()
{
var attractor=new Attractor(canvas);
for(var i=0;i<25;i++)
new attractor.Accelerator();
for(var i=0;i<150;i++)
new attractor.Ball();
attractor.id_int=setInterval(function(){attractor.run()},60/1000); //run at 60 fps
}());
Más simulaciones
Sistema aleatorio con 10 atractores y 20 objetos móviles:
reiniciar simulación aleatoria
Dedicatoria
Esta entrada está dedicada a Virgy; que me inspiró para escribir el simulador tras una conversación; y que finalmente he aprovechado (al margen de su fin original) para hacer esta entrada 🙂
Leer másLa visión en los perros
Llevo tiempo preguntándome cómo se ve el mundo desde los ojos de Arthur. Algunas personas piensan que los perros ven en blanco y negro, pero lo cierto es que sí distinguen algunos colores. En concreto, perros, lobos, zorros… y la mayoría de los cánidos y mamíferos poseen una forma de visión dicromática llamada deuteranopía; y es también, de hecho, una forma de daltonismo en humanos.
Los colores
Para entender nuestra visión: en la retina humana (y en la de muchos primates) podemos encontrar bastones (células de visión por intensidad) y tres tipos de conos (células responsables de la visión en color); cada tipo de cono especializado en una longitud de onda diferente: “azules” sensibles en un rango de longitudes centradas en el azul (unos 430nm), “verdes” especializados en longitudes de onda similares a la del verde (unos 520nm) y “rojos” para longitudes de onda largas cercanas al rojo (~650nm).
Dado que estas son las longitudes de los colores primarios aditivos (rojo – verde- azul); la composición de la información receptada por cada tipo de cono completa un espectro como el siguiente:
En los animales con deuteranopía los conos responsables de la recepción del verde no están presentes o no son funcionales. La composición de la información receptada presenta en estos casos un espectro como el siguiente:
Niveles de gris
Otra diferencia importante entre la visión humana y la de un perro es la diferencia mínima perceptible entre factores de gris (AR). Para esta medición se suele usar la ley de Weber-Fechner; una regla psicofísica que establece una relación cuantitativa entre la magnitud de un estímulo físico y cómo éste es percibido. El resultado de la fracción de Webber para humanos se estima en 0.11, y para perros en 0.22; esto significa que los perros son sólo capaces de captar la mitad de niveles de grises que nosotros.
Por ejemplo, los perros verían los siguientes cuadrados del mismo color:
Si tú ves el mismo color en ambos cuadrados, puede deberse a que tu monitor sea de gama baja y no sea capaz de mostrar correctamente los colores (o que esté mal configurado en parámetros de brillo, contraste y temperatura). Descartamos la opción de que seas un perro, claro.
Agudeza visual
Dicho en pocas (pero muy imprecisas) palabras; la agudeza visual es la resolución espacial del sistema visual. Al tratarse de un término de resolución es fácil imaginar que está directamente relacionado con la densidad de conos y bastones de la retina. Esta densidad varía según la zona; así, perdemos resolución en los límites de nuestro campo visual (donde apenas podemos percibir objetos) y alcanzamos el máximo en el punto donde centramos la vista.
El máximo de agudeza visual para humanos está entre 50 y 60 CPD. En los perros se sitúa entre 6.5-9 y 11.6 CPD.
Simular cómo es la visión con una agudeza visual diferente a la nuestra es complicado; intervienen muchos factores y es alterada circunstancialmente por el entorno. Por ejemplo el desenfoque por profundidad de campo se dispara (no podemos simularlo en una imagen 2D) al enfocar; y la cantidad de luz recibida también afecta de forma virtual; ya que los perros tienen mayor densidad de bastones que los humanos, pupilas más grandes y un tapetum lucidum refractivo; así en la oscuridad la agudeza visual de los perros acabaría superando a la del humano.
Aunque esta aproximación sea muy imprecisa; podemos ilustrar la pérdida de resolución con la siguiente figura; asumiendo un factor de /4 para los perros.
AV /1
AV /4
Al lío. Dedos al código.
Con todo este resumen; la idea es escribir una pequeña función javascript que dada una imagen de ejemplo la transforme y haga una aproximación acertada de cómo la percibiría un perro.
La función tiene 3 tareas independiente a realizar:
1) Simulación de deuterapía: No es tan sencillo como pueda parecer, ya que los niveles de luz y ponderaciones entre colores cambian gradualmente según la curva de sensibilidad de los conos que sí están presentes (rojo-azul). En este paper [1] encontramos los factores medidos experimentalmente -para humanos; aunque la diferencia con perros es mínima (comparamos con [2])- que necesitaremos para la conversión.
2) Discretizamos algunos niveles de luz. No vamos a implementar todas las curvas; sólo una aproximación. Algo de información adicional en [3]. Hacemos una reducción de la intensidad a niveles pares.
3) Desenfoque gausiano a la imagen. Para hacerlo correctamente el desenfoque variaría su grado arbitrariamente según la zona; y para ello necesitaríamos una imagen tridimensional, información sobre la exposición y el punto de enfoque marcado. Como no tenemos nada de eso, hacemos un desenfoque suave y global que simule el máximo grado de agudeza visual en todos los puntos. El desenfoque podemos hacerlo eficientemente (usando la GPU) mediante un filtro en CSS3 (fuera del código javascript).
function DogVision() //v0.1 - LaNsHoR @ 2014
{
//===================================================
//Factors
var gamma=1.0 /* normal value: 1.8;
-0.8 correction for dogs */
var by=127, k=[9591, 23173, -730];
for(var luts={}, invluts={}, i=0; i<256; i++)
{
luts[i]=0.992052*Math.pow(i/255, gamma)+0.003974;
invluts[i]=255*Math.pow(i/255, 1/1.8);
}
//===================================================
//Functions
//-- RGB normalization
function normalize(value)
{return value<0?0:value>255?255:value;}
//-- Dog light reduction
function even(value) {return value%2?value-1:value;}
//===================================================
this.dogPixel=function(r,g,b)
{
var r_lin=luts[r], g_lin=luts[g], b_lin=luts[b];
var r_blind=((k[0]*r_lin)+(k[1]*g_lin))/by;
var b_blind=((k[2]*r_lin)-(k[2]*g_lin)+32768*b_lin)/by;
r_blind=normalize(r_blind);
b_blind=normalize(b_blind);
var red=invluts[Math.round(r_blind)];
var blue=blue=invluts[Math.round(b_blind)];
return [even(red), even(red), even(blue)];
}
}
Y ya está! Podemos empezar a jugar con los resultados 🙂
Ejemplos
El resultado aparecerá en el rectángulo de borde azul.
Conclusión de andar por casa
Por todo lo dicho y visto; los perros en ambientes de luz diurna tienen una visión mucho más precaria que la humana. Sin embargo, con condiciones de baja luz su percepción de siluetas es mucho mejor que la nuestra (fijaos en los ejemplos). Además, su mayor número de bastones les hacen mejores detectando pequeños movimientos a grandes distancias (algo que no podemos percibir en imágenes estáticas).
No está nada mal. Después de todo; su sentido principal es el olfato (la mayoría de neuronas de su corteza cerebral están dedicadas a este sentido; la mayoría de las nuestras, a la visión). Aún con una visión mediocre, si sumamos sus excelentes olfato y oído su nivel de percepción sensorial global es superior al nuestro. Y desde luego, perfectamente adaptado a lo que necesitan.
Dedicatoria
Esta entrada y el código fuente de la función que hemos escrito están dedicados a Quarkyta… una de las preciosas perritas de Virgy, que nos dejó hace unos días 🙁
Referencias
[1] Digital Video Colourmaps for Checking the Legibility of Displays by Dichromats
[3] The spectral luminosity curves for a dichroma tic eye and a normal eye
Leer másIntersección Rayo-Triángulo
Hagamos un poco de abstract:
La intersección rayo-triángulo en R3 es un problema clásico que a menudo nos encontramos cuando buscamos soluciones a un sistema que depende de la topología de una malla triangular que interactúa con otra entidad en la escena. Especialmente en aplicaciones y variaciones de raytracing tipo mapeado fotónico; esto es, para renders offline.
frame de Big Buck Bunny, renderizado con Blender
Hay muchas técnicas para calcularla, una casi omnipresente y eternamente recomendada es la Intersección Rayo-Triángulo Rápida de Mínimo Almacenamiento de Moller y Trumbore (conocida como MT97); basada en el cambio de coordenadas baricéntricas del triángulo en cuestión.
El paper del último link de MT97 explica bien el algoritmo y proporciona una implementación. Probablemente, es la que necesitáis para vuestro problema.
¿Probablemente?
Lo mejor de MT97 es que nos da una valiosa lección de filosofía. Su defecto, como el de tantas cosas en la vida, es el de su adopción masiva sin planteamiento.
Más que en ningún otro sitio, en este tipo de problemas hay que ser especialmente crítico; sobretodo en las líneas clave que van a suponer más del 70% del código de ejecución (no fuente) del programa.
Así, hace algunos años cuando me encontré con la necesidad de implementar una solución para calcular la intersección rayo/triángulo hice una comparativa (cosas de la intuición) entre MT97 y mi opción alternativa. Sorprendentemente encontré que para mi problema era mucho más rápida (computacionalmente) esta última.
La primera ventaja que intuí antes de hacer la prueba venía por su naturaleza paralelizable. Programándola usando instrucciones SSE y registros XMM* en gran parte del proceso obtenía hasta tres operaciones por instrucción.
El algoritmo es el siguiente:
- Calcular la intersección del rayo con el plano que define el triángulo (problema mucho más sencillo y bien conocido -se da en matemáticas de bachillerato <o al menos yo lo dí>-).
- Sumar las áreas de los tres triángulos formados al unir los vértices del triángulo original con el punto de intersección (si el rayo es aproximado o totalmente paralelo al plano, terminamos aquí).
- Sumar las áreas de los tres triángulos formados y calcular el área del triángulo original.
Se puede demostrar matemáticamente que si la suma de las tres áreas de los formados es igual a la del original, entonces el punto de intersección está dentro del original. Dando como trivial el punto 1, nuestro problema se reduce pues a calcular el área de un triángulo en R3.
La función implementada y con comentarios en español quedaría:
/*
Los registros XMM* son de 128bits, utilizo 4 floats
para llenarlos en una estructura V4.
*/
struct V4
{
float a;float b;float c;float d;
V4(float x=0){a=b=c=d=x;}
V4(float x, float y, float z, float i){a=x;b=y;c=z;d=i;}
};
/* P*: Vértice *(x,y,z) - R: Resultado */
void inline AreaOfATriangle(V4& P1, V4& P2, V4& P3, V4& R)
{
V4 DIV(0.5);
//Nota: Las coordenadas .d deben ser 0
//TODO: Setear a 0 estos últimos 32bits para P1, P2 y P3
asm volatile
(
"MOVUPS %0, %%XMM1\n\t" //Sube x,y,z de P1 a XMM1
"MOVUPS %1, %%XMM2\n\t" //Sube x,y,z de P2 a XMM2
"MOVUPS %2, %%XMM3\n\t" //Sube x,y,z de P3 a XMM3
"MOVUPS %%XMM2, %%XMM4\n\t" //Copia P2 en XMM4
"SUBPS %%XMM1, %%XMM4\n\t" //Resta P2-P1 (XMM4)
"MOVUPS %%XMM3, %%XMM5\n\t" //Copia P3 en XMM5
"SUBPS %%XMM1, %%XMM5\n\t" //Resta P3-P1 (XMM5)
//--- Producto vectorial de P2P1 y P3P1 (determinante 3x3)
"PSHUFD $0xC9, %%XMM4, %%XMM6\n\t" //P2P1 en XMM6 (Y,Z,X)
"PSHUFD $0xC9, %%XMM5, %%XMM7\n\t" //P3P1 en XMM7 (Y,Z,X)
"MULPS %%XMM4, %%XMM7\n\t" //P2P1 * P3P1(Y,Z,X)
"MULPS %%XMM5, %%XMM6\n\t" //P3P1 * P2P1(Y,Z,X)
"SUBPS %%XMM7, %%XMM6\n\t" //Resta XMM6 y XMM7
//--- Módulo del vector resultado
"MOVUPS %%XMM6, %%XMM7\n\t" //Resultado a XMM7
"MULPS %%XMM6, %%XMM7\n\t" //Cuadrado de componentes XMM7
"HADDPS %%XMM7, %%XMM7\n\t" //Suma parcial de coordenadas
//... (A+B,C+D,A+B,C+D)
"HADDPS %%XMM7, %%XMM7\n\t" //Fin la suma
//... ((A+B)+(C+D),...)
"SQRTPS %%XMM7, %%XMM7\n\t" //Raíz cuadrada del resultado
"MOVUPS %3, %%XMM6\n\t" //Sube 0.5 a XMM6
"MULPS %%XMM6, %%XMM7\n\t" //Divide el resultado entre 2
"MOVUPS %%XMM7, %4" //Pone el resultado en R {R(0)}
:"=m"(P1), "=m"(P2), "=m"(P3), "=m"(DIV), "=m"(R)
:"m"(R)
);
return;
Como siempre lo mejor es hacer pruebas en las máquinas que vayan a ejecutar el código. Gran parte de la mejora en velocidad viene del hecho de reducir los accesos a memoria (todo lo que sea pasearse por el bus de la placa base es siempre como subirse a una tortuga…) y no salir de los registros de la CPU.
De esto hace ya algunos años, la extensión SSE4.2 incluye una instrucción para calcular el producto vectorial (nos ahorramos 5 instrucciones en el código dado). Si la integráis con la intersección plano/triángulo y la adaptáis a una CPU contemporánea de gama alta, vuestro raytracer podrá casi volar… 😉
Leer más