Tutorial viejito: Un Par de Trucos al hacer Guis en Matlab 7

Este tutorial data de mis años de estudiante hace casi dos décadas atrás, enjoy!

Matlab 7 permite poner imágenes en las ventanas hechas con la herramienta Guide, aunque no directamente pues no posee un control Picture o Image como en Visual Basic o C++ Builder. Lo malo es que los archivos de Ayuda de Matlab no explican cómo poner una imagen a una ventana (y si lo hacen yo no pude encontrar dónde), después de una tarde de examinar y diseccionar el Matlab descubrí la manera de hacerlo.

Y aquí viene lo que quiero explicar:

 
CÓMO PONER UNA IMAGEN EN UNA INTERFAZ GRÁFICA

Primero la imagen que queremos colocar en la ventana debe estar en formato BMP y en la misma carpeta donde se guardarán los archivos de la Interfaz Gráfica (el *.fig y el *.m).

Luego se debe crear una ventana en blanco, esto se hace escribiendo "Guide" en la ventana de comandos y presionando enter, entonces se abrirá la siguiente ventana:

Elegimos "Blank GUI" con lo que se abrirá esta ventana:

Después se jala un control Axes (el botón con la curva azul):

Mi ventana la he guardado como "unt1.fig" lo cual genera un archivo "unt1.m", ahora nos vamos para allá. Lo que debemos buscar es la función OpeningFcn que se activa apenas se abre la ventana:

Para el Matlab la ventana tiene por defecto el nombre de "figure1" y el Axes el de "axes1". Estos nombres pueden cambiarse, pero en este caso los vamos a usar así como están.

Lo que debe hacerse es añadir el siguiente código (no hay que hacerle caso a lo que está en verde, son comentarios y pueden borrarse):

Lo que hace este código es cargar la imagen (cuyo nombre es "dsp.bmp") a una matriz con la función "imread", por precaución se cambia a formato entero de 8 bits (aunque una imagen en BMP ya está en este formato).

La tercera línea hace que la imagen dependa de axes1. Las siguientes dos líneas son para configurar axes1: la propiedad "Visible" debe estar en "off" o de otro modo se verán los ejes coordenados a los costados de la imagen.

Lo demás es para que la imagen tenga el mismo tamaño que axes1. Los tres puntos simplemente son para indicarle al Matlab que la función continúa en la siguiente línea.

Al correr nuestra ventana debe verse así:

El ancho y alto de la imagen puede variarse variando el ancho y alto de axes1.

 

CÓMO HACER LA PARTE MÁS OSCURA PARA EL TÍTULO:

Para esto se deben jalar un control tipo "panel" y uno tipo "static text". A ambos se les da click derecho y se escoge "Send to back", lo cual los pondrá detrás de la imagen y no la taparán:

Al hacer doble click sobre cualquiera de ellos aparecerá la ventana de propiedades con las respectivas propiedades que pueden usarse para cada control. Por ejemplo, el panel y el texto estático tienen una propiedad llamada "BackgroundColor", la cual no tiene figure1, pero ésta posee la propiedad "Color". Es importante saber esto, pues no todos los controles tienen las mismas propiedades.

El panel tiene por defecto el nombre de "uipanel1" y el texto estático el de "text1".

Nos volvemos a ir a la función OpeningFcn y añadimos (sí, borré los comentarios del Matlab 7):

El color de la ventana por defecto es el que le da el sistema operativo. Si queremos que la ventana tenga el color que le asigna el windows y no uno fijo, y que el recuadro donde irá el título sea un poco más oscuro, es necesario tener los valores del color asignado (son 3, para rojo, verde y azul). Estos valores se cargan a un vector llamado "cc".

Para oscurecerlo un poco se le resta un número entre 0 y 1 al vector. El valor resultante se asigna a la propiedad "BackgroundColor" de text1 y uipanel1.

El resultado es éste:

Y con otras propiedades de pantalla (en Windows 2000, este es un tutorial MUY viejo):

Se le pueden arreglar algunas propiedades, como el tamaño y color de la letra de text1 y borrar el título de uipanel1 si se desea:

(no podía dejar mi sistema con ese amarillo tan molesto!)

Filtrado E Histograma De Una Imagen Para Matlab 6.5 ó 7

Otra entrada migrada de la vieja web :)


Una imagen digitalizada no es más que una matriz de números, donde cada número representa el valor de un píxel. Se diferencia de una señal digitalizada en que una señal es un vector, o varios vectores separados (como el caso de la música estéreo, compuesta por dos vectores de números que representan los canales derecho e izquierdo).

Para el procesamiento de una imagen existen varias técnicas, las cuales se usan para obtener información no visible, resaltar bordes, suavizarlos, quitar detalles irrelevantes, agrupar objetos, eliminar el fondo, etc.

Una técnica muy usada es el filtrado, el cual puede ser lineal o no lineal. El filtrado lineal es simplemente convolucionar la imagen con una matriz predefinida. (la convolución es una operación de sumas y multiplicaciones que se usa tanto en señales como en imágenes). Un filtro lineal es el "Filtro de Promedio" o "Mean Filter". Se usan para suavizar, detectar o resaltar bordes, eliminación de ruido, etc.

El filtrado no lineal se compone de técnicas más complicadas. El más básico es el Filtrado de Mediana o "Median Filter", que consiste en coger una pequeña porción de la imagen (generalmente 3x3 píxeles), evaluar el valor de los píxeles y cambiar sus valores al de aquél que tenga un valor central. Es útil para eliminar ruido impulsivo.

Otra técnica muy usada es la "Ecualización del Histograma". El Histograma de una imagen es el ploteo de los valores de sus píxeles. Una imagen en blanco tendrá todos sus valores iguales a 255, si la mitad es negra, en la gráfica del histograma aparecerán dos líneas iguales a ambos extremos: en los valores correpondientes al 0 y al 255. Una imagen de escala de grises tendrá en su histograma "x" píxeles con el valor 0, "y" píxeles con el valor 1, etc.

Así, el histograma es la representación de la densidad de probabilidad de cada valor de gris para esa imagen.

Tanto el histograma, como el histograma ecualizado, son vectores.

Ecualizar el Histograma es hacerlo lo más llano y separado posible. Esto hace que los píxeles se distribuyan más ampliamente por todo el rango de valores (del 0 al 255) y que en la imagen ecualizada se resaltarán detalles que antes no eran evidentes.

Para generar una imagen con el histograma ecualizado se requieren varios pasos:

1. Calcular el histograma de la imagen

2. Normalizar el histograma (dividirlo entre el número total de píxeles)

3. Calcular el histograma acumulado (ir sumando los píxeles desde el valor 0 al 255, esto originará una gráfica creciente)

4. Se aplica el algoritmo (como se trata de matrices, debe estar dentro de dos bucles anidados):

Para valores de la imagen original distintos de cero:
imagen_ecualizada(i,j) = histograma_acumulado(imagen_original(i,j))

Para valores de la imagen original iguales a cero:
imagen_ecualizada(i,j) = histograma_acumulado(imagen_original(i,j)+1)

Donde lo que va entre paréntesis es el índice, o índices, de cada píxel de la imagen o de cada valor del histograma. Así  "imagen_original(i,j)" viene a indicar el índice del vector del histograma acumulado correspondiente. De esta manera se le asigna a cada píxel de la imagen ecualizada (o imagen con el histograma ecualizado) la densidad de probabilidad acumulada correspondiente al valor del píxel de la imagen original.
Como este algoritmo está hecho para Matlab, y como Matlab no maneja índices iguales a cero, se considera que el histograma va de 1 a 256, en lugar de 0 a 255. Así la densidad de probabilidad del valor cero será la que está en el índice 1 en el vector del histograma.

El algoritmo que presento akí está considerado para matrices con valores tanto positivos como negativos. Es un archivo zip que incluye ejemplos de filtraje, ecualización y eliminación de ruido.

Link de donde aprendí la teoría para poder hacer esto:

http://homepages.inf.ed.ac.uk/rbf/HIPR2/filtops.htm

Matlab: Resolviendo la ecuación de onda bidimensional en el caso de una membrana de forma piramidal

El código fuente para resolver una ecuación de onda bidimensional está dado en el archivo twodimwavedirbc.m disponible en esta web: MATH 5343 Numerical Solution for Partial Differential Equations, Spring 2011.

La idea es resolver el ejercicio 12a del capítulo 12.2 del libro "Introduction to Numerical Ordinary and Partial Differential Equations Using MATLAB" de Alexander Stanoyevitch (de este libro proviene el scrip twodimwavedirbc.m):

Es decir, es una ecuación de onda cuya forma inicial es una pirámide. La dificultad aquí es dibujar la pirámide, yo la deduje como explico a continuacion:

La pirámide mostrada:

está formada por cuatro planos limitados por las rectas Y=X y Y = -X/2+4:

Las regiones donde Z es una función de X o de Y (dependiendo del plano donde se proyecta el eje Z) se muestran a continuación:

Para X<2, Y<2 y la región donde X<Y, la variable Z está definida por la función Z= Y/2.
Esto se hace para cada región de cada uno de los cuatro planos que definen la pirámide.

El script resultante es:

function z = pyramid(x, y)
z=0;

if y<=2 && x <= 2
   if y < x
       z=y/2;
   else
       z=x/2;
   end
end

if y<=2 && x > 2 && x <= 4
   if y < 4-x
       z=y/2;
   else
       z=-x/2+2;
   end
end

if y > 2 && x <= 2 && y <= 4
   if y < 4-x
       z=x/2;
   else
       z=-y/2+2;
   end
end

if y > 2 && y <= 4 && x > 2 && x <= 4
   if y < x
       z=-x/2+2;
   else
       z=-y/2+2;
   end
end

end 

Modifiqué un poco el archivo twodimwavedirbc.m para hacerlo más flexible:

function [x, y, t, U] = twodimwavedirbc2(phi, nu, a1, a, b1, b, T, h, c)
% a1: punto inicial del eje x
% b1: punto inicial del eje y

% solves the two-dimensional wave problem u_tt = c^2(u_xx+u_yy)
% on the rectangle {a1<= x <= a,  b1<= y <= b}, with u(x,y)=0
% on the boundary of the rectangle.

% variables de entrada: 
% phi=phi(x,y) = función de la forma inicial de la onda
% nu=nu(x,y) = función inicial de la velocidad onda

% a= right endpoint of x, b = right endpoint of y, 
% T = final time solution will be computed.  h=common gap on
% x, y-grids,   k= t-grid gap, c = speed of wave.
% wave.
% Output variables: x = row vector for first space variable, y =
% row vector for second space variable,  t = time grid row vector
%(starts at t=0, ends at t=T, has Nt equally spaced values), 
% U = (Nx)by(Ny)by(Nt) matrix of solution approximations at
% corresponding grid points (where Ny = number of y-grid points)
% x grid will correspond to first entries of U, y
% grid values to second entries of U, and t grid to third entries of U.
% CAUTION:  This method will only work if h is chosen so that the x
% and y grids can have a common gap size, i.e., if h = a/h and
% b/h must be integers.
% The time grid gap is chosen so that mu = 1/2, this guarantees the
% Courant-Friedrichs-Levy condition holds and simplifies the
% main finite difference formula.

k=h/sqrt(2)/c; %k is determined from mu = 1/2

if (b/h-floor(round(b/h))>10*eps)||(a/h-floor(round(a/h))>10*eps)
    fprintf('Space grid gap h must divide evenly into both a and b \r')
    fprintf(' Either try another input or modify the algorithm')
    error('M-file will exit')
end

Nx = round((a-a1)/h)+1; %number of points on x-grid
Ny = round((b-b1)/h)+1; %number of points on y-grid
Nt = floor(T/k)+1; %number of points on t-grid
U=zeros(Nx, Ny, Nt);
x=a1:h:a;
y=b1:h:b;
t=0:k:T;
% Recall matrix indices must start at 1.  Thus the indices of the
% matrix will always be one more than the corresponding indices that
% were used in theoretical development.

% Note that by default, zero boundary values have been assigned to all
% grid points on the edges of the rectangle (and, for the time being,
% at all grid points).
%Assign initial time t=0 values and next step t=k values.
for i=2:(Nx-1)
     for j=2:(Ny-1)
    U(i,j,1)=feval(phi,x(i),y(j));
    U(i,j,2)=.25*(feval(phi,x(i-1),y(j))+...
      feval(phi,x(i+1),y(j))+feval(phi,x(i),y(j-1))+...
      feval(phi,x(i),y(j+1)))+ k*feval(nu,x(i),y(j));     
     end
end
%Assign values at interior grid points
for ell=3:Nt %letter ell looks too much like number one
    U(2:(Nx-1),2:(Ny-1), ell) = ...
     +.5*(U(3:Nx,2:(Ny-1), ell-1)+ U(1:(Nx-2),2:(Ny-1), ell-1)...
     + U(2:(Nx-1),3:Ny, ell-1) + U(2:(Nx-1),1:(Ny-2), ell-1))...
     - U(2:(Nx-1),2:(Ny-1), ell-2);
end

Y el script principal es:

clc; clear all; close all;
phi = @(x, y)sin(2*x) *sin(2*y);
nu=@(x, y)0;
[x, y, t, U] = twodimwavedirbc2('pyramid', nu, 0, 4, 0, 4, 4, 0.1, 1);
[a b c] = size(U);
%ans = 26 26 46
for ell=1:c
     surf(x,y,U(:, :,ell));
     axis([0, 4, 0, 4 -1, 1]);
     M(:,:,ell) = getframe;
end
movie(M, 2, 4)


Lo he probado en Matlab, y acá está corriendo lo mejor que puede en Octave 4.4.1 a pesar de que me dice que la función Movie no está implementada:

Resolviendo una ecuación parabólica unidimensional en Matlab/Octave

El siguiente ejercicio se encuentra en el libro "Computational Partial Differential Equations Using MATLAB" by Yi-Tung Chen, Jichun Li, capítulo 2, ejercicio 3:

Una pista de cómo resolver este ejercicio está en la sección 2.21:

Los autores dan un código base sobre el qué trabajar. Luego de recibir un poco de ayuda de la gentita de Stackoverflow, este es el código para resolver el problema:


 %---------------------------------------------------------------
% para1d.m:
%    use the explicit scheme to solve the parabolic equation
%    u_t(x,t) = u_{xx}(x,t),          xl < x < xr, 0 < t < tf
%    u(x,0) = f(x),                   xl < x < xr
%    u_x(0,t) = gl(t), u_x(1,t) = gr(t),  0  < t < tf
%---------------------------------------------------------------
clc;
close all;
clear all;                  % clear all variables in memory

xl=0; xr=1;                 % x domain [xl,xr]
J = 31;                     % J: number of division for x
dx = (xr-xl) / J;           % dx: mesh size
tf = 0.1;                   % final simulation time
Nt = 51;                    % Nt: number of time steps
dt = tf/Nt/4;

mu = dt/(dx)^2;

if mu > 0.5         % make sure dt satisy stability condition
    error('mu should < 0.5!')
end

% Evaluate the initial conditions
x = xl : dx : xr;              % generate the grid point
f = x;

% store the solution at all grid points for all time steps
u = zeros(J+1,Nt);
u_ex = zeros(J+1,Nt);

% boundary condition at left side
gl = 0;
% boundary condition at right side
gr = 0;

% Find the approximate solution at each time step
for n = 1:Nt
    t = n*dt;         % current time
     
    if n==1
      for j=2:J
       % first time step
        u(j,1) = f(j) + mu*(f(j+1)-2*f(j)+f(j-1));
      end
    end
  
    if n > 1
        for j=2:J
         % interior nodes
           u(j,n) = u(j,n-1) + mu*(u(j+1,n-1) - 2*u(j,n-1) + u(j-1,n-1));
        end
    end
  
    % hint 2.21 Los valores en la frontera es mejor calcularlos al final
    % pues estoy cogiendo índices que dependen de los nodos interiores
    u(1, n) = u(2, n) - dx * gl;
    u(J+1,n) = u(J, n) + dx * gr;
  
    % calculate the analytic solution
    for j=1:J+1
      xj = xl + (j-1)*dx;    
      suma = zeros(1000 , 1);

      for k= 1:1000
        suma(k) = 4/(((2*k-1)*pi)^2);
        suma(k) = suma(k) * exp(-((2*k-1)^2)*t*pi^2) * cos((2*k-1)*pi*xj);
      end

      u_ex(j, n)= double(0.5 - sum(suma));
    end
end

% Plot the results
tt = dt : dt : Nt*dt;
figure(1)
colormap(cool);     % draw gray figure
surf(x, tt, u');     % 3-D surface plot
xlabel('x')
ylabel('t')
zlabel('u')
title('Numerical solution of 1-D parabolic equation')

figure(2)
surf(x, tt, u_ex');     % 3-D surface plot
xlabel('x')
ylabel('t')
zlabel('u')
title('Analytic solution of 1-D parabolic equation')

maxerr=max(max(abs(u-u_ex))),


El resultado es:

Probado en Octave 4.4.1 y Matlab 2008 R2.

Algoritmos para calcular la factorización QR de una matriz

La factorización QR se utiliza para descomponer una matriz (que puede representar un sistema de ecuaciones a resolver) en dos matrices, con las cuales es más fácil de operar y hallar soluciones al sistema de ecuaciones.
Más información en la WIkipedia.

Los siguientes códigos están basados en los algoritmos mostrados en el libro "Numerical Linear Algebra" de Trefethen y Bau.

Algoritmo de Gram-Schmidth Modificado (devuelve la factorización QR reducida):

El código en Matlab es:

function [Q, R] = gs(A)

[m, n] = size(A);

R = zeros(n, n);

V = A;

Q=zeros(m, n);

for i =1:n

    R(i,i)= norm(V(:,i));

    Q(:,i)= V(:,i)/R(i,i);

   

    for j=i+1:n

       R(i,j)= (Q(:,i)')*V(:,j);

       V(:,j)=V(:,j) - R(i,j)*Q(:,i);

    end

end

"A" es una matriz de rango completo, la función devuelve dos matrices: Q (la matriz ortonormal) y R (la matriz triangular).

Problemas con este algoritmo: Sólo da la factorización QR reducida, y sólo funciona para matrices de rango completo (es decir: ninguna columna es múltiplo de, y/o resultado de, sumar otras columnas de la mariz entre sí).

Modifiqué este algoritmo para que funcionara también con matrices de "rango deficiente" (lo opuesto a "rango completo"), y que devolviera la factorización completa; esto último se hace completando las columnas de la matriz Q, las cuales son ortogonales entre sí. En palabras del curso de Álgebra Lineal: "completando la base ortonormal del espacio vectorial generado por las columnas de Q". Dicho de otra forma: aplicando el Teorema de la Base Ortonormal Incompleta.

El código en Matlab es:

function [Q, R] = gsC(A)

clc;

[m, n] = size(A);

R = zeros(m, n);

V = A;

Q=zeros(m, n);

k=1; % inicializo eta variable, sólo la uso para matrices de rango deficiente

for i =1:n

    R(i,i)= norm(V(:,i));  

   

    if abs(R(i,i)) > 1.0e-10

        Q(:,i)= V(:,i)/R(i,i);

    else  % salió cero! La matriz A es de rango deficiente

        r= zeros(1, m); % este va a ser el vector que reemplazará a la fila problemática de A

              

        Q(:,i)= r;

        Q(k, i)=1;

        k=k+1;

        Q2 = Q; % Q2 tendrá los valores de Q, en Q irá la nueva base ortonormal

        % ortogonalizo otra vez aplicando Gram-Schmidt a la matriz Q para

        % hallar un vector con el qué reemplazar a la fila problemática y

        % seguir teniendo una base ortonormal

        for j = 1 : i-1

            r(j) = (Q2(:,j)')*Q2(:,i);

            Q(:,i)=Q(:,i) - r(j)*Q2(:,j);

            Q(:,i)=Q(:,i)/norm(Q(:,i));

        end

    end

   

    for j=i+1:n

       R(i,j)= (Q(:,i)')*V(:,j);

       V(:,j)=V(:,j) - R(i,j)*Q(:,i);

    end

end

% hasta aquí Q y R son la descomposición reducida

% acá hallo la completa, Gram Schmidt otra vez

if m > n   

    I = eye(m, m - n);

    Q=[Q I];

    V=Q;

    r= zeros(1, m);

   

    k=0;

    for i=n+1:m

        for j = 1 : n + k

            r(j) = (V(:,j)')*V(:,i);

            Q(:,i)=Q(:,i) - r(j)*V(:,j);

            Q(:,i)=Q(:,i)/norm(Q(:,i));

        end

        k=k+1;

    end

end

"A" es una matriz cualquiera, la función devuelve dos matrices: Q (la matriz ortonormal) y R (la matriz triangular). 

Otra forma de calcular la factorización QR (y la que utiliza internamente Matlab) es usando la factorización de Householder:

El código en Matlab es:

function [Q, R] = houseHolder(A)

clc;

[m n] = size(A);

R= A;

Q = eye(m,m);

n1=min(m,n);

for i= 1:n1

    x = R(i:m, i);

    s=0;

   

    if x(1)==0

        s = 1;

    else

        s= sign(x(1));

    end

   

    e = zeros(m - i + 1, 1);   

    e(1)=1;   

    u = s * norm(x) * e + x;

    u = u / norm(u);

    R(i:m, i:n) =  R(i:m, i:n) - 2*(u*u')*R(i:m, i:n);

    Q(:,i:m) = Q(:,i:m) - Q(:,i:m)*(2*u*u');   

end

Igual que en los casos anteriores, "A" es una matriz cualquiera, la función devuelve dos matrices: Q (la matriz ortonormal) y R (la matriz triangular).

Los tres algoritmos han sido probados en Matlab 2008a, así que deben correr en cualquier versión posterior, y no tan posterior :)

Diferencias entre las versiones del Matlab

Matlab 6.5: Puede instalarse en Windows 98 SE y en Windows 2000 y en computadoras con 128Mb de RAM y procesador de 500MHz.

Matlab 7: Puede instalarse en Windows 2000 y en Windows XP y en computadoras con 512Mb de RAM y procesador de un núcleo de 2GHz. No corre en computadoras con procesadores de doble núcleo. Tampoco ejecuta los programas hechos con el Guide de Matlab 6.5 ni con Simulink del Matlab 6.5

Matlab 7.6: Puede instalarse en Windows XP y Windows Vista. Corre en computadoras con procesador de doble núcleo. Puede ejecutar los programas hechos con el Guide del Matlab 7. No es recomendable correrlo en computadoras con menos de 1Gb de RAM.

Página oficial de Matlab: http://www.mathworks.com/

Alternativas gratuitas:
Octave
FreeMat
Scilab