API de procesado de audio para JAVA

A la hora de realizar un prototipo en JAVA sobre un proyecto que me encuentro realizando, me he encontrado con que no hay mucho material disponible escrito para dicho lenguaje de fácil uso para el procesado de señal. Por ello, yo mismo he tenido que programar métodos para realizar funciones que necesitaba de procesado. Mi proyecto versa sobre el tratamiento del audio por lo que, básicamente, el código escrito está orientado a ello. Muchas de las funciones son muy concretas, pues no quise complicarme escribiendo algo muy genérico pero, otras tantas, sí se pueden usar de forma genérica (también está la opción de modificar según la necesidad).

Podéis descargar la API como fichero con extensión .java haciendo clic aquí.

A continuación presento una descripción breve sobre las funciones incluidas, ya que la clase escrita no está documentada.

- Suavizado de Pitch

La función admite como entrada un vector de números reales en doble precisión representantes de un contorno de frecuencia de pitch en el tiempo (también podría usarse con otros propósitos). Como salida devuelve el mismo contorno suavizado. Este suavizado se basa en umbrales heurísticos supuesto que el contorno de pitch es el logaritmo en base e de la secuencia de frecuencias de pitch en hertzios. La realización de este método se encuentra motivada por el hecho de que, en muchas ocasiones, durante la detección del pitch a partir del espectro o funciones de autocorrelación, existen armónicos más energéticos que el tono fundamental, lo que se traduce en picos erróneos dentro de la secuencia. Se considera errónea una ráfaga de valores cuando dicha ráfaga es inferior a 10 cifras consecutivas y la pendiente generada en su entorno es superior a 0.5.

- Cálculo de la PSD

Calcula la PSD (módulo de la FFT pero no en dB) entre $0$ y $2\pi$ de una secuencia de muestras dada como un vector de números reales en doble precisión. La secuencia de salida tiene la misma longitud y se calcula según:

$Y(k) = \left | \sum_{n=0}^{N-1}x(n)e^{-i\frac{2\pi nk}{N}} \right |$

- Filtro de Mediana de Orden 5

Muy útil para la eliminación de ruido aditivo impulsivo. Realmente puede modificarse en su interior el orden, siempre y cuando este sea impar. La muestra k-ésima filtrada se obtiene a partir del cálculo de la mediana como:

$y(k) = mediana(x(k - (n-1)/2),...,x(k + (n+1)/2))$

- Cuantizador Musical

Aproxima una secuencia de valores de pitch frecuencial dados en forma logarítmica a semitonos de la escala musical. Esto sólo se encuentra habilitado en el rango frecuencial que va del Sol2 al Sol5.

- Cálculo de la Energía

Calcula la energía de una secuencia de muestras dadas en un vector de datos reales de doble precisión. Dicha energía se obtiene de la forma:

$E(dB) = 10log_{10}\sum_{n=0}^{N-1}x(n)^{2}$

- Interpolador de Orden 5

Este interpolador aumenta la resolución en un factor 5 de un conjunto de muestras dadas en el formato ya mentado. Este orden de interpolación es inamovible ya que el funcionamiento del método se basa en la inserción de ceros y en la aplicación de un filtro FIR simétrico especial (cuyos coeficientes sólo son válidos para este orden de interpolación).

- Enventanador de Hamming

Enventana una secuencia de muestras con Hamming. La función de la ventana de Hamming es:

$w(n) = a_{0}-a_{1}\cos(\frac{2\pi n}{N-1})$

- Cálculo de la Media

Calcula la media de un vector de datos de entrada como:

$\mu = \frac{1}{N}\sum_{n=0}^{N-1}x(n)$

- Cálculo de la Varianza

Calcula la varianza de un vector de datos de entrada como:

$\sigma^2 = \frac{1}{N-1}\sum_{n=0}^{N-1}(x(n)-\mu)^2$

- Normalizador

Normaliza un vector de muestras de entrada por el máximo del vector en valor absoluto.

- Filtrado

Método para el filtrado de un vector de datos reales en doble precisión de entrada. Devuelve otro vector del mismo tamaño al que se le ha aplicado un filtrado dados los coeficientes del correspondiente filtro FIR o IIR. Dado un filtro de orden (p,q), la salida n-ésima filtrada se obtiene como:

$y(n) = b_{0}x(n)+...+b_{q}x(n-q) - a_{1}y(n-1)-...-a_{p}y(n-p)$

- Cálculo de la Autocorrelación

Dada una secuencia de muestras de entrada, este método calcula su autocorrelación haciendo uso del estimador sesgado, el cual es de uso preferente debido a su inferior varianza, lo que se traduce en una mejor convergencia de las muestras de autocorrelación lejanas. Se devuelve un vector de salida de la misma dimensión que el de entrada, el cual contiene las muestras correspondientes a la parte positiva del eje de abscisas. El estimador sesgado de la autocorrelación para la parte positiva del eje de abscisas se define como:

$r_{x}(k) = \frac{1}{N}\sum_{n=0}^{N-1-k}x(n)x(n+k)$

- Análisis LPC

Se encarga de obtener hasta un total de 13 coeficientes del predictor lineal que modela, como un filtro todo-polos, la respuesta del tracto vocal. Estos coeficientes se obtienen a partir de la resolución de las ecuaciones de Yule-Walker a partir del cálculo de muestras de autocorrelación dado un vector de muestras reales en doble precisión de entrada.

- Detector de Pico Máximo

Esta función detecta la presencia de un pico máximo absoluto en una secuencia de valores de entrada. Es la misma función que ya presenté en otra entrada escrita en MatLab.

- Remuestreador de Orden 5

No sólo es un remuestreador de orden 5 sino que, únicamente, funciona como remuestreador para pasar de una frecuencia de 8kHz a 1.6kHz. Esto es así ya que, al principio, se aplica un filtrado FIR de fase lineal cuyos coeficientes son estáticos y almacenados para este propósito. Una vez filtrado el vector de entrada, se selecciona una de cada 5 muestras.

Detector de pitch (II)

Vuelvo a presentar un nuevo detector de pitch, esta vez, basado en autocorrelación. Dicho detector está programado en MatLab y está basado en el algoritmo SIFT. Su diagrama de bloques se observa a continuación.

Primero, el audio de entrada se preprocesa remuestreando desde 44100Hz (asumimos dicha frecuencia de muestreo para el audio de entrada) a 8000Hz, habiendo ponderado previamente ambos canales estéreo para, finalmente, normalizar la señal.

Seguidamente filtramos la señal mediante un filtro tipo Butterworth paso-bajo con frecuencia de corte 800Hz (correspondiente, aproximadamente, a un Sol-5 en la escala musical) para, a continuación, aplicar un downsampling correspondiente a la nueva frecuencia máxima presente en la señal, lo que se traduce, según el criterio de Nyquist, en un remuestreo a 1600Hz. Fragmentamos en tramas el audio filtrado y, para cada trama, aplicamos el proceso descrito a partir de ahora. Se toma el frame y con él se realiza un estudio LPC (Linear Predictive Coding), con el fin de modelar la respuesta espectral referente al tracto vocal (formantes) mediante un filtro todo-polos. Una vez obtenidos los coeficientes de dicho filtro, se filtra la trama mediante otro filtro FIR cuyos coeficientes son los anteriores del modelo LPC. Con ello, intentamos eliminar la información propia del tracto vocal de la señal con el fin de mantener sólo la información de interés: el pitch. Por último se calcula la función autocorrelación de la señal filtrada y se interpola para ganar resolución temporal. El pico máximo detectado a la salida, se corresponderá con el período de pitch deseado.

La función que implementa el anterior proceso toma como datos de entrada el tamaño del frame (en segundos, con un valor típico de 40ms) y el solapamiento (también en segundos, con un valor típico de 10ms) aparte de, claro está, el propio audio cuyo pitch deseamos estimar.

También se incluye un detector de actividad de voz basado en hangover con los parámetros por defecto fijados en la función (cambiar si fuera necesario). La salida es la secuencia de pitch, a la que se le ha aplicado el logaritmo neperiano (si se desea obtener la frecuencia en Hz's, cambiar desde el código o aplicar la exponencial al vector de salida), de tal forma que no se incluyen las estimaciones para las tramas cuyo flag proporcionado por el detector de voz se corresponda con silencio.

La gráfica adjunta muestra un ejemplo del pitch estimado para la canción Hush de Deep Purple tarareada.

IMPORTANTE: La función detPico.m se encuentra en la entrada justamente anterior.

% Detector de pitch
% Autor: Iván López Espejo

function pitch = detPitch(x,wsize,sol)


% Preprocesado.
fm1 = 44100;
fm = 8000;
x = (x(:,1) + x(:,2))/2;
x = resample(x,fm,fm1)
x = x/max(abs(x));

% Parámetros del detector de actividad de voz.
THRESHOLD = 0.01;
MIN_SPEECH_FRAME_HANGOVER = 4;
HANGOVER = 1;
meanEn = 0;
nbSpeechFrame = 0;
ind = 1;

% Diseño de filtro paso-bajo Butterworth con frecuencia de corte a 800Hz.
fc = 800;
Nf = 2; % Orden del filtro.
Wn = 2*fc/fm;
[b,a] = butter(Nf,Wn,'low');
x = filter(b,a,x); % Filtramos paso-bajo.
x = resample(x,2*fc,fm); % Downsampling.
% Análisis por frame.
wsize = round(wsize * 2 * fc);
sol = round(sol * 2 * fc);
% Definimos las constantes de la ventana.
a0 = 0.53836;
a1 = 0.46164;
Npred = 12; % Número de coeficientes del predictor lineal.
b0 = 1;
fint = 5; % Factor de interpolación.
pitch = [];
for m = 1:(wsize-sol-1):length(x)
    if wsize > length(x) - m
        N = length(x) - m + 1;
    else
        N = wsize;
    end
    % Detección de la actividad de voz.
    frameEn = dot(x(m:m+N-1),x(m:m+N-1))/N;
    if frameEn - meanEn > THRESHOLD
        VAD(ind) = 1;
        nbSpeechFrame = nbSpeechFrame + 1;
        if nbSpeechFrame > MIN_SPEECH_FRAME_HANGOVER
            hangOver = HANGOVER;
        end
    else
        if hangOver == 0
            VAD(ind) = 0;
        else
            VAD(ind) = 1;
            hangOver = hangOver - 1;
        end
    end
    ind = ind + 1; % Incremento del contador.
    % Enventanado del frame mediante Hamming.
    n = 1:N;
    win = a0 - a1*cos(2*pi*n/(N-1));
    frame = x(m:m+N-1).*win';
    % Método de la autocorrelación para modelado de los coeficientes del
    % predictor lineal asociado a la trama.
    acoef = lpc(frame,Npred);
    % Usamos los anteriores coeficientes como coeficientes de un filtro
    % FIR.
    eps = filter(acoef,b0,frame);
    reps = xcorr(eps,'biased'); % Cómputo de la autocorrelación.
    % Interpolación para la obtención de mayor resolución temporal.
    Pm = interp(reps,fint);
    Pm = Pm(floor(length(Pm)/2):length(Pm));
    % Detectamos el pico máximo en la función de autocorrelación.
    picomax = detPico(Pm);
    pause();
    if VAD(ind-1) == 1
        pitch = [pitch (log(2*fc*fint/(picomax-1)))];
    end
end
% Representamos el resultado.
plot(pitch,'-*r')
grid on
xlabel('Evolución temporal')
ylabel('Frecuencia de pitch logarítmica')
title('PITCH estimado')

Detector de actividad de voz (VAD)

En esta entrada se recoge la implementación de un detector de actividad de voz a través de una red neuronal artificial entrenada mediante propagación hacia atrás.

Este tipo de redes neuronales es ideal para el reconocimiento de patrones, de modo que entrenamos el sistema con vectores de características representativos de pasajes de silencio (ruido de fondo) y pasajes de voz. Para ello caracterizamos una trama de audio mediante 12 coeficientes cepstrales MFCC (los primeros 13 coeficientes salvo el primero, ya que este está relacionado con la energía de la trama y en este caso no es representativo) y un coeficiente de entropía, computada como:

 $H = -\sum_{k=1}^N p_{k}\log p_{k}$

Donde $N$ es el número de muestras del módulo del espectro de la trama en consideración y $p_{k}$ representa el valor k-ésimo de la función densidad de probabilidad calculada a partir de dicho espectro, definiéndose como:

$p_{k} = \frac{S_{k}}{\sum_{j}^N S_{j}}$

Con $S$ el módulo del espectro de la trama.

La red neuronal que he construido para solventar el problema utiliza tramas de 128 muestras (donde consideramos la señal aproximadamente estacionaria) para un audio muestreado a 8kHz y está diseñada para usar 10 capas ocultas. Dicha red está entrenada con 20 vectores de características representativos de silencio y otros 22 vectores representativos de voz. El resto de parámetros se extrae del código MatLab incluido a continuación.

Obsérvese un ejemplo de un fragmento de audio al que se le eliminan los silencios mediante esta implementación de VAD. El funcionamiento es bastante bueno. El programa va consumiendo tramas de 128 muestras a 8kHz, clasificándolas como voz o silencio gracias a la red neuronal previamente entrenada. Recomiendo sacar del código la parte del entrenamiento de la red neuronal, a fin de no estar inicializándola constantemente cada vez que se llame a la función.

 Audio de entrada.

Audio de salida con silencios eliminados.

Esta herramienta que he construido está pensada para funcionar previamente a la detección del pitch, a fin de poder eliminar las partes sin voz para no obtener falsas estimaciones.

Se puede ver el código a continuación, así como descargar mi fichero con patrones de entrenamiento haciendo clic aquí.

% Detector de actividad de voz basado en red neuronal
% Autor: Iván López Espejo.

function [y,sil] = vadet(x,fm)

% Remuestreo a 8kHz.
x = resample(x,8000,fm);
x = x/max(abs(x));
% Definición de algunos parámetros.
N = 128;
% Inicialización y entrenamiento de la red neuronal.
nump = 2; % Número de clases.
% Patrones de entrenamiento.
load Patrones
P = Patrones;
% Salidas.
T = zeros(2,42);
T(1,21:42) = ones(1,22);
T(2,1:20) = ones(1,20);
S1 = 10; % Número de capas ocultas.
S2 = nump; % Número de capas de salida (= al número de clases).
[R,Q] = size(P);
epochs = 10000; % Número de iteraciones.
goal_err = 5e-5; % Error objetivo.
a = 0.7; % Define el rango de variables aleatorias.
b = -0.7;
W1 = a + (b-a) * rand(S1,R); % Pesos entre la entrada y las neuronas ocultas.
W2 = a + (b-a) * rand(S2,S1); % Pesos entre las neuronas ocultas y las de salida.
b1 = a + (b-a) * rand(S1,1); % Pesos entre la entrada y las neuronas ocultas.
b2 = a + (b-a) * rand(S2,1); % Pesos entre las neuronas ocultas y las de salida.
n1 = W1 * P;
A1 = logsig(n1);
n2 = W2 * A1;
A2 = logsig(n2);
e = A2 - T;
error = 0.5 * mean(mean(e.*e));   
nntwarn off
for itr = 1:epochs
    if error <= goal_err
        break;
    else
         for i = 1:Q
            df1 = dlogsig(n1,A1(:,i));
            df2 = dlogsig(n2,A2(:,i));
            s2 = -2 * diag(df2) * e(:,i);                  
            s1 = diag(df1)* W2' * s2;
            W2 = W2 - 0.1 * s2 * A1(:,i)';
            b2 = b2 - 0.1 * s2;
            W1 = W1 - 0.1 * s1 * P(:,i)';
            b1 = b1 - 0.1 * s1;

            A1(:,i) = logsig(W1*P(:,i),b1);
            A2(:,i) = logsig(W2*A1(:,i),b2);
         end
            e = T - A2;
            error = 0.5 * mean(mean(e.*e));
            disp(sprintf('Iteración: %5d mse: %12.6f%',itr,error));
            mse(itr) = error;
    end
end
threshold = 0.9; % Umbral del sistema (umbral más alto = mayor precisión).
TrnOutput = real(A2>threshold);
disp('Resultado del entrenamiento: ')
disp(TrnOutput)
% Inicialización del audio de salida.
y = [];
sil = [];
for j = 1:N:length(x)-N
   coefs = calCoeficientes(x(j:j+N-1))';
   % Clasificamos el frame en voz o silencio.
   n1 = W1 * coefs;
   A1 = logsig(n1);
   n2 = W2 * A1;
   A2test = logsig(n2);
   resultado = real(A2test>threshold);
   if resultado(1) == 1 && resultado(2) == 0
       y = [y; x(j:j+N-1)];
   else
       sil = [sil; x(j:j+N-1)];
   end
end

Detector de pitch

Como pequeña parte de mi Proyecto Fin de Carrera, en el cual ya estoy sumergido y que desvelaré más adelante, he requerido programar en MatLab un detector del pitch para una secuencia de audio de entrada. Con esta herramienta pretendo modelar la secuencia de notas musicales desconocidas contenidas en un audio de entrada.

El detector de pitch programado en MatLab se basa en la idea de máxima semejanza espectral. 

Se construye en primer lugar una base de datos de espectros ideales de las distintas notas plausibles de detectar mediante convolución de una ventana espectral con un tren de impulsos que simboliza las posiciones del tono fundamental y sus armónicos.

El usuario entonces podrá controlar como parámetros el tamaño de la ventana y el solapamiento con el entorno, determinando así el tamaño de los fragmentos temporales que aproximaremos como estacionarios en el audio de entrada. A cada uno de estos fragmentos enventanados (mediante hamming por defecto) se le calcula el módulo de la FFT y se realiza un producto escalar de vectores entre dicho módulo y cada uno de los espectros ideales de la base de datos. El espectro ideal, tal que multiplicado escalarmente por el módulo de la FFT del fragmento enventanado del audio de entrada devuelva el máximo valor, será finalmente el correspondiente con la nota detectada, habiéndose minimizado así el error cuadrático:

$E = ||H(\omega)-G(\omega)||^{2} = ||H(\omega)||^{2} + ||G(\omega)||^{2}-2H(\omega)G(\omega)$

Donde $H(\omega)$ es el espectro ideal y $G(\omega)$ es el espectro real.

Se han realizado distintas pruebas y los resultados óptimos se obtienen utilizando un único armónico aparte del tono fundamental en la construcción del espectro ideal. También se han fijado otros parámetros como un umbral energético para establecer si existe nota musical o el fragmento se corresponde con silencio.

El código para MatLab se encuentra a continuación. El éxito en las pruebas realizadas es del 100% y su funcionamiento óptimo. La ayuda se encuentra detallada en la propia función.

Nota: He limitado la detección de notas en el rango Sol2 - Sol5 (3 octavas), correspondiente a mi registro vocal (ya se explicará el porqué más adelante).

% DETECTOR DE PITCH
%
% Autor: Iván López Espejo
%
% y = detectorPitch(x,wsize,sol,fm)
%
% Donde 'x' es el audio mono a analizar de entrada, 'wsize' el tamaño de la
% ventana o fragmentos en los que se segmentará el audio de entrada (valor
% típico en el rango 50ms-250ms) dado en segundos, 'sol' el solapamiento
% entre fragmentos adyacentes también dado en segundos (valor típico en el
% rango 10ms-50ms) y 'fm' la frecuencia de muestreo del audio de entrada
% dada en Hz's.

function y = detectorPitch(x,wsize,sol,fm)

% Lo damos en término de muestras.
wsize = round( wsize * fm );
sol = round( sol * fm );

% Obtención de los espectros ideales de referencia.
% k = floor((fm/(2*freq)) - 1);
k = 1;
l = 1;
h = [0.2 0.55 0.8 0.95 1 0.95 0.8 0.55 0.2];
o = 2;
for n = 8:12
        freq = 440*(2^((o-4) + ((n-10)/12)));
        vec = zeros(1,fm/2);
        for arm = 1:(k+1)
            vec(round(freq*arm)) = 1;
        end
        vec = conv(vec,h);
        vec = vec(5:length(vec)-4);
        sideal(l,:) = vec;
        l = l + 1;
end
for o = 3:4
    for n = 1:12
        freq = 440*(2^((o-4) + ((n-10)/12)));
        vec = zeros(1,fm/2);
        for arm = 1:(k+1)
            vec(round(freq*arm)) = 1;
        end
        vec = conv(vec,h);
        vec = vec(5:length(vec)-4);
        sideal(l,:) = vec;
        l = l + 1;
    end
end
o = 5;
for n = 1:8
        freq = 440*(2^((o-4) + ((n-10)/12)));
        vec = zeros(1,fm/2);
        for arm = 1:(k+1)
            vec(round(freq*arm)) = 1;
        end
        vec = conv(vec,h);
        vec = vec(5:length(vec)-4);
        sideal(l,:) = vec;
        l = l + 1;
end

% Calculamos el espectro de cada trozo enventanando según hamming.
% Definimos las constantes de la ventana.
a0 = 0.53836;
a1 = 0.46164;
% Umbral de detección de nota.
umbral = 2;
l = 1;
for j = 1:(wsize-sol-1):length(x)
    if wsize > length(x) - j
        N = length(x) - j + 1;
    else
        N = wsize;
    end
    n = 1:N;
    vent = a0 - a1*cos(2*pi*n/(N-1));
    trozo = x(j:j+N-1).*vent';
    % Calculamos la FFT del trozo.
    spc = abs(fft(trozo,fm));
    spc = spc(1:fm/2);
    spc = spc/max(spc);
    for n = 1:1:length(sideal(:,1))
        vecAux(n) = sum(spc'.*sideal(n,:));
    end
    [maxim,pos] = max(vecAux);
    y(l) = 0;
    if maxim > umbral
        y(l) = pos;
    end
    l = l + 1;
end

% Bloque para devolver la secuencia de notas.
ind = 1;
for j = 1:length(y)-1
    if y(j) ~= y(j+1)
       y2(ind) = y(j);
       ind = ind + 1;
    end
end
y2(ind) = y(length(y));
vNotas = ['Sol2 ';'Sol#2';'La2  ';'La#2 ';'Si2  ';'Do3  ';'Do#3 ';'Re3  ';'Re#3 ';'Mi3  ';'Fa3  ';'Fa#3 ';'Sol3 ';'Sol#3';'La3  ';'La#3 ';'Si3  ';'Do4  ';'Do#4 ';'Re4  ';'Re#4 ';'Mi4  ';'Fa4  ';'Fa#4 ';'Sol4 ';'Sol#4';'La4  ';'La#4 ';'Si4  ';'Do5  ';'Do#5 ';'Re5  ';'Re#5 ';'Mi5  ';'Fa5  ';'Fa#5 ';'Sol5 '];
ind = 1;
for j = 1:length(y2)
    if y2(j) ~= 0
        y3(ind,:) = vNotas(y2(j),:);
        ind = ind + 1;
    end
end
% Imprimimos la secuencia de notas detectada por pantalla.
y3

Corrección de ojos rojos

A veces, cuando realizamos una fotografía con flash, podemos apreciar el indeseado efecto de los ojos rojos. Esto es debido a que la pupila se encuentra, normalmente, muy dilatada en el momento de realizar la fotografía con flash (ya que esta tiene lugar en un ambiente oscuro y, por tanto, la pupila ha de captar mayor luminosidad para integrar una imagen del entorno correctamente). La luz del flash entonces entra a través de la pupila y se refleja produciendo ese característico color rojo propio de los capilares sanguíneos situados en dicha zona.

Flashes modernos, previo a la realización de la fotografía, emiten un primer destello cuya finalidad es que la pupila se contraiga. Con esta misma premisa, un truco casero previo a la realización de la foto, consistiría en someterse a una fuente intensa de luz durante unos instantes a fin de conseguir la contracción de la pupila.

Si no podemos evitar el efecto, algunas cámaras digitales y programas de retoque fotográfico incluyen alguna funcionalidad para mitigarlo. De esto se trata la presente entrada, en la que desarrollamos un sistema semi-automático para intentar corregir el efecto de los ojos rojos.

A continuación detallamos los pasos de la solución particular propuesta.

Esta solución se encuentra inspirada someramente por el artículo de Jalal acerca de la detección de ojos. Observé que, ante el color rojo intenso del efecto producido por el flash sobre los ojos, las componentes cromáticas, según el modelo YCbCr, presentaban unas características bien diferenciadas y es que, los ojos en la componente Cb tendían al negro y en la componente Cr al blanco. A partir de esta idea, tratamos de aislar y localizar los ojos en la imagen. Asistamos paso a paso la solución mediante un ejemplo. 

Consideremos la siguiente fotografía de partida:

Transformando al modelo de color YCbCr, obtenemos que los canales Cb y Cr son, respectivamente:

Si invertimos el canal Cr, los ojos quedarán también bastante más oscuros que el resto de la imagen, obteniendo:

Ahora, si sumamos Cb + Cr’, conseguiremos un contraste tal que los ojos serán mucho más oscuros que el resto de la imagen, de tal modo que obtenemos:

Aplicando una función de stretching sobre el anterior resultado conseguimos un mayor contraste entre los ojos y el fondo, obteniendo lo siguiente:

Justamente en este punto entra en juego el primer parámetro de caracterización de la función, y este es el nivel de umbralización de la anterior imagen con stretching. Para un valor umbral de 128, obtenemos:

Este parámetro es crítico y su nivel óptimo impredecible, pues es dependiente de la característica aleatoria de cada imagen. Si el valor es demasiado bajo, se tiende a perder potencial información sobre ojos en la imagen y, si es demasiado alto, pueden considerarse otras partes de la imagen como falsos ojos rojos, siendo ambas situaciones indeseables (podría aplicarse un algoritmo de umbralización óptimo como es el método de Otsu para la  estimación del umbral tal que se maximice la varianza inter-clase). En la mayoría de los casos, un valor intermedio en un rango con bastante tolerancia debe ser más que suficiente. En el caso que nos ocupa, como se observa, el valor intermedio comentado de 128 (umbral standard) nos proporciona un resultado bastante bueno, pues ya tenemos delimitados nuestros ojos rojos.

Lo único que nos queda ya por hacer es determinar su posición dentro de la imagen. Para tal fin, usaremos proyecciones de histogramas. Para mayor comodidad de análisis con los mismos, invertiremos de nuevo la imagen anterior, a fin de observar los ojos en blanco y quedar estos caracterizados en las proyecciones de histograma como lóbulos de un cierto ancho:

A continuación calculamos las proyecciones de histograma, obteniéndose para las proyecciones horizontal y vertical, respectivamente:

Se observa en ambas imágenes dos picos, los cuales caracterizan las coordenadas del punto central de cada uno de los ojos. Para no llevar a equívocos, pues ahora precisamos determinar la posición de los picos y el ancho de los lóbulos, suavizamos las proyecciones de histograma mediante una función gaussiana 1-D, en la cual entra en juego el segundo parámetro de la rutina, y este es la varianza de la función. Conseguiremos un mayor suavizado para una mayor varianza (la longitud de la función gaussiana con la que convolucionamos las proyecciones de histograma es constante y es de 11 muestras). Un valor que funciona bien es σ^2 = 16. De hecho, este parámetro es potencialmente extraíble de la función, ya que funciona bien para el valor constante antes comentado, aunque, no obstante, he decidido dejarlo. A continuación se muestran las proyecciones de histograma suavizadas:

Como es normal, si usamos un valor bajo de varianza, obtendremos un resultado más abrupto, lo cual, nos puede llevar a equívocos a la hora de determinar el número y la posición de los ojos. Por el contrario, si el valor de varianza es excesivamente alto (o lo fuera el del número de muestras de la función gaussiana con la que convolucionamos), podríamos perder información acerca de los picos y, por tanto, del número y posición de los ojos.

Mediante un algoritmo de detección de pendiente, conseguimos contar el número de picos sobre las proyecciones de histograma (lo cual nos revela información acerca del número de ojos y las coordenadas de sus centros). Mediante otro algoritmo de detección de pendiente, conseguimos calcular el número de valles a fin de caracterizar el ancho de los lóbulos (esto nos da información acerca del radio del ojo, tanto horizontal, como vertical). Para el ejemplo que nos ocupa obtenemos los siguientes datos mediante el algoritmo que implementé en MatLab:

                                   Posición de picos Posición de valles

Proyección Horizontal     83, 103                71, 93 y 115

Proyección Vertical         114, 242              102, 127, 230 y 255

Pero, ahora surge una duda. Normalmente, el número de picos de la proyección horizontal es igual al de la proyección vertical y, estos son iguales al número de ojos rojos en la fotografía. Sin embargo, tenemos N^2 posibles ojos rojos debido a las combinaciones posibles entre picos de ambas proyecciones (en nuestro caso particular, hasta 4, donde realmente hay 2). Para solventar esta ambigüedad, sumamos todos los niveles de gris de cada uno de los rectángulos cuyos lados quedan definidos por la posición de los valles y cuyos centros por las distintas combinaciones de los picos en la última imagen con la región de los ojos en blanco. Si la combinación de picos no se corresponde con un ojo, la suma de este rectángulo será eminentemente 0. Sin embargo, si se corresponde, la suma será un valor elevado. El punto de corte para considerar si una combinación de picos determina o no el centro geométrico de un ojo, viene definido porque, al menos, el 25% de píxels del rectángulo que comprende al hipotético ojo, sea blanco. Si es así, en la imagen RGB se pinta una cruz con centro el centro geométrico del ojo y cuya longitud de brazos viene determinada por el ancho de los lóbulos (es decir, por la posición de los valles).

Sabido esto, para la anterior imagen obtenemos:

Una vez detectada la posición de los ojos rojos, únicamente habremos de corregir su tonalidad. Aquí entra en juego el último de los parámetros. Este parámetro se basa también en establecer un umbral para la imagen Cb + Cr’, de forma que se analiza cada uno de los rectángulos que comprende aun ojo sobre dicha imagen según los valores antes obtenidos de picos y valles. Si el nivel de gris i-ésimo en el rectángulo que comprende a cada uno de los ojos es menor que el umbral, se considera ojo rojo y se castiga a los niveles RGB con factores 0.10, 0.12 y 0.12, respectivamente (es decir, se tiende a llevar al negro, que sería lo correcto para el color de la pupila). Si es superior el nivel de gris i-ésimo, se entiende que no forma parte del ojo rojo y, por tanto, se deja intacto.

Este parámetro es muy importante y participa, junto con el primero, de forma directa en la calidad final de la imagen retocada. Si el umbral es muy bajo, corremos el riesgo de no corregir el color rojo de los ojos. Sin embargo, si es muy alto, tenderemos a llevar al negro todo el rectángulo que abarca al ojo.

El resultado final para un umbral de corrección de 200, resulta ser:

La forma de utilizar la función implementada en MatLab es la siguiente y se puede descargar haciendo clic aquí:

y = correctorOjosRojos(imagen,nUmb,var,cte)

Donde en y se obtiene la matriz imagen corregida, siendo imagen la imagen de entrada con el defecto de los ojos rojos, nUmb el umbral para binarización aplicado tras el primer stretching, var la varianza de la función gaussiana de suavizado y cte el punto de corte para consideración de píxel rojo o no durante la corrección final.

Reconocimiento automático de partituras

Mi amigo Jonathan y yo hemos programado un pequeño código que intenta ser un reconocedor automático de partituras basado en el procesamiento de la imagen (algo análogo al OCR para el reconocimiento de caracteres).

La idea es que una partitura escaneada pueda ser digitalizada para ser interpretada y procesada su información por un computador (por ejemplo, para conseguir una reproducción automática de la misma).

Sólo reconoce partituras sencillas, funcionando del siguiente modo (caracterizamos cada uno de los bloques funcionales):

Preprocesamiento (OMR.m)

En este primer apartado adaptamos la partitura para su tratamiento. Una vez obtenemos la imagen de la misma, supuesta su entrada con el modelo de color RGB, la transformamos a niveles de gris para seguidamente invertirla. La misión de la inversión es la de un tratamiento más sencillo mediante la contemplación de la información contenida en la partitura como niveles de gris relativamente grandes frente al fondo tendente a nivel cero. Finalmente aplicamos binarización para que los píxels de interés se mantengan en blanco frente al fondo negro.

Corrección de la inclinación (corregirInclinacion.m)

Esta etapa está incluida en el preprocesado de la imagen y como su nombre indica consiste en corregir la inclinación que pueda presentar la imagen de partida. Este paso es de vital importancia dado que los métodos de segmentación y reconocimiento desarrollados no operarían correctamente ante la presencia de este defecto.

El método seguido para llevar a cabo la corrección de la inclinación se basa principalmente en la transformada de Hough. En primer lugar se aplica la transformada de Hough con una resolución en la componente theta de 0.05º. El espacio de parámetros de salida usado es el polar. La ventaja que aporta (aparte de los problemas que nos evitamos con las líneas verticales) es que obtenemos directamente una aproximación del ángulo de inclinación de la imagen. El porqué de optar por el uso de la transformada de Hough es por las características inherentes de la imagen que vamos a emplear. Dado que estamos ante la imagen de una partitura musical que presenta numerosos pentagramas, los máximos de la transformada de Hough se corresponderán con las líneas que componen dichos pentagramas. En consecuencia uno de los máximos de la transformada de Hough nos proporcionará la ecuación de la recta en coordenadas polares asociada a una de las líneas de los pentagramas y a partir de la misma podremos obtener directamente una estimación del ángulo de inclinación. Posteriormente, por si acaso la imagen presenta aún cierta inclinación, se opta por un ajuste de tanteo. Este consiste en aplicar inclinaciones a la imagen en los márgenes de error que tenemos tras aplicar el paso de la transformada de Hough. La realimentación para determinar si las correcciones que estamos realizando están surtiendo efecto o no (y en definitiva para obtener convergencia hacia la solución) se basan en el estudio de la proyección horizontal del histograma de la imagen resultado.

Segmentación de pentagramas (OMR.m)

A continuación segmentamos la imagen separando pentagramas y título de la partitura. Esto lo logramos con ayuda de la proyección horizontal de histograma. Dado que el fondo es completamente negro, las zonas inter-pentagrama se caracterizarán en la proyección horizontal de histograma por una secuencia continua de ceros por lo que, mediante un algoritmo de búsqueda de ceros basado en comparación con el valor anterior y posterior, establecemos el punto medio entre dos pentagramas para que estos sean separados.

Segmentación de símbolos (segSimb.m)

Ahora, para cada uno de los pentagramas, así como para el título, la función programada en MatLab se encarga de segmentar cada uno de los símbolos que componen la línea con la ayuda de la proyección vertical de histograma de cada una de las líneas por separado. En el caso de un pentagrama, en su proyección vertical de histograma se podrá observar un nivel mínimo de continua correspondiente a la contribución de las cinco líneas del pentagrama. Este nivel, que conocemos como el mínimo no-cero, es detectado para cada uno de los pentagramas individualmente de modo que, con un algoritmo similar al anteriormente presentado para la segmentación de la proyección horizontal de histograma, segmentamos los símbolos de cada una de las líneas detectadas (con la diferencia del uso de un algoritmo de búsqueda del mínimo no-cero comentado). Estos límites intersimbólicos podrían aprovecharse para la segmentación directa del elemento en la imagen. No obstante, se prefirió realizar una caracterización del símbolo mediante el estudio del contorno de la proyección vertical de histograma del mismo, ya que podemos limitar su extensión mucho más fácilmente respecto de la extensión del símbolo en el interior del pentagrama (comentar en este sentido que la caracterización del símbolo mediante el estudio de seis momentos estadísticos arrojó peores resultados respecto del método finalmente implementado).

Reconocimiento de símbolos (segSimb.m y funClasif.m)

A continuación se procede a la clasificación de cada uno de los símbolos detectados según el siguiente patrón. Ya hemos dicho que realmente estamos segmentando la proyección de histograma vertical del símbolo. Pues bien, con esta información hacemos lo siguiente:

- Eliminación de la componente de continua detectada como el mínimo no-cero anteriormente comentado.

- Normalización de los valores de la proyección vertical de histograma del símbolo.

- Obtención de 24 coeficientes del módulo de la FFT de la proyección vertical de histograma del símbolo.

- Comparación con cada uno de los patrones base almacenados por clase en la base de datos. La selección de la clase a la cual pertenece el símbolo se realiza en base a la minimización del error lineal calculado como la acumulación de las diferencias en valor absoluto entre cada par de coeficientes de los módulos de las FFT’s del símbolo en estudio y del patrón de la base de datos con el cual se compara.

El estudio del símbolo en base a la FFT se basa, junto con la normalización, en la importancia de caracterizar con un número standard de coeficientes el símbolo, de modo que se contemple la propiedad de la invarianza a la escala.

La función de MatLab asume que el símbolo que se estudia se va a encontrar en la base de datos, por lo que se selecciona mediante minimización del error y no se establece ningún umbral necesario tipo tolerancia para poder clasificarse o no.

Otra limitación importante es que la función también asume que el título sólo ocupa la primera línea segmentada mediante la proyección horizontal de histograma. Esto se ha pensado así debido a la forma de segmentar un símbolo dentro de un pentagrama. Recordemos que el símbolo se lleva sobre el nivel cero (eliminación de la componente de continua introducida por las cinco líneas de pentagrama). Pues bien, el programa hace una excepción sobre el cálculo del nivel mínimo de continua para la primera línea ya que, si se selecciona un mínimo no-cero en la línea del título, podemos perder y distorsionar la información no proveyendo un resultado coherente.

Comentar que sobre los cinco sencillos ejemplos provistos se obtienen unos resultados excelentes con este método de clasificación.

Los símbolos contemplados en la base de datos son los siguientes. Notar que, para mayor precisión en la selección, por cada símbolo se incluyen varias caracterizaciones en función de si este puede aparecer reflejado, invertido (piénsese en una figura como una blanca)…

El reconocimiento de más símbolos aparte de los incluidos en la base de datos sería inmediatamente extensible, pudiendo cubrir con este método partituras de mayor complejidad si fuese necesario.

Como se observa, tampoco se ha incluido el reconocimiento de armaduras por simplicidad (sólo se admiten partituras cuya tonalidad sea de Do mayor o su relativo menor La menor) ni alteraciones accidentales. Por tanto, las partituras que se proporcionan están todas transportadas a la tonalidad de Do mayor.

Detección de tono (deteccionTono.m)

Una vez tenemos segmentada la imagen y se han reconocido las figuras musicales se procede a la detección de tono donde se requiera. El método seguido para la detección de tono está basado en el uso de la proyección horizontal del histograma. Dado el trozo de la imagen hallamos su proyección horizontal (sobre el eje y). A continuación se sigue un método en el que se determina tanto la posición de las líneas del pentagrama como las líneas virtuales (se pueden definir 13 líneas desde el Re 7 al Sol 3, 5 visibles y 8 virtuales). Lo que hacemos es recorrer la proyección horizontal hasta determinar las líneas del pentagrama. Como las líneas del pentagrama estarán representadas en blanco y ocuparán todo el trozo de la imagen se corresponderán con máximos. Así determinamos sus posiciones y las distancias de separación entre cada par de líneas. Haciendo un promedio de la separación entre pares de líneas adyacentes visibles podemos obtener una estimación de la separación (aplicando redondeo por estar en un espacio discreto) entre líneas. Con esta separación estimada podemos determinar las posiciones de las 8 líneas virtuales restantes.

Una vez determinadas todas las posiciones vamos recorriendo la proyección del histograma situándonos en las líneas estimadas y en el punto medio de ambas. En cada posición haremos la suma de todos los elementos del histograma en un entorno igual a la mitad de la separación estimada.

Finalmente en aquella posición (de las 26 posibles, líneas y puntos medios entre dos líneas) donde se obtenga una suma mayor se considera que se encuentra la figura musical y en consecuencia podremos establecer.

Los distintos ficheros y partituras de prueba se pueden descargar aquí. El fichero demo.m ejecuta un ejemplo, obteniéndose el resultado en un fichero de texto con los distintos símbolos reconocidos y las frecuencias auditivas de las notas musicales.

Limpieza de imágenes biomédicas

La limpieza de imágenes biomédicas consiste en obtener una imagen de la muestra de interés lo más pura posible. Esta imagen, normalmente se obtiene aplicando operaciones matemáticas sencillas, basadas en transformaciones globales, sobre las imágenes obtenidas a través del microscopio con el portaobjetos con y sin muestra.

La función que finalmente he desarrollado es muy sencilla y en ella no se puede controlar ningún parámetro, sino que devuelve un resultado unívoco para las imágenes de entrada. Si la función del parámetro está bien definida, es interesante no pasar por alto su contribución en un problema de índole aleatoria como es este, donde no conocemos las características de las imágenes de entrada (no con total determinismo).

En primer lugar restamos píxel a píxel el fondo a la imagen de la muestra contaminada con el ruido. Esto lo hacemos así y no a la inversa debido a que la muestra contrasta con el fondo por estar a contraluz, es decir, la muestra es más oscura. Si lo hiciésemos al contrario perderíamos mucha información siendo, en algunos casos, imposible dar una solución medianamente aceptable al problema.

Una vez restadas las imágenes, obtendremos un resultado bastante limpio, como es lógico, aunque con poco contraste. Por ello, le aplicamos stretching, cuya función de transferencia adjunto en la siguiente figura, así como el código de la función implementada en MatLab.

H(l) = 255(l - mín)/(max - mín)

Donde max y mín son los niveles de gris máximo y mínimo presentes en la imagen.

% STRETCHING
% Autor: Iván López Espejo
%
% y = stretch(x)
%
% Donde en 'y' se obtiene la imagen tras aplicar la operación, siendo % 'x' la imagen de entrada.
function y = stretch(x)
% Algunos parámetros necesarios.
minimo = double(min(min(x)));
maximo = double(max(max(x)));
filas = length(x(:,1));
columnas = length(x(1,:));
% Definición de la función de transferencia y aplicación.
for j = 1:1:filas
for k = 1:1:columnas
y(j,k) = round((255/(maximo-minimo))*double(x(j,k))-minimo*(255/(maximo-minimo)));
end
end

Finalmente sólo nos queda invertir los niveles de gris de la imagen resultante a fin de obtener el objeto de interés en negro sobre el fondo blanco. El siguiente código implementa todo el sistema haciendo uso de la función anteriormente descrita. Además, las siguientes figuras recogen algunos resultados obtenidos sobre ejemplos reales.

% FUNCIÓN PARA LA LIMPIEZA DE IMÁGENES RUIDOSAS
% Autor: Iván López Espejo
function y = Biolimp2(img,fon)

% Diferencia del fondo y la imagen.
img2 = fon - img;
% Stretching del resultado.
y2 = stretch(img2);
% Inversión del resultado y muestra por pantalla.
for j = 1:1:length(y2(:,1))
for k = 1:1:length(y2(1,:))
y(j,k) = 255 - y2(j,k);
end
end
subplot(2,1,1), imshow(img), title('Imagen original')
subplot(2,1,2), imshow(y), title('Imagen limpia')

Pueden existir muchas variantes del anterior método mucho más complejas como, por ejemplo, extraer estadísticamente algún parámetro de la imagen para delimitar el tipo de crecimiento (entre exponencial y logarítmico) de una función de stretching en forma de S, que “mejorase” el contraste entre el fondo y la muestra biomédica. No obstante, dado que el análisis final será personal y cualitativo, estas pequeñas diferencias no aportarán nada o prácticamente nada, por lo que pienso que es preferible una función sencilla y eficiente a algo más complejo que será desaprovechado finalmente por el ojo humano (en principio).