Interpolación (1)

A veces conocemos el valor de una función f(x) en unos puntos dados siendo pero no tenemos una expresión que nos permita calcular f(x) en otro punto cualquiera.

La misión es obtener una estimación de f(x) para un x cualquiera de modo que la curva que dibujemos sea "suave" entre los puntos y en su alrededor. La forma en que aproximemos las funciones tiene que ser general porque no sabemos en principio con que función nos estamos enfrentando. Las aproximaciones mas comunes o al menos las que aprendí en mis clases de calculo son (1)Polinomios (2)Cocientes de polinomios y (3)Funciones trigonométricas.

Hay funciones que tienen un "mal comportamiento" con la interpolación, para saber si el comportamiento ha sido bueno o malo es útil tener una estimación del error, hay que notar que estamos presumiendo que las funciones son suaves y continuas, lo que es perfecto para el movimiento real.

La idea básica es encontrar una función que encaje en los puntos dados y luego evaluarla para otro punto deseado, la mayoría de las aproximaciones empiezan con un punto y van haciendo correcciones al ir agregando otros nuevos puntos, esto nos da una complejidad de

La misión es obtener una estimación de f(x) para un x cualquiera de modo que la curva que dibujemos sea "suave" entre los puntos y en su alrededor. La forma en que aproximemos las funciones tiene que ser general porque no sabemos en principio con que función nos estamos enfrentando. Las aproximaciones mas comunes o al menos las que aprendí en mis clases de calculo son (1)Polinomios (2)Cocientes de polinomios y (3)Funciones trigonometricas.

Hay funciones que tienen un "mal comportamiento" con la interpolación, para saber si el comportamiento ha sido bueno o malo es útil tener una estimación del error, hay que notar que estamos presumiendo que las funciones son suaves y continuas, lo que es perfecto para el movimiento real.

La interpolación local con un numero finito de vecinos no es continua en las derivadas, si necesitamos que sea continua en las derivadas tenemos que usas la no localidad, por ejemplo con splines, es decir entre dos puntos dados definimos un polinomio que los interpola.

Pequeñas Mejoras

A raíz de las ultimas pruebas descubrí algunos fallos, el mas evidente era el que perdiese a los objetos que se movían con una velocidad "no entera", para solucionarlo se me ocurrió agregar un acumulador de velocidad que funciona de la siguiente forma:

velxAcum=velx+velxAcum-(int)(trunc(velx+velxAcum));
velyAcum=vely+velyAcum-(int)(trunc(vely+velyAcum));

De forma que guarda la porción del movimiento que por no ser entera no se ha podido aplicar a un solo paso, pero que al sumarse con otras partes no enteras dara finalmente con una solución mejor que si las ignoramos.

Este es el resultado:




Los saltos que se observan son debido a que el centroide se calcula desde las partículas que se dispersan aleatoriamente. No esta mal teniendo en cuenta que el objeto es invisible al filtro durante 80 frames y aun así la estimación cae sobre el objeto :)

Probando el filtro

Finalmente solucione el problema del circulo y la espiral, la idea es que la x no tiene porque avanzar un pixel en cada iteración, así que la solución que se me ocurrió es rellenar los espacios en que la y avanza mas de un pixel interpolando de alguna forma los pasos intermedios en que la x estaría estática.

Una vez conseguido eso pasamos a las pruebas:

Prueba 1: El rectángulo




Valoración: Satisfactoria, el filtro hace exactamente lo que tiene que hacer y no pierde el objeto aunque "desaparezca", cuando el objeto está rojo el filtro no lo reconoce, esto ayuda a que se pueda evaluar el funcionamiento ya que yo si lo puedo ver (es un super poder que tengo).

Prueba 2: La Recta




Valoración: Mal, muy mal, lo he puesto a posta para que lo haga mal, el problema es que la ecuación de esta recta es y = 50 + 1.2x, dado que el filtro esta trabajando sobre cantidades discretas pierde ese 0.2 a cada iteración, por eso se retrasa, esta en mi lista de TODO arreglar ese aspecto del filtro.

Prueba 3: El Circulo




Valoración: Regular, el filtro funciona como era de esperar, el problema es que esto no es suficiente para seguir al objeto, sigue la recta perpendicular al circulo, esto no pasaría o mejoraría si usase de alguna forma la derivada de la función... esto aun son suposiciones.

Prueba 4: La Espiral




Valoración: Similar al circulo, este es solo un caso ligeramente mas general al del circulo pero quedaba bonito y ademas me gustan las espirales que pasa. Se ve en el video que en algunos casos funciona mejor que en otros, esto depende de la zona de la circunferencia en que nos encontremos ya que en las zonas que estan proximas a los ejes del plano se recorren pixeles sucesivos en una misma linea, es decir solo hay movimiento en el eje x o en el eje y pero no en ambos.

Prueba 5: La Mosca




Valoración: Este es el prototipo de caso en que la realimentación no funciona, ya que el movimiento cambia rápidamente de dirección de forma quasi-aleatoria. No se me ocurre ninguna forma de hacer que esto funcione así que de momento tendrá que esperar.

Prueba 6: Función Sinusoide




Valoración: En este caso existe el mismo problema que sobre la recta, dado que estamos aproximando reales con enteros tenemos el problema de perdida de precisión, aun así no esta mal como primera aproximación.

                      uuuuuuu
                 uu$$$$$$$$$$$uu
              uu$$$$$$$$$$$$$$$$$uu
             u$$$$$$$$$$$$$$$$$$$$$u
            u$$$$$$$$$$$$$$$$$$$$$$$u
           u$$$$$$$$$$$$$$$$$$$$$$$$$u
           u$$$$$$$$$$$$$$$$$$$$$$$$$u
           u$$$$$$"   "$$$"   "$$$$$$u
           "$$$$"      u$u       $$$$"
            $$$u       u$u       u$$$
            $$$u      u$$$u      u$$$
             "$$$$uu$$$   $$$uu$$$$"
              "$$$$$$$"   "$$$$$$$"
                u$$$$$$$u$$$$$$$u
                 u$"$"$"$"$"$"$u
       uuu        $$u$ $ $ $ $u$$       uuu
      u$$$$        $$$$$u$u$u$$$       u$$$$
       $$$$$uu      "$$$$$$$$$"     uu$$$$$$
     u$$$$$$$$$$$uu    """""    uuuu$$$$$$$$$$
     $$$$"""$$$$$$$$$$uuu   uu$$$$$$$$$"""$$$"
      """      ""$$$$$$$$$$$uu ""$"""
                uuuu ""$$$$$$$$$$uuu
       u$$$uuu$$$$$$$$$uu ""$$$$$$$$$$$uuu$$$
       $$$$$$$$$$""""           ""$$$$$$$$$$$"
        "$$$$$"                      ""$$$$""
          $$$"                         $$$$"

Algo de Geometria I

Mi tarea de hoy ha sido programar pequeñas funciones en C++ que sean capaces de devolver coordenadas en pantalla que representen ciertas funciones.

Sen(x)

El problema con esta función es que esta definida entre -1 y 1, por lo que hace falta escalar la imagen y discretizarla.

y=a+(int)round(b*sin(gradosAradianes(i)));



Rectangulo:
Esta funcion es tan simple como concatenar 4 bucles que lo describan.


Circunferencia:




Al despejar la y obtenemos lo siguiente:



Esto debería generar un circulo pero como se puede comprobar no lo genera:




Esto se debe a que el avance de las x es constante cuando en realidad no tiene por que ser así, es decir que aunque describe un circulo correctamente hay zonas en que se mueve mas rápido de lo que cabria esperar debido al efecto de la discretización.

Algunas imágenes contienen oclusiones totales o cambios de color, esa característica se usara para probar el filtro.

Todas las gráficas fueron sacadas de la Wikipedia

Diseño Modular

Después de muchas vueltas y mas vueltas me he decidido por un diseño para el nuevo sintetizador de vídeo mejorado (mágico, plus, maestro del universo!), la idea es que el usuario tenga que escribir un pequeño fichero describiendo el tipo de video que desea y se genere automáticamente.

Un ejemplo del fichero de usuario sería:

NFrames=340
Width=400
Height=400
FPS=1
ObjectSize=40
Function=/home/zenko/workspace/senoidal/Debug/senoidal
args=100 200 0 340
Color=/home/zenko/workspace/colores/Debug/colores
Sections=200 255 255 255 50 0 0 0 90 255 255 255

Donde cada cosa es exactamente lo que parece según la descripción.

Function es la ruta a un ejecutable que recibe 4 argumentos a, b, inicio y fin, con esto ejecuta la función:

for(i=primero;i < ultimo;i++)
{
 y=a+b*sin(gradosAradianes(i)));
}

Color es la ruta de un ejecutable que devuelve ternas R G B y recibe n parámetros de tal forma que si dividimos los parámetros en grupos de 4 tenemos:

1) Numero de ternas
2) R
3) G
4) B

y nos devuelve una lista de N ternas.

Finalmente el script toma la salida de ambos programas y pega los ficheros en uno solo tal que guarde la sintaxis esperada por el generador.

#Numero de frames del vídeo
4
#Ancho del vídeo
400
#Alto del ideo
400
#FPS para el vídeo
1
#Tamaño del objeto
40
#x y R G B para cada frame
1 102 255 255 255
2 103 255 255 255
3 105 255 255 255

De esta forma cada vez que queramos generar un vídeo con una función extraña o exótica solo hará falta escribir la función y todo el resto del código sera reutilizable.

Sintetizador de video

Una vez conseguido el algoritmo básico de seguimiento se inicia la etapa de pruebas en donde necesitaré muchos vídeos de distintos tipos, es una buena practica tener un sintetizador que me cree vídeos a medida para comprobar el seguimiento.

En principio las funciones escogidas son:

a*Sen(x)+b
a*x+b
(x-h)^2 + (y-k)^2 = r^2
Espiral de Arquímedes

Esta será la primera batería de pruebas para comprobar el rendimiento del filtro.

Prueba usando 100*Sen(x)+200



El vídeo es bastaaaante lento, pero es parametrizable así que se puede poner a la velocidad deseada.

Falta escribir algunas de las funciones y modelizar las oclusiones, una vez conseguido esto necesito también alguna forma de tomar estadísticas del seguimiento del objeto.

Evitando la dispersion desde un centro

Un efecto extraño que teníamos en anteriores versiones es que la dispersión generaba un cuadrado casi perfecto, se debe a que las todas las partículas se dispersan desde un mismo centro, esto no es tan buena idea ya que se supone que estamos perdiendo precisión con respecto a cuando tenemos medida, y en ese momento la nube es mas dispersa.

Para solucionarlo he cambiado la forma en que se calcula la velocidad y los centroides, el centroide se sigue calculando a partir de las n partículas en vez de moverlo por la inercia, las partículas se dispersan con la función aleatoria estándar como de costumbre pero la velocidad no se actualiza de nuevo hasta obtener una nueva medida.

Esto que parece ser la idea mas lógica y obvia no se me ocurrió hasta tiempo después de estarlo pensando!

Seguimiento exitoso

Después de mucho trabajo arduo y todo eso que se suele decir por fin el filtro ha seguido un objeto correctamente mezclando las dos aproximaciones anteriores, si hay medida usa su esquema estándar de filtro de partículas y cuando no hay medida usa el esquema secundario realimentado.

Durante el tiempo que no tiene medida se va aumentando la incertidumbre linealmente en función del tiempo que ha pasado desde la ultima medida, siendo la incertidumbre inicial la dispersión que se estaba usando en el pf.

Al aumentar la incertidumbre desde la estimación de la realimentación obtenemos un cuadrado en vez de un circulo como cabria esperar para un radio constante, esto se debe a la forma en que está programado.

Usando 150 partículas



El cuadrado que se genera se debe a dos cosas, la primera es que todas las partículas se dispersan desde el mismo centro, la segunda es que todo el área que cubren es equiprobable ya que no tenemos medida.

Usando 10 particulas


Una prueba curiosa es ver como funciona el filtro con tan solo 10 partículas, para mi sorpresa termino encontrando al objeto y siguiéndolo, las fluctuaciones del centroide (cuadrado rojo) mientras no tiene medida son debidas a que se esta calculando el centroide real de las partículas no el estimado debido a la inercia, Como las partículas se dispersan con una función aleatoria normal si usamos muchas partículas ambos centroides convergen, aunque para pocas partículas da un poco de mala pinta.

El problema del aleatorio no tan aleatorio

#include 
#include 

int main()
{
    int i,cont=0,aux=0;
    for(i=0;i<300;i++)
    {
        aux=rand()%7-3;
        if(aux>0)
            cont++;
        else if (aux<0)
            cont--;
    }
    printf("cont= %d\n",cont);
}

Al ejecutar este código en un bucle durante 300 iteraciones me encuentro con que cont vale -13, es decir que han habido 13 números negativos mas que positivos, si tuviese un aleatorio realmente aleatorio la cantidad de negativos y positivos debería ser igual, el problema es que al salir mas números negativos que positivos el centro de la distribución sufre un retroceso.

Si estamos siguiendo un objeto con partículas ese retroceso significa perder el objeto, por lo que ahora la opción es aumentar el grado de incertidumbre sin usar ninguna función aleatoria, ya que perdemos incertidumbre de una forma constante desde que perdemos el objeto.

El siguiente objetivo es crear una función Esparcir() o algo similar que distribuya uniformemente x partículas en una zona de espacio.

O incluso mejor mover el centroide y luego dispersar las partículas a partir de el sin volverlo a actualizar si no tenemos una medida!

La informacion es poder

La inercia solo funciona si nuestro sistema inercial tiene suficiente información, es decir el sistema tiene una "creencia" de la velocidad a la que se mueve la partícula que está siguiendo, esta creencia sera mas o menos informada y según sea mas o menos podrá simular el efecto de movimiento con mayor o menor precisión.

Tras 5 iteraciones de dispersión normal:



Tras 10 iteraciones de dispersión normal


Tras 15 iteraciones de dispersión normal se obtiene información suficiente



Lo importante de esto es notar que no podemos pedir una buena estimación si antes no hemos proporcionado cierta información para que siga una cierta inercia.

Actualización: He subido los vídeos que ilustran lo anterior para subsanar el error de ayer que subí tres vídeos iguales.

Ligeramente atascado con los modelos de dispersión

He estado trabajando con nuevos modelos de dispersión para las partículas, y viendo las simulaciones me di cuenta de algunos errores en el anterior enfoque, ni la dispersión lineal ni la inercial funcionan como esperaba, así que habrá que hacer algunos cambios.

Dinámica Inercial

Cada partícula guarda su propia inercia, esto es un "error" porque la inercia del movimiento se ve aproximada por la inercia de la estimación, la inercia de las partículas individuales no aproxima nada y se pierde el objeto.

Ademas al usar la inercia de las partículas individuales se tiene que la velocidad de dispersión es constante, aunque la velocidad del movimiento no lo sea.

Para solucionar esto hay que guardar un historial de estimaciones por grupos de dispersión.

Una vez conseguido que el historial nos devuelva la velocidad media de desplazamiento de las ultimas N iteraciones conseguimos un comportamiento bastante mejor:




El siguiente cambio es conseguir un método que represente la perdida de seguridad que se tiene cuando desaparece una partícula.

Dispersiones Multiples

Hoy la tarea es ver como funcionan diferentes tipos de dispersiones, la dispersión normal que es dado un radio máximo usamos una función aleatoria normal para esparcir las partículas en dicho radio, la segunda opción es seguir la recta de mínimos cuadrados que hemos creado en anteriores entregas y finalmente habrá partículas que únicamente tengan inercia, es decir que continúen con la velocidad que ya tenían.

El ultimo enfoque no me gusta tanto porque al cambiar el movimiento del objeto no hay forma de que la dispersión lo encuentre, sin embargo si el objeto desaparece nos daría una idea de por donde esta en dicho momento porque las cosas en el mundo real tienen cierta inercia.

He hecho una simulación, la partícula roja muestra la dispersión estándar, la partícula azul sigue la recta de mínimos cuadrados y la verde es inercial.



Cada movimiento funciona con 100 partículas independientes, aunque como veremos en un futuro este no es el mejor enfoque posible por lo que ya comente.

La recta de mínimos cuadrados

lo que queremos es aproximar una distribución de N puntos por una recta de forma que la recta este lo mas cerca posible de todos los puntos, lo primero es decir a que nos referimos por cerca, esto es que la distancia entre el todos los puntos y la recta sea mínima.




Tenemos que encontrar entonces los valores de a y b para que el valor siguiente sea minimo:



Es decir

En este caso lo que tenemos que encontrar es a y b, para eso supongamos que las derivadas parciales son 0.





Es decir que:



Y substituyendo obtenemos a y b:







Mucho mejor explicado aunque mucho mas en ingles

El Historial

La idea del filtro de partículas es (en nuestro caso) seguir un objeto a través de un vídeo, por lo que es útil tener una estimación puntual para cada iteración, si guardamos esa iteración podemos realizar muchos "bonitos" cálculos sobre la trayectoria del objeto.

Con esa idea y también la del mínimos esfuerzo me puse a buscar una lista enlazada para guardarlo y me encontré con la libreria estándar de c++ que tiene una interfaz genial para las listas enlazadas y nos evita tener que escribir una cada vez que la necesitamos.

El uso es de lo mas simple:

Declarar una lista:

list < CvPoint >  historial;

Insertar un objeto:

historial.push_back(punto);

Iterar en la lista:

list < CvPoint > ::const_iterator actual = historial.begin();
while(actual!=historial.end())
actual++;

Si después de cada Resampling guardamos el centroide del conjunto de partículas tendremos un historial para hacer todos los cálculos que queramos.

La presentacion

El jueves 13 (13/12/2007) se celebro una reunión para explicar la parte hecha del proyecto y salieron varias ideas que he estado implementando poco a poco:

1. El filtro de partículas debería dar una "estimación" de donde esta el objeto que se esta buscando tras cada iteración, he llamado a esto Centroide, es básicamente una partícula que nos dice la mejor estimación que tiene del objeto para cada frame.
2. Las partículas deberían incorporar una "velocidad" para poder guardar un historial de lo que han hecho en el ultimo movimiento.
3. Para las siguientes pruebas voy a usar unas 100 partículas para que haya menos posibilidades de que ninguna caiga dentro del objeto.
4. Ahora tengo un sintetizador de vídeo que usa opencv, cortesía de Juanjo Pantrigo.

Bueno basándome en esas ideas y aunque con no demasiado tiempo por la cantidad de trabajo que he tenido la ultima semana he preparado un filtro funcional, de momento es muy sencillo y la única función discriminante que usa es: todos los pixeles deben estar a 0 salvo los del objeto.

Como he aprendido a usar el vlc para grabar vídeos desde el escritorio ahí va una de las pruebas que he hecho al filtro:




Aunque no se vea muy bien un cuadro rojo sigue al blanco durante el vídeo de una forma bastante aproximada :P

Implementacion parcial

Esta semana he conseguido una implementación parcial de lo que seria el conjunto de partículas y lo he probado con una imagen fija, la idea esta conseguida pero faltan bastantes mejoras.

La operación mas básica es lanzar las partículas sobre posiciones aleatorias, existe el problema de que ninguna caiga sobre el objeto que estamos buscando, para lo que se me ocurren dos soluciones, iterar hasta que alguna caiga o lanzar un numero mas grande de partículas.


Una vez están todas las partículas iniciadas, debemos discriminar las que han caído sobre el objeto de estudio de las que no, para esto hace falta una función que nos diga que es objeto y que no, en el caso que he probado es muy fácil porque lo que no es objeto esta todo en negro, pero puede no ser tan evidente y habría que hacer por ejemplo un vector de características del objeto. Como la función puede variar y la idea es que el filtro sea reutilizable he definido la función así:

void MultFuncion(IplImage* img,float (*evaluacion)(IplImage* img, Particula p));

Se le pasa un puntero a la función que se encarga de reconocer el objeto y esta función es responsabilidad del usuario del filtro (al menos por ahora).

Lo anterior solo cambia el peso de las partículas, pero lo que necesitamos es que las partículas "mas aptas" tengan hijos de forma proporcional a su peso, esto es lo que hace Resampling, tal como explique en otro post anterior.


Nos falta aplicar la dinámica, por el momento esta definida como sigue:

void Dinamica(int desviacionMax,int ancho, int alto);

La idea es cambiarla para que sea capaz de simular diferentes distribuciones, y para poder hacer experimentos con computación basada en palabras.

No entiendes realmente algo a menos que seas capaz de explicarselo a tu abuela.

Albert Einstein

Fin de la transmisión

Lanzando particulas

La idea de hoy es la siguiente, refinar las partículas para poder lanzarlas de forma aleatoria sobre una imagen, hay algunas cosas que hay que tener en cuenta:

• Las partículas deberían tener el mismo tamaño
• No se debe dibujar ninguna partícula hasta que haya terminado el estudio de la imagen!
• Las partículas deben estar enteramente dentro de la imagen

Aquí un resultado improvisado:

Ahora necesito encontrar una función que sea capaz de dar una medida, dada una partícula que tanto se parece su contenido a la imagen que se está buscando, que genere una "firma", algo así como una cadena MD5 que al ser comparada con la del objeto original nos de un porcentaje de concordancia.

Sigo en la búsqueda...

Programando Particulas

Mi tarea últimamente ha sido programar las partículas, aunque no he tenido mucho tiempo esta semana por mis "múltiples ocupaciones", en una primera aproximación tengo que hacer tres clases, una para la partícula, otra para un conjunto de partículas que implemente las operaciones normales de conjuntos de partículas y otra de filtro de partículas:


El diagrama es algo chapucero, solo es para dar una idea de lo que estoy intentando hacer.

Aquí una simulación con varias partículas usando valores aleatorios:


Algo que tengo que pensar es como representar el conocimiento sobre la posición siguiente de la partícula conocida la actual.

El mundo de las particulas (4)

La ultima de las operaciones básicas sobre partículas se llama resampling (¿remuestreo?), esto viene de la idea de que una distribución de probabilidad p(x) se puede representar de muchas formas, pero algunas son mas eficientes que otras, esto tiene su definición y explicación matemática (que me salto así porque si).

La idea clave es que para que una representacion sea eficiente la mayoría de las partículas deberían tener pesos iguales, y el resampling es una forma de conseguirlo.

Supongamos que tenemos un conjunto de partículas S, después de aplicar el resampling obtendremos otro S' que sera como sigue:



S' sera otro conjunto de partículas con:

con probabilidad


Donde escoger es independiente cada vez.

Tambien se puede hacer resampling determinista, es decir que el número de partículas escogidas es directamente proporcional al peso de las antiguas sin ninguna aleatoriedad, eso presenta algún problema matemático con el que como no nos encontraremos haremos como que no importa.

Supongamos que tenemos un conjunto de partículas S como de costumbre y vamos a definir uno S' después de hacer un resampling determinista, antes tenemos que definir los pesos acumulados como sigue:



Entonces obtendremos un S' cuyas partículas serán:

donde es el j mas pequeño tal que


Mas gráficamente:


El resultado de la operación seria algo similar a esto:

El mundo de las particulas (3)

Aplicando la dinamica, esta es la segunda operación básica sobre una partícula, la idea es que tenemos alguna idea de donde va a estar la partícula en un instante t+1 si sabemos donde estaba en el instante t, esto lo expresamos con una distribución de probabilidades como de costumbre.

Tenemos nuestro conjunto de partículas y una función de densidad:



La operación nos devuelve otro conjunto de partículas:



La idea es que las partículas se dispersen siguiendo la función de densidad que hemos dado, esa función representa lo que sabemos del movimiento, si no sabemos nada sera la distribución uniforme aunque parece una buena idea usar una Gausiana, habrá que probar que funciona mejor.



En el vídeo se aplica una dinámica sobre las partículas para simular un fluido, impresiona bastante, es para otro uso pero la idea es la misma o al menos eso parece porque solo he podido ver el vídeo y no la implementación.