Octave plot custom grid

Necesitaba graficar algo, pero quería cambiar el espaciado de la rejilla a 0.25

P = dlmread("serial500.txt"," ");
plot(P(:,1),P(:,2),"*")
set(gca(),"xtick",[0:0.25:1])
set(gca(),"ytick",[0:0.25:1])
grid
xlabel("x")
ylabel("y")
title("500 pares de nums")

Pastilla para hacerse invisible

Hoy al fin probé la famosa pastilla para hacerse invisible, es impresionante la rapidez con que hace efecto. Os dejo un video, haber si alguien puede explicarme cómo es que funciona esta revolucionaria mezcla de ingredientes exóticos.

Implementación de k-means en OpenCV

k-means o k-medias es un algoritmo de agrupamiento (clustering), lo que significa que es usado para dividir un conjunto de datos en grupos de tal forma que los datos pertenecientes a un grupo tengan propiedades en común. También se dice que este algoritmo es de aprendizaje no supervisado, pues es capaz de clasificar datos sin tener un supervisor que durante la etapa de entrenamiento le indique la clase a la que pertenece cada dato.
Se puede usar este algoritmo para segmentar una imagen por color, clasificando los pixeles que tienen colores parecidos en un mismo grupo.

k-means funciona de la siguiente manera:

• Elige aleatoriamente k centros (uno por cada grupo)
• Se etiqueta cada dato como perteneciente al grupo cuyo centro es más cercano
• Se calcula el centro promedio de cada grupo y se repite el paso anterior.
• El algoritmo se detiene cuando el promedio de cada centro de grupo no cambia o se alcanza el número máximo de iteraciones.

El centro de cada grupo y las distancias quedan definidos por el problema.
En el caso de las imagenes RGB se define el centro como una tripleta que contiene un valor para cada color:rojo, verde y azul. La distancia que se usará en este caso es la euclídea.
Sólo hay que imaginar los colores como una coordenada en el espacio 3D y la distancia entre dos puntos como el largo de la recta que los une.

El resultado del algoritmo dependerá de los centros iniciales, si estos son escogidos aleatoriamente probablemente se obtendrán resultados diferentes en cada ejecución.

OpenCV trae implementado este algoritmo, sin embargo yo programé mi propia implementación porque el algoritmo es sencillo y me pareció interesante hacerlo.
Un poco de código:

  //calculo de centroides y matriz de clasificacion
for(it=0; it < iterations;it++) {
    for(i=0;i<size;i++) {
        red[i] = 0;
        green[i] = 0;
        blue[i] = 0;
        count[i] = 0;
    }
    for(j=0; j < image->height; j++) {
        for(i=0; i < image->width; i++) {
            tcolor = getColor(image,i,j);            
            dist = distance(centroid[0], tcolor );
            label = 0;
            for(k=1; k<size; k++) {
                temp = distance(centroid[k], tcolor );
                if( temp < dist) {
                    dist = temp;
                    label = k;
                }
            }
            ++count[label];
            red[label]   += tcolor.r;
            green[label] += tcolor.g;
            blue[label]  += tcolor.b;
            if(it == (iterations-1)) {
                labeled[i][j] = (uchar)label;
            }
        }
    }
    //actualizar centroides
    for(i=0;i<size;i++) {
        if(count[i]==0) continue; //no actualizar centroide si no se hay colores que promediar
        centroid[i].r = (uchar)(red[i]/count[i]);
        centroid[i].g = (uchar)(green[i]/count[i]);
        centroid[i].b = (uchar)(blue[i]/count[i]);    
    }
                    
}
//dejar con color de centroide mas cercano a cada pixel de la imagen a mostrar
for(i=0; i < image->width; i++) {
    for(j=0; j < image->height; j++) {
        image2->imageData[j*image->widthStep + i*image->nChannels + 2] = centroid[labeled[i][j]].r;
        image2->imageData[j*image->widthStep + i*image->nChannels + 1] = centroid[labeled[i][j]].g;
        image2->imageData[j*image->widthStep + i*image->nChannels ] = centroid[labeled[i][j]].b;
    }
}

centroid[] : contiene un color por cada clase
size: cantidad de clases
labeled[][]: almacena los resultados de la clasificación, un número de 0 a (size-1) que relaciona cada pixel con su clase.
red[],green[],blue[],cout[]: variables usadas para calcular el color promedio

En el código anterior, dados los centroides de cada clase se procede a identificar cada pixel con la clase cuyo color se encuentre más cercano a su color.
Posteriormente se actualizan los centroides de cada clase, para ello se promedian los colores de los pixeles que fueron etiquetados como pertenecientes a la clase.
El proceso se repite hasta completar la cantidad de iteraciones especificada en la variable iterations.
Código fuente completo: kmeans-2010.tar.gz

Resultados al aplicar a la imagen superior izquierda el algoritmo con 5 clases y 90 iteraciones:

Filtro de Sobel en OpenCV

Una forma más elaborada de detectar bordes es el filtro de Sobel, dado que usa una mayor cantidad de píxeles de la vecindad para calcular el borde.

#include <opencv/highgui.h>
#include <opencv/cv.h>

void sobel(IplImage *img) {
  int i,j,indice;
  uchar *pixels = (uchar*)(img->imageData);
  uchar *new_pixels = (uchar*)malloc(img->widthStep*img->height*img->nChannels*sizeof(uchar));

    uchar sup;
    uchar sup_izq;
    uchar sup_der;
    uchar inf;
    uchar inf_izq;
    uchar inf_der;
    uchar izq;
    uchar der;
    int A;
    int B;
    double SUM;
  //el origen esta en esquina superior izquierda del pantalla
 for( i=1 ; i < (img->width-1) ; i++) {
   for( j=1 ; j < (img->height-1) ; j++) {
         indice = j*img->widthStep + i*img->nChannels;
         sup = pixels[indice - img->widthStep]; //pixel superior al actual
            inf = pixels[indice + img->widthStep]; //pixel inferior al actual
            izq = pixels[indice - img->nChannels]; //pixel izq al actual
            der = pixels[indice + img->nChannels]; //pixel der al actual
            sup_izq = pixels[indice -  img->nChannels - img->widthStep]; 
            sup_der = pixels[indice +  img->nChannels - img->widthStep];
            inf_izq = pixels[indice -  img->nChannels + img->widthStep]; 
            inf_der = pixels[indice +  img->nChannels + img->widthStep];
            A = sup_izq + 2*sup + sup_der - inf_izq - 2*inf - inf_der;
            B = sup_izq + 2*izq + inf_izq - sup_der - 2*der - inf_der;
            SUM = A*A + B*B;            
            new_pixels[indice] = 255-(uchar)sqrt(SUM);

   }
 }
 memcpy(pixels, new_pixels, img->widthStep*img->height*img->nChannels*sizeof(uchar));
 free(new_pixels);
  
}

int main() {
 CvCapture *cap = cvCreateCameraCapture(0);
 IplImage  *img;
 IplImage  *imgGray;
 int notCreated=1;
 while(1) {
  img = cvQueryFrame(cap);
  cvNamedWindow("hahaha",CV_WINDOW_AUTOSIZE);
  if(notCreated) {
   imgGray    = cvCreateImage( cvSize(img->width,img->height), IPL_DEPTH_8U, 1 );
   notCreated = 0;
  }
  cvCvtColor(img, imgGray, CV_RGB2GRAY);
  sobel(imgGray);  
  cvShowImage("hahaha",imgGray);
  char key = cvWaitKey(40);
  if(key == 27)
   break;
 }
    cvReleaseImage(&img);
    cvReleaseImage(&imgGray);
 cvReleaseCapture(&cap);
 cvDestroyWindow("hahaha");
    return 0;
}

Efecto fotocopia en tiempo real

He estado jugando un poco con la librería de visión artificial de Intel, OpenCV. Programé el mismo efecto fotocopia que postie tiempo atrás, pero esta vez aplicandolo al video capturado por la webcam en tiempo real.
Esto no es difícil de lograr, sin embargo es interesante ver los resultados.

Aca el código:

#include  <opencv/highgui.h>
#include  <opencv/cv.h>

void fotocopia(IplImage *img, int difMax) {
    int i,j,difh,difw,indice;
    for( i=0 ; i < (img->width-1) ; i++) {
        for( j=0 ; j < (img->height-1) ; j++) {
            indice = j*img->widthStep + i*img->nChannels;
           /*diferencia brillo en pixeles de una misma columna
            widthStep es el tamaño de una fila en bytes, 
             ya que la imagen esta en un arreglo unidimensional.
              Al sumar widthStep estamos pasando a la siguiente fila
            */
            difh = ( (uchar*)(img->imageData) )[indice] - ( (uchar*)(img->imageData) )[indice+img->widthStep];
            if( difh < 0) {
                difh = difh * -1; 
            } 
        
            //diferencia brillo en pixeles de una misma fila
            
            difw = ( (uchar*)(img->imageData) )[indice] - ( (uchar*)(img->imageData) )[indice+1];
            if( difw < 0) {
                difw = difw * -1;
            }            
     
            //si hay mucha diferencia, es un contorno, hay que pintar el pixel negro
            if ( (difh + difw) > difMax) {        
                ( (uchar*)(img->imageData) )[indice] = 0 ;
            } else {
                ( (uchar*)(img->imageData) )[indice] = 255 ;
            }
        }
    }        
}

int main() {
    CvCapture *cap = cvCreateCameraCapture(0);
    IplImage  *img;
    IplImage  *imgGray;
    cvNamedWindow("hahaha",CV_WINDOW_AUTOSIZE);
    img = cvQueryFrame(cap);
    imgGray    = cvCreateImage( cvSize(img->width,img->height), IPL_DEPTH_8U, 1 );
    while(1) {
        img = cvQueryFrame(cap);
        cvCvtColor(img, imgGray, CV_RGB2GRAY);
        fotocopia(imgGray,5);        
        cvShowImage("hahaha",imgGray);
        char key = cvWaitKey(40);
        if(key == 27)
            break;
    }
    cvReleaseCapture(&cap);
    cvReleaseImage(&img);
    cvReleaseImage(&imgGray);
    cvDestroyWindow("hahaha");
}

Para compilar:

gcc `pkg-config --cflags opencv --libs opencv` fotocopia.c -o fotocopia

Aca un video de cómo se queda el efecto:

Editando imagenes con Octave

Me encanta el filtro fotocopia que trae el Gimp para dejar una foto como un dibujo, me dí cuenta de que es muy intuitivo hacer este efecto en Octave. Para ello sólo hace falta usar la función gradiente que toma una matriz y calcula la diferencia que hay entre dos elementos contiguos en cada fila y en cada columna.
Por ejemplo:

M =
1 2
3 4

[diffX, diffY] = gradient(M)

diffX =
1 1
1 1

diffY =
2 2
2 2

Se aprecia que diffX es la diferencia entre los elementos de las filas y diffY es la diferencia entre los elementos de las columnas.

Aplicando este concepto a las imagenes en escala de grises :

function fotocopia(imagen, salida)
im            = imread(imagen);
[diffX,diffY] = gradient(im);
fotocopia     = (diffX.^2 + diffY.^2).^(1/2);
fotocopia     = fotocopia * 3;
fotocopia     = 255 - fotocopia;
imwrite(fotocopia,salida);
endfunction

El operador . es para que se haga la operación a cada elemento de la matriz (elevar al cuadrado elemento por elemento)
La multiplicación por 3 es para realzar los bordes
La resta es para que el fondo quede blanco y el contorno oscuro (255 = blanco, 0 = negro).

fotocopia("plum.jpg","fotocopia.jpg");

Recuerdo

El año 2005 programé un juego de naves en C con la librería ncurses, usando caracteres ASCII.
Aprendí muchas cosas y sentí gran satisfacción al ver mi jueguito terminado. No creo que el código sea un ejemplo, pero fue para mí una pequeña victoria.
Os dejo el link para descargar el código y el ejecutable.
PD: Las balas no se veen porque estan mal programadas, su velocidad depende del procesador.

Puntero void*

Un puntero void* es un puntero genérico que se puede usar para hacer funciones que pueden recibir y/o devolver distintos tipos de datos según sea el caso. O sea que usando este tipo de puntero puedes hacer por ejemplo una función que trabaje con enteros si se le pasan punteros a enteros o con una estructura si se le pasan punteros a la estructura.
En el siguiente ejemplo se hace un casting a otro tipo de puntero para que así el compilador sepa cuantos bytes debe sumar a la dirección base, en este caso el casting es a char* por lo tanto el compilador sumará tamElem*k*sizeof(char) bytes a la dirección base.
Esta función recibe de argumento la función de comparación y busca en un arreglo el argumento llamado valor, en caso de hallarlo retorna el indice del elemento en el arreglo, si no encuentra al elemento en el arreglo retorna -1.

int busquedaLineal(void *valor, void *dirBase, int tamElem, int cantElem, int (*comparacion) (void*,void*)) {
    int k;
    for(k=0; k <  cantElem; k++) {
        void *elemActual = (char*)dirBase + tamElem*k;
        if( comparacion(valor,elemActual) == 0 ) {
            return k;
        }
    }
    return -1;
}

Editando sonido con Octave

Octave es un software muy parecido a Matlab, pero libre. Trae consigo un sin fin de funciones matemáticas, sobre todo para trabajar con matrices.
Una de las muchas cosas interesantes que se pueden hacer con este software es cargar un archivo de audio en una matriz y editarlo.

Un archivo de audio tiene tres propiedades :

Muestras por segundo (Hz)
Número de bits por muestra
Cantidad de canales.

Usando la función wavread se puede leer un archivo WAV. Un sonido monofónico se guardará en una matriz de una columna, un sonido stereo una matriz de dos columnas, etc.

Sonido original:



Efecto de eco :

#archivo entrada, archivo salida, sample rate(ej. 44100,22050,etc)
function ecoAudio(nombre_wav, salida_wav, samplerate_hz)
 wav_matrix = wavread(nombre_wav); #leer wav
 desp = int32( size(wav_matrix)(1)/20 );
 m2         = shift(wav_matrix,desp);
 wav_matrix = wav_matrix + m2;
 wavwrite(wav_matrix,samplerate_hz,salida_wav);
endfunction

Aplicando efecto de eco:




Cambio de rapidez :

#archivo entrada, archivo salida,sample rate,rapidez
function rapidezAudio(nombre_wav, salida_wav, samplerate_hz, rapidez)
 wav_matrix = wavread(nombre_wav);
 wavwrite(wav_matrix, samplerate_hz*rapidez, salida_wav);
endfunction

Reproduciendo a la mitad de su velocidad:



Reproduciendo al doble de su velocidad:


Reproducir al revés :

#archivo entrada, archivo salida, sample rate
function revertirAudio(nombre_wav, salida_wav,samplerate_hz)
 wav_matrix = wavread(nombre_wav);
 wav_matrix = flipud(wav_matrix); #revertir filas
 wavwrite(wav_matrix,samplerate_hz,salida_wav);
endfunction

Reproduciendo al revés:

Usando un entero como un arreglo de caracteres

Si tenemos una variable de tipo int por lo general ocupará 4 bytes de memoria, los que podremos usar literalmente a nuestro antojo. Por ejemplo podemos usar este espacio como un arreglo de 4 elementos de 1 byte cada uno.

Ejemplo uso de un entero como un arreglo de caracteres:

#include <stdio.h> 

int main() {

    int num;
    num=0;
    ( (char *)&num )[0] = 7;
    ( (char *)&num )[1] = 3;
    printf("%i %i\n",( (char*)&num )[0],( (char*)&num )[1]);
    *( (char *)&num ) = 2;
    *( (char *)&num + 1) = 4;
    printf("%i %i\n",( (char*)&num )[0],( (char*)&num )[1]);
    printf("%d\n",num);
}

Lo anterior imprime :

7 3
2 4
1026

La variable num queda como 1026 porque:
| 0000 0000 |0000 0000 | 0000 0100 | 0000 0010 |
| 0 | 0 | 4 | 2 |

La idea de este ejemplo es comprender mejor cómo funciona el lenguaje porque en la práctica usar un entero como un arreglo de char no hace que el código sea más simple.

Programación y ajedrez

El verano del 2008 comencé a programar un motor de ajedrez en C, que probablemente será un chimpacé al lado de cualquier otro motor. Aún no lo termino, pues le dedico tiempo solo los veranos y unos escazos días.
En los primeros programas de ajedrez el tablero se representaba usando un arreglo, luego se comenzó a usar una técnica que revolucionó el mundo del ajedrez, los bitboards.
Un bitboard es una representación del tablero usando un entero de 8 bytes, dado que cada byte tiene 8 bits se pueden representar los 64 casilleros del tablero con un bitboard.
La ventaja de esta representación es que los cálculos son mucho menos costosos en memoria y CPU porque un bitboard pesa mucho menos que un arreglo de 64 carácteres y se usan operaciones a nivel de bit como el AND,OR,XOR y shift para muchos de los cálculos.

Verificar si un casillero del tablero se encuentra ocupado por una pieza usando bitboards:

typedef uint64_t bitboard;
...
/*  square number 0-63 (a1-h8) */ 
int isOccuped(bitboard board, int square) {
    
    if( board & (1<< square) ) {
        return 1;    
    }
    return 0;
}

Para revisar si una casilla del tablero esta ocupada usando un arreglo tendriamos que recorrer el arreglo e ir revisando casilla por casilla, lo cual resulta más costoso y menos elegante.

Aritmética de punteros

Estaba indeciso acerca de que escribir en mi descuidado blog y decidí que a partir de hoy comenzaré a postear algunos fragmentos de código.
Sé que es muy fácil encontrar buenos ejemplos acerca del tema de hoy, pero más de alguno lo habrá olvidado y lo recordará por casualidad al llegar aquí.
C permite sumar o restar números enteros a un puntero, desplazando la dirección a la que este apunta la cantidad de bytes que ocupa el tipo de dato del puntero multiplicada por el número entero que se le ha sumado o restado.

#include <stdio.h> 
#include <stdlib.h> 

int main() {

    int* intArray;
 /*  (int *) indica al compilador que en la memoria reservada 
     se almacenaran enteros
     malloc (memory allocate) reserva (sizeof(int)*3) bytes de memoria 
     sizeof(int) retorna la cantidad de bytes que ocupa un int 
     (en mi caso 4 bytes)
     inicialmente intArray apunta al primer byte de memoria reservada
 */
    intArray  = (int *) malloc(sizeof(int)*3); 
   *intArray = 1;  //equivalente a intArray[0] = 1
//desplazar direccion apuntada por intArray sizeof(int) bytes
    intArray  = intArray + 1; 
   *intArray = 2;  
    intArray  = intArray +1; 
   *intArray = 3; 
    printf("%i %i %i\n",*(intArray - 2),*(intArray - 1),*intArray);
 /* si llamamos inmediatamente a free(intArray) obtendremos un error 
    en tiempo de ejecucion pues intArray ya no apunta al primer byte 
    de memoria reservada*/

 //intArray debe apuntar al 1er byte de la memoria reservada 
    intArray = intArray - 2; 
    free(intArray);

}

Muestra en pantalla:
1 2 3