Leer imagenes en formato .img en Matlab

Es muy común que datos científicos como imágenes medicas o astronómicas se encuentren con la extensión .img. Al menos hasta la versión de 2016, Matlab no tiene una función para extraerlas pero en esta entrada veremos como hacerlo leyendo directamente los datos binarios con fopen(). Para este ejemplo usaremos los datos topográficos del altímetro láser de la sonda espacial Mars global Surveyor que pueden ser descargados aquí. Es muy importante saber de antemano lo siguiente::

• Dimensiones de la imagen 
• Formato y precisión de los datos

Esta información es necesaria para poder redimensionar y obtener correctamente el array de datos leídos con la función fread(). En este caso usaremos el mapa topográfico de la superficie de Marte a una resolución de 4 pixeles por grado (archivo megt90n000cb.img). Estos datos tienen un archivo de documentación asociado llamado megt90n000cb.lbl donde podemos leer que la imagen tiene dimensiones de 1440 por 720 px con datos en formato big-endian int16. Sabiendo esto escribimos el siguiente programa:

%Abrir archivo .img
fid = fopen('megt90n000cb.img','rb');
dim = [1440 720];
img = fread(fid,dim,'int16','ieee-be');
img = img';

%Gráfica 2D
figure(1)
imagesc(img)
xticklabels = 0:30:360;
xticks = linspace(1, size(img, 2), numel(xticklabels));
set(gca, 'XTick', xticks, 'XTickLabel', xticklabels)
yticklabels = -90:30:90;
yticks = linspace(1, size(img, 1), numel(yticklabels));
set(gca, 'YTick', yticks, 'YTickLabel', flipud(yticklabels(:)))
colormap(hot), grid('on'), title('Mars Orbiter Laser Altimeter Data')
xlabel('Longitud(E)'),ylabel('Latitud')

%Gráfica 3D (Valle Marineris)
figure(2)
marineris = img(325:478,1079:1322);
surf(marineris,'EdgeColor','none')
xlim([0 200]), ylim([0 150])
title('Valles Marineris')

Plots:

Máquina de Estados Finitos en Python

Una máquina de estado o autómata finito (Finite State Machine en inglés) es un modelo computacional bastante útil en robótica y automatización. Para este ejemplo haremos una implementación por software de una máquina de estados en Python que puede ser fácilmente adaptada para algún proyecto en una Raspberry Pi. Esta implementación es bastante sencilla, para FSM's más robustas basadas en clases recomiendo checar aquí y aquí. El diagrama de estados para FSM de esté ejemplo será el siguiente:

Este tipo particular de FSM es llamada Máquina de Moore ya que la transición de estado sólo depende del estado actual y del valor de la entrada. Esta implementación en Python aprovechará una característica peculiar de los diccionarios que permite asociar una función con una llave [ref]. De esta manera los estados quedan definidos como funciones dónde podemos incluir todas las acciones que nos interese realizar durante dicho estado (o durante las transiciones). El programa completo:

# -*- coding: utf-8 -*-
"""
Created on Thu Feb  4 23:36:28 2016

@author: Rodolfo
"""
from time import sleep
from random import randint

#Variable global
estado = 'i'

#Estados
def EDOi(entrada):
    global estado
    print('Estado Inicial')
    #Transiciónes
    sleep(2)
    if entrada == 0:
        estado = 'i'
    if entrada == 1:
        estado = 0
        print(u'Transición hacia 0...')

def EDO0(entrada):
    global estado
    print('Estado 0')
    #Transiciones
    sleep(2)
    if entrada == 0:
        estado = 1
        print(u'Transición hacia 1...')
    if entrada == 1:
        estado = 2
        print(u'Transición hacia 2...')
        
def EDO1(entrada):
    global estado
    print('Estado 1')
    #Transiciones
    sleep(2)
    if entrada == 0:
        estado = 0
        print(u'Transición hacia 0...')
    if entrada == 1:
        estado = 2
        print(u'Transición hacia 2...')

def EDO2(entrada):
    global estado
    print('Estado 2')
    #Transiciones
    sleep(2)
    if entrada == 0:
        estado = 1
        print(u'Transición hacia 1...')
    if entrada == 1:
        estado = 0
        print(u'Transición hacia 0...')

#Finite State Machine (FSM)   
def FSM(entrada):
    global estado
    switch = {
       'i':EDOi,
        0 :EDO0,
        1 :EDO1,
        2 :EDO2,
    }
    func = switch.get(estado, lambda: None)
    return func(entrada)

#Programa Principal
while True:    
 FSM(randint(0,1))
 sleep(2)

Salida del programa en la consola en Raspberry:

En caso de que te interese una implementación en lenguaje C, aquí tengo otra entrada dónde explico como hacerlo.

Streaming de datos por TCP en Raspberry Pi + LabVIEW

 Si se desea hacer un streaming de datos desde un sistema de instrumentación remota como una estación meteorológica, sismológica o un telescopio robótico, puede resultar muy útil y económico tener un servidor TCP en un Raspberry Pi y visualizar los datos en tiempo real con LabVIEW desde cualquier lugar. Esto puede hacerse muy fácilmente utilizando el modulo socket para Python.

Un socket es una abstracción de programación para la representación de conexiones y permiten realizar una comunicación bidireccional. En Python se crea de la siguiente manera:

import socket
s = socket.socket()    #Instaciamiento para el objeto 's'
host = ''                      #Host local por defecto
port = 8006                #Puerto
s.bind((host,port))     #Asignación de dirección y puerto a la instancia
s.listen(5)                  #Número máximo de conexiones entrantes 

Para este ejemplo simularemos 2 señales que usaremos para enviar por TCP. En una aplicación real estas señales serían adquiridas mediate ADC's u otros dispositivos conectados a los pines GPIO de la tarjeta:

t = np.linspace(0,2*np.pi,100)
y0 = np.float16(np.sin(t))
y1 = np.float16(0.5*np.sin(t)*np.sin(2*t))

Se puede observar que estoy usando la función float16() para reducir la precisión de los datos para ahorrarle un poco de trabajo al cliente en LabVIEW pero pueden omitirla si requieren más dígitos significativos. Ahora, para parte principal del programa tenemos:

print u"Esperando conexión"
s, addr = s.accept()
print u"Conexión desde: ", str(addr)
i = 0
try:
 while True:
     s.send('CH0'+str(y0[i])+'\n')
     s.send('CH1'+str(y1[i])+'\n')
     i += 1
     if i>=99:
         i = 0
     time.sleep(0.01)
except KeyboardInterrupt:
s.close()

 Observen que en el argumento del método send() estoy concatenando las etiquetas de 'CH0' y 'CH1' para poder separar las señales después con LabVIEW. (recuerden que '\n' es una secuencia de escape que indica nueva linea).

El diagrama a bloques del cliente para LabVIEW es bastante sencillo:

(Click para agrandar)

En el bloque TCP Open Conection colocamos la dirección IP* de la Raspberry y el puerto que establecimos para el socket en el programa de Python. Después de hacer la lectura utilizamos los bloques Match Patern para identificar las etiquetas de los canales y separar sus valores. Finalemnte se convierten las cadenas aisladas a numeros de tipo double para enviarlos al graficador. El sistema funcionando se ve así:

(Descargar archivos)

*Aquí estoy usando una ip local. En la práctica se requerirá utilizar una IP pública. Hay varias formas de hacerlo: #1, #2 , #3

Referencias:
Python Advanced Tutorial 6, DrapsTV
Raspberry Pi + Arduino + LabVIEW, Alfredo Cruz

Modelo lineal de un motor DC

El modelo lineal de un motor eléctrico de DC consiste en 2 ecuaciones diferenciales acopladas: el modelo eléctrico y el modelo mecánico. Debido a que ambos modelos están relacionados podemos escribir un modelo general el cual nos permitirá obtener una función de transferencia. En este ejemplo encontraremos la función de transferencia que relacione voltaje (entrada) con posición angular (salida). Primero, observemos el diagrama del circuito equivalente de la armadura del motor y el diagrama de cuerpo libre de rotor:

Para obtener la ecuación diferencial para el modelo eléctrico consideramos la ley de voltajes de Kirchoff:

 La fuerza contra-electromotriz se genera al iniciar el movimiento del rotor debido a que el campo magnético fijo del estator induce un voltaje en el devanado de la armadura (este voltaje es negativo con respecto al voltaje de entrada). La fcem es proporcional a la velocidad angular del rotor por lo que la constante Kb puede determinarse experimentalmente graficando el voltaje en las terminales del motor contra la velocidad angular del rotor (se verá que la relación no es realmente lineal en un intervalo grande pero puede usarse sólo la región lineal como una aproximación para el modelo).

Para obtener el modelo mecánico consideramos la segunda ley de Newton para movimiento angular:

Dónde J y b son el momento de inercia del rotor y el coeficiente de amortiguamiento por fricción respectivamente. Vemos que el torque es proporcional a la corriente en el motor. Asumiendo que no hay perdidas electromagnéticas, por fricción ni por calor, la potencia mecánica en el rotor debe igualar a la potencia eléctrica de manera que :

Por lo que el valor para el factor de fcem encontrado experimentalmente puede ser usado como valor para Kt. Observamos ahora que las ecuaciones (1) y (2) están relacionadas por la función de corriente.Para facilitar la sustitución obtendremos primero la transformada de Laplace de ambas ecuaciones:

Despejando I(s) de (4) y sustituyendo en (3) obtendremos el modelo unificado para el motor DC. Teniendo ya una sola ecuación obtenemos la función de transferencia posición/voltaje:

Ejecutando el siguiente código en Matlab asignando algunos valores a las constantes obtenemos su respuesta al escalón unitario en lazo cerrado:

J = 3.2284E-6;
b = 3.5077E-6;
K = 0.0274;
R = 4;
L = 2.75E-6;
s = tf('s');
P_motor = K/(s*((J*s+b)*(L*s+R)+K^2))
sys_cl = feedback(P_motor,1)
step(sys_cl)

Nota: Esto es lo que se conoce como identificación de sistema de caja blanca. Si no conocen ningún parámetro físico del motor, pueden usar el procedimiento de identificación de caja negra que describo aquí.

Referencias:
DC Motor Position: System Modeling, University of Michigan
DC Motor, MathWorks
DC Motor, How It Works? [video subtitulado en español], Learn Engineering