Curva de Koch con Numpy y Pylab

Esta vez venimos con código para generar la famosa Curva de Koch 
http://es.wikipedia.org/wiki/Copo_de_nieve_de_Koch , este es solo la curva..no el copo ...
pero el sistema es el mismo. 
Como verán el código es genérico porq tiene toda la intención de usarse 
con otras gramáticas de Lin. .. pero para eso falta un poco de pulido.
Los puntos clave pasan por la recursión, el cursor y los niveles :).
Espero que les sea útil.
#! /usr/bin/python2
# -*- coding: utf-8 -*-

import numpy as np
import pylab as pl

class Cursor:
    def __init__(self,x,y, angulo):
        self.x = x
        self.y = y
        self.ang = angulo

class SistemaL:
     def __init__(self,gramatica, inicial, distancia, nivel):
        self.gramatica = gramatica
        self.inicial = inicial
        self.distancia = distancia
        self.nivel = nivel

     def iniciar(self):
        self.cursor = Cursor(0,0,0)
        pl.plot([self.cursor.x, self.cursor.x+ self.distancia],
                [self.cursor.y, self.cursor.y])

     def generar(self, nivel, cursor, distancia):
        if nivel == 0:
           a = (cursor.x, cursor.y)
           b = (cursor.x + distancia*np.cos(cursor.ang),
                cursor.y + distancia*np.sin(cursor.ang))
           print a, "-->", b
           pl.plot([a[0],b[0]],[a[1],b[1]])
           cursor.x = b[0]
           cursor.y = b[1]
        else:
           print "nivel = ", nivel
           for i in self.gramatica[self.inicial]:
               if i == 'F':
                    self.generar(nivel-1, cursor, distancia/3)
               elif i == '+':
                    cursor.ang = cursor.ang - np.pi/3
               elif i == '-':
                    cursor.ang = cursor.ang + np.pi/3

if __name__ == "__main__":
     miSistema = SistemaL(gramatica = {'F':"F-F++F-F"},inicial = 'F',
                          distancia = 10000, nivel = 5)
     miSistema.iniciar()
     miSistema.generar(miSistema.nivel, miSistema.cursor,
                       miSistema.distancia)
     print "termine con todos los niveles"
     pl.show()

Yyyy unas imágenes de como sería nuestro resultado ...

El pase de diapositivas requiere JavaScript.

Utilizando python-perceptron para posicionar vertebras ( o lo que gustes …)

Usando la biblioteca http://code.google.com/p/python-perceptron/ con unos minimos cambios que se pueden encontrar en mi repositorio de gitorious .. quizas solo eran necesarios para mi puntualmente … y como esta la biblioteca les sirve perfectamente …

En esta ocasión levantamos los datos de un entrenamiento en particular  dado por un experto en el tema en una planilla de calculos …

los cargamos en nuestra red .. la entrenamos y la colocamos a clasificar …  y para darnos una mejor idea implementamos con matplotlib un bonito gráfico


#! /usr/bin/python2
# -*- coding: utf-8 -*-

import perceptron
import xlrd
import numpy as np
from pylab import *

#---##---## ---##---##---##---##---##---##---##---##---##---##---##---
def cargarPatronEntrenamiento(path):

      book = xlrd.open_workbook(path)
      sheet = book.sheet_by_index(0)

      patron = {}
      entrada = []

      for i in range(1,sheet.nrows):
             entrada = []
             for j in range(1,sheet.ncols):
                  entrada.append([sheet.cell_value(i,j)])
             patron[i-1]= (np.array(entrada))

return patron
#---##---## ---##---##---##---##---##---##---##---##---##---##---##---
def entrenar(network, patron):

"""Con la matriz de entrenamiento proporcionada entrenar para modificar la
matriz de pesos"""

      for j in range(50): #TODO ...
           for i in patron.iterkeys():
               network.learn(patron[i], i)

#---##---## ---##---##---##---##---##---##---##---##---##---##---##---
if __name__ == "__main__":

      per = perceptron.Perceptron(35,23)
      patron = cargarPatronEntrenamiento("../data/entrenamiento.xls")
      entrenar(per, patron)
      forPlot = []

      for i in patron.iterkeys():
              output = per.classify(patron[i])
              forPlot.append(float(output[1]))

      ylim(0.20,1.0)
      xticks( arange(0,23), ('C2', 'C3', 'C4', 'C5', 'C6', 'C7', 'C8',
                             'C9', 'C10', 'C11', 'C12', 'C13', 'D1', 
                             'D2', 'D3', 'D4', 'D5','D6', 'D7', 'D8', 
                             'D9', 'D10', 'D11') )
#Estos son valores fijos dados por el experto en tema

      plot(arange(0,23),[0.80,0.80,0.60,0.60,0.80,0.80,0.80,0.80,
                        0.60,0.60,0.6,0.60,0.80,0.80,0.80,0.80,0.80,
                        0.70,0.70,0.50,0.50,0.50,0.50,],'x-',
                        label='Confianza proporcionada por humano')

      plot(arange(0,23),forPlot,'o-',label='Estimados por la red' )

      y = []
      for i in range(23):
          y.append(0.40)

#este es el limite de confianza que suponemos ... (un poco bajo para que queden bien las estadisticas :P)

      plot(arange(0,23),y,'-',label=u'Límite de confianza')

      xlabel(u'Vértebras')
      ylabel('Porcentaje de confianza')
      title(u'Predicciones de posicionamiento de vértebras')
      legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)

      show()

Y la salida quedaría algo asi …

Log 1 – Cargando BN con pyAgrum

A continuación veamos un ejemplo de una BN con su respectivas tablas de probabilidad…

Teniendo una idea general de la red bayesiana que queremos realizar (usemos la imagen como guía), pasemos a la implementación…


from pyAgrum import BayesNet, LabelizedVar, LazyPropagation

#Creamos nuestra red bayesiana vacia
myBN = BayesNet("Noche Calurosa")

#Cargamos los nodos

n = myBN.add(LabelizedVar('n','noche calurosa?',2))
c = myBN.add(LabelizedVar('c','cerveza?',2))
h = myBN.add(LabelizedVar('h','helado?',2))
a = myBN.add(LabelizedVar('a','llamar amigos?',2))

#LabelizedVar(nombre,comentario,cantidad de variables)

#insertamos los arcos
arcos = [(n,c),(n,h),(c,a),(h,a)]
for a in arcos:
 myBN.insertArc(a)

#Ahora cargamos nuestras tablas de probabilidades.

myBN.cpt(n)[:] = [0.5,0.5]

myBN.cpt(c)[:] = [[0.5,0.5],
 [0.1,0.9]]

myBN.cpt(h)[:] = [[0.8,0.2],
 [0.3,0.7]]

#otra forma de cargar la tabla puede ser:

bn.cpt(a)[0,0,:] = [1,0]

bn.cpt(a)[0,1,:] = [0.05,0.95]

bn.cpt(a)[1,0,:] = [0.90,0.10]

bn.cpt(a)[1,1,:] = [0.010,0.99]

Este código respecta a la creación de nuestra red bayesiana, más adelante veremos como realizar la inferencia y trabajar con
la evidencia.

Log 0 – Manipulación Básica de Imagenes

La suma de valores a cada pixel puede ser utilizado para obtener,
Ajuste de Contraste, sumando un valor a cada pixel de una imagen incrementa su valor, por lo tanto aumenta su brillo.


from PIL import Image

myIm = Image.open("tucan.jpg")

myImResta = myIm.point(lambda x: x-100)
myImSuma = myIm.point(lambda x: x+100)

myImResta.save("contraste-100.png")
myImSuma.save("contraste+100.png")

En esta imagen vimos suma y resta de valores constantes.

Otra operación interesante es el ‘blending‘ .. que consiste en sumar dos imagenes (ponderando una mas que la otra con un alpha constante).

from PIL import Image

fox = Image.open("fox.jpg")
buddy = Image.open("fox2.jpg")

# para realizar blending deben tener el mismo tamano
buddy.resize(fox.size)

out = Image.blend(fox, tucan, 0.40)
# out = image1 * (1.0 - alpha) + image2 * alpha

out.save("myFoxyBlending.jpg")

En cuanto a operaciones lógicas una de las mas simples, como lo es la Negación , podemos inferir facilmente que esta operación invierte la representación de la imagen.
En imagenes binarias, el negro se convierte en blanco y viceversa.
En imagenes de color o escala de grises, el proceso consiste en aplicar a cada pixel $(MAX – p(x,y))$, donde $MAX$ es el valor máximo posible.

from PIL import Image, ImageChops

myImC = Image.open("tucan.jpg") #imagen color

myGreyTucan = myImC.convert("L")
myBWTucan = myImC.convert("1")

out = ImageChops.invert(myImC) #Invertido Color
out2 = ImageChops.invert(myGreyTucan) #Invertido Color
out3 = ImageChops.invert(myBWTucan) #Invertido Blanco y Negro

También podemos jugar con los espacios de color …

La representación de los colores en una imagen digital se logra usando una combinación de uno o más canales.

La representación que utilizamos para almacenar los colores, especificando numero y ‘naturaleza’ del canal es conocido como espacio de color.

 from PIL import Image
 from pylab import subplot,bar,savefig

 myIm = Image.open("tucan.jpg")
 myHis = myIm.histogram()

 subplot(333)
 bar(range(256), myHis[0:256], ec = 'r')
 subplot(336)
 bar(range(256), myHis[256:512], ec = 'g')
 subplot(339)
 bar(range(256), myHis[512:768], ec = 'b')

 savefig("myHIsto.png")
 

RGB es un espacio tridimensional, por lo que podemos imaginarlo como un cubo, con ejes R,G y B, todos los ejes poseen el mismo rango (escalado de 0-255 para 1 byte por canal, asignando 24 bit para la representación de la imagen).

El color negro lo encontramos en el origen del cubo $RGB(0,0,0)$ , el cual corresponde a la ausencia de color y el color blanco en la esquina opuesta $RGB(255,255,255)$ indicando la máxima cantidad de los tres colores.
El espacio de color RGB esta basado sobre la porción del espectro visible de los humanos.

Otro espacio con el que podemos experimentar es HSV,

El espacio de color ‘apreciable’ es una forma alternativa para representar imagenes en su verdadero color de una forma más natural para el humano que el espacio RGB.
En este espacio nos concentramos en el :

  • Matiz, es la longitud de onda dominante del color, (rojo, azul o verde).
  • Saturación, muestra el nivel de “pureza” del color, la cantidad de “luz blanca” mezclada con el color.
  • Valor, es el brillo del color , luminancia.

La representación HSV de una Imagen , al igual que RGB es un arreglo tridimensional, todo pixel de la imagen $f(x,y)$ posee una tupla del tipo $(h,s,v)$.
Podemos ver en el siguiente ejemplo la división de estos tres canales.

import cv
myIm = cv.LoadImage("tucan.jpg")
myHSV = cv.CreateImage(cv.GetSize(myIm),8,3)
cv.CvtColor(myIm, myHSV, cv.CV_RGB2HSV)

cv.SaveImage("myHSV.png", myHSV)

#Con esto tenemos la trasnformacion de RGB a HSV
#Ahora dividimos los canales

h_myHSV = cv.CreateImage(cv.GetSize(myIm), 8, 1)
s_myHSV = cv.CreateImage(cv.GetSize(myIm), 8, 1)
v_myHSV = cv.CreateImage(cv.GetSize(myIm), 8, 1)

cv.Split(myHSV, h_myHSV, s_myHSV, v_myHSV, None)

cv.SaveImage("h.png",h_myHSV)
cv.SaveImage("s.png",h_myHSV)
cv.SaveImage("v.png",h_myHSV)

 y como en el caso anterior podemos realizar un merge y ver que se vuelve a la imagen inicial.
cv.Merge(h_myHSV, s_myHSV, v_myHSV,None, myIm)
cv.SaveImage("myMergeHSV.png",myIm)

Estos y los siguientes logs forman parte de ejemplos básicos  (de un trabajo) para que no queden en /dev/null 🙂