ES2389982T3

ES2389982T3 - Procedimiento y sistema para el análisis de singularidades en señales digitales

Info

Publication number: ES2389982T3
Application number: ES08842547T
Authority: ES
Inventors: Antonio María TURIEL MARTÍNEZ
Original assignee: Consejo Superior de Investigaciones Cientificas CSIC
Current assignee: Consejo Superior de Investigaciones Cientificas CSIC
Priority date: 2007-10-26
Filing date: 2008-10-24
Publication date: 2012-11-05
Anticipated expiration: 2028-10-24
Also published as: US20110075903A1; BRPI0816572A2; WO2009053516A1; EP2214125A4; ES2322120A1; CN101911099A; EP2214125B1; RU2010119650A; JP2011501308A; EP2214125A1; ES2322120B1

Abstract

El método propuesto comprende determinar para cada punto de la señal un entorno que comprende los primeros vecinos y calcular paracada punto de la señal una medida de reconstructibilidad, o medida de singularidad, ponderando las contribuciones de dicho entorno local del punto. El valor de la señal en cada punto se infierea partir de los valores de la señal en los puntos del entorno local utilizando una fórmula de reconstrucción. La medida de singularidad incluye la diferencia entre el valor de la señal en el punto y el valor estimado por su entorno local. Sobre dicha medida desingularidad se realiza una transformación logarítmica con el fin de obtener una medida independiente de la amplitud de la señaldigital muestreada y que varíe de una manera controlada bajo cambios de resolución. Se propone un sistema con medios para obtenerpara cada punto dicha medida de singularidad y para realizar la referida transformación logarítmica y otros cálculos.

Description

Procedimiento y sistema para el análisis de singularidades en señales digitales

5 CAMPO DE LA INVENCIÓN

Esta invención se refiere al análisis de señales digitales de fenómenos físicos, en otras palabras, señales que están muestreadas a intervalos regulares, aplicando para la puesta en práctica del procedimiento propuesto un análisis de ondículas que permite identificar algunos subconjuntos de puntos particulares en el espacio a una escala y posición

10 que son más informativos acerca de la señal que otros.

Asimismo, la invención se refiere a un sistema para la implementación del procedimiento propuesto.

A lo largo de toda esta memoria descriptiva se entenderá como señal digital cualquier colección estructurada de

15 datos muestreada uniformemente que se pueda representar mediante una matriz multidimensional cuyas posiciones se denominan puntos de la señal.

La invención aporta técnicas y herramientas útiles para el procesado, reconstrucción y compresión de señales digitales en base a información parcial acerca de su gradiente y, en particular, operando sobre la base de medidas

20 basadas en gradientes obtenidos por incrementos finitos. Dichas técnicas y herramientas pueden implementarse de manera ventajosa por medio de algoritmos automáticos materializados en programas de ordenador que pueden ejecutarse en entornos computacionales.

La invención, dada la gran eficiencia que proporciona principalmente en tareas de reconstrucción de señales

25 digitales, en particular aquellas que representan imágenes, encuentra aplicación en numerosos campos entre los cuales cabe citar como aplicaciones específicas la compresión de señales digitales (incluyendo la compresión de imagen) y la evaluación de líneas de flujo en señales referidas a fluidos (incluyendo la determinación de líneas de corriente en imágenes de fenómenos físicos); y como aplicaciones más generales la detección de estructuras y el reconocimiento de patrones en imágenes de entornos reales, tales como imágenes fotográficas, geofísicas,

30 biomédicas y de otros tipos.

La invención concierne señales definidas en cualquier número de dimensiones, aunque una vez descrito el procedimiento para un número determinado de dimensiones (por ejemplo dos), resultará bastante evidente para un experto en la técnica generalizarlas para señales definidas en cualquier número de dimensiones. Por este motivo, y

35 en aras de simplicidad, muchas de las ecuaciones y derivadas presentadas a lo largo de esta memoria descriptiva se han escrito para señales 2D, es decir bidimensionales, susceptibles de constituir elementos tales como imágenes. Sin embargo, también se han obtenido resultados útiles en otros números de dimensiones y, en particular, en el procesamiento de señales 1D.

40 ANTECEDENTES DE LA INVENCIÓN

Las patentes US-A-5901249, US-A-6141452 y US-A-6865291 se refieren a técnicas de compresión de señales digitales utilizando análisis de ondículas.

45 La patente US-A-6434261 describe un procedimiento para la detección y segmentación de imágenes digitales con el fin de localizar objetivos en dichas imágenes en base a una determinación de un umbral adaptativo para realizar un análisis de ondículas de las imágenes digitales que se descomponen en diferentes canales de escala.

La patente US-A-7181056 concierne a un procedimiento para la detección automática de regiones de interés en una

50 imagen digital representativa de al menos una porción de un tejido biológico, en la que se genera una representación basada en ondículas de las regiones a explorar.

La patente US-A-7062085 se refiere a un procedimiento para detectar aspectos en regiones de imágenes en color donde se hace referencia a unas características de textura materializadas mediante coeficientes derivados de una

55 transformada de ondícula basada en un análisis de multi-resolución de la imagen digital en color.

La solicitud de patente US-A-2005/0259889 se refiere a un procedimiento para la eliminación de ruido de una imagen de rayos X que comprende la aplicación de una transformación de ondícula compleja a la imagen portadora de un motivo, operando con los coeficientes de ondícula para reducir el ruido.

La solicitud de patente WO-A-2004/068410 concierne a un procedimiento para la detección de puntos de interés en una imagen digital que implementa una transformación de ondícula, asociando una imagen sub-muestreada con una imagen origen.

5 El análisis de singularidad (véase la referencia [14]) que implementa el concepto de caracterizar el comportamiento local de una función f(x) estimada en Rm y definida sobre Rd alrededor de cada uno de sus puntos de dominio x de acuerdo con el denominado exponente de singularidad de Hölder, o exponente Hurst, representado por h(x), es muy útil para muchas tareas de procesado de señales, y en particular, es muy relevante para fines de compresión y como

10 herramienta de reconocimiento de patrones, y, dependiendo del contexto, puede ser utilizado también para revelar información sobre la evolución y dinámica de señales complejas.

La patente US-A-6745129 concierne a un procedimiento basado en ondículas para el análisis de singularidades en datos sísmicos en base al procesamiento de una serie temporal representativa de un registro del fenómeno. El 15 objeto de esta patente es calcular el exponente de Hölder sobre registros sísmicos a través de una transformada de ondícula continua. Usando este procedimiento, al realizar el análisis de la señal (como se muestra en la figura 2b de dicha patente), se producen inestablilidades que repercuten tanto en la resolución espacial, así como en la calidad de la determinación del exponente de Hölder de cada punto (véase el análisis de esta cuestión en la referencia [11]). Esta problemática imposibilita de hecho la utilización del procedimiento del documento US-A-6745129 para tareas 20 de reconstrucción de señales digitales a diferencia de las propuestas del procedimiento de esta invención. Esta invención proporciona una determinación más precisa de los exponentes de singularidad, tanto en cuanto a su posición, como en cuanto a su valor. La diferencia de precisión entre esta invención y el documento US-A-6745129 es debida al uso de medidas de gradiente (que elimina las fluctuaciones indeseables asociadas a las ondículas complejas, véase la referencia [17]) y también por incorporar dicha medida un indicador del grado de

25 reconstructibilidad. De acuerdo con lo anterior, esta invención permite además reconstruir una señal con gran calidad en base a información parcial, al contrario del procedimiento de la patente US-A-6745129 (véase la referencia [11]).

Dentro del campo del análisis de señales basado en ondículas, aplicado en particular al procesamiento de señales

30 digitales, uno de los procedimientos utilizados más conocidos es el denominado Máximos del Módulo de la Transformada de Ondícula (conocido como WTMM) que está determinado por los máximos locales de las proyecciones de ondículas. Mallat y Zhong (véanse las referencias [4], [5] y [6]) conjeturaron que este conjunto puede ser utilizado para reconstruir de manera completa la señal. Con posterioridad se ha verificado que el conjunto conduce a una señal atenuada y que se han de introducir diversos coeficientes empíricos con el fin de poder

35 reproducir las amplitudes correctas de la señal. Desde la publicación del documento de Mallat y Zhong ha habido múltiples intentos de obtener reconstrucción de alta calidad en base a WTMM. En cualquier caso, lo más interesante acerca del procedimiento WTMM es que en el caso de imágenes la mayor cantidad de líneas están concentradas alrededor de los bordes y contornos; y puesto que desde hace años se conoce (véase la referencia [7]) que los bordes y contornos contienen la mayoría de la información de una escena visual, el WTMM se ha evidenciado como

40 un buen candidato para extraer información perceptual (bordes) utilizando un algoritmo automático canónico en base a dicho procedimiento.

Otra rama de investigación también centrada en el uso del WTMM se inició por Arneodo y colaboradores (véanse las referencias [1] y [2]) que reconocieron la capacidad de que el WTMM pueda tratar con señales multiescala,

45 centrando sus estudios en sistemas en los que se conocía que presentaban propiedades de invariancia de escala, tales como flujos turbulentos y sistemas caóticos.

El principal inconveniente de todas las propuestas basadas en el WTMM es la imposibilidad de extraer de manera sistemática los máximos cuando éstos se acumulan (topológicamente), situación que ocurre siempre cuando se trata

50 de señales reales, como se analiza en [17]. En [10] se estudió el problema y se propusieron soluciones parciales.

Al margen del uso del WTMM, se conocen en este campo de la técnica las investigaciones recientes del presente inventor, A. Turiel, relativas al análisis de singularidades en base a medidas de gradiente. En estos trabajos se define una medida de gradiente μ asociada a una señal s(x) sobre un conjunto arbitrario A como la integral del

55 módulo de gradiente sobre ese conjunto:

(1)

Esta investigación de A. Turiel muestra que cuando se trabaja sobre datos reales, discretizados y con ruido, es necesario operar con transformadas de ondículas de las medidas de la manera siguiente: dada una ondícula Ψ, la 5 transformada de ondícula de la medida de gradiente μ a una escala r y en un punto x viene dada por la expresión:

donde d es la dimensión de la señal.

10 La transformada de ondícula de la medida μ permite determinar el exponente de singularidad local, ya que el término dominante cuando la escala r es pequeña depende de r como una ley de potencias (véase la referencia [14])

15 Con la introducción de medidas de gradiente es posible mejorar la resolución espacial de los exponentes de singularidad (véanse las referencias [11] y [14]). De este modo, en lugar de tener un exponente de singularidad cada diez píxeles o de precisar un control de las oscilaciones debidas a la ondícula, es posible asignar un exponente de singularidad a cada píxel con una dispersión mínima del punto. Se aprecia que existen diferencias en las

20 capacidades de resolución de las diferentes ondículas, por lo que se ha reconocido la necesidad de buscar una ondícula optimizada apta para tratar con datos discretizados.

Un elemento importante en la construcción de ondículas con capacidad de resolución optimizada es el concepto de reconstrucción de señales a partir de información parcial acerca de su gradiente. Los planteamientos teóricos e 25 implementaciones prácticas acerca de un algoritmo de reconstrucción de gradiente aparecen introducidos en la referencia [12], que presenta un análisis acerca de la estructura de señales multifractales (véanse las referencias [8],

[3] sobre la estructura multifractal de los fluidos turbulentos). Para señales con estructura multifractal, el conjunto asociado al vértice superior de la jerarquía es bien conocido, al menos desde un punto de vista teórico, y se conoce como la Variedad Más Singular (MSM), que es el conjunto que comprende los puntos con valores más singulares (es

30 decir, más negativos) de h(x).

La referencia [12] que se ha mencionado anteriormente sostiene la teoría de que la MSM contiene suficiente información para reconstruir de manera completa la señal, analizando la reconstrucción de imágenes, aunque las fórmulas son válidas para cualquier número de dimensiones. La fórmula de reconstrucción que se obtuvo en [12]

35 para dominios infinitamente grandes es como se indica a continuación:

(4) en la que

-: s es una señal dada,

es la variedad más singular o MSM de dicha señal s

-

es el vector de gradiente esencial de s, es decir, el gradiente restringido a la MSM, cuyas componentes son

, ; ;

es el nucleo vectorial de reconstrucción universal, que en el espacio de Fourier

viene dado por la expresión:

10 y

de convolución, es decir,

, donde * se refiere al producto de convolución ordinario de funciones.

15 A partir de los estudios realizados por este inventor (véase la referencia [12]) tomando en consideración las medidas de gradientes se ha concluido que, de existir, hay únicamente un algoritmo posible para reconstruir señales usando el gradiente en base a la MSM en dominios infinitamente grandes. Asimismo se sabe que la MSM conduce a muy buenas reconstrucciones con este algoritmo (véanse las referencias [11] y [12] para imágenes, y [15] , [16] para

20 series temporales).

El artículo "Reconstructing Images from their Most Singular Fractal Manifold", [TURIEL, Antonio y col.; IEEE Transactions on image processing, IEEE Service Center, Piscataway, NJ, US, vol. 11, Nº 4, 1 de abril de 2002; XP0111074255] describe un procedimiento para reconstruir una imagen del mundo real usando el gradiente sobre

25 su MSM (variedad más singular).

BREVE DESCRIPCIÓN DE LA INVENCIÓN

Esta invención se define por las reivindicaciones adjuntas y propone un enfoque para el análisis de singularidades 30 en señales digitales de fenómenos físicos, que comprende

a) determinar para cada punto de la señal un entorno de primeros vecinos C; y b) calcular para cada punto x de la señal una medida de reconstructibilidad o capacidad de reconstrucción proporcionada por el entorno C (que se denominará, indistintamente, "medida de singularidad”) en base al entorno

35 que se ha mencionado anteriormente, construida a partir de la inferencia del valor de la señal en x de acuerdo con el valor de la señal en los puntos de dicho entorno de x usando la función de reconstrucción explicada en la referencia

[12] e indicada en la fórmula (4) anterior, pero adaptada al entorno, de manera que se obtiene una medida de singularidad que contiene la diferencia entre el valor medido del punto y el valor inferido a partir de su entorno.

40 El procedimiento propuesto comprende además ventajosamente una tercera etapa c) en la que se realiza al menos una transformación logarítmica sobre dicha medida de reconstructibilidad que suprime la dependencia de la medida del número total de puntos de la señal, obteniendo de este modo un exponente de singularidad para cada punto de la señal.

45 En una realización mejorada, el procedimiento propuesto comprende las siguientes etapas: a1) obtener una función derivada estable de una señal digital, que está muestreada a intervalos regulares; b1) obtener para cada punto de dicha señal digital, una medida de singularidad de la función en ese punto, ponderando las contribuciones de un entorno local del punto (usando una función de reconstrucción como se ha indicado anteriormente) y el valor de la derivada en todos los puntos de dicho entorno local, y

5 c1) realizar al menos una transformación logarítmica sobre dicha medida de singularidad con el fin de obtener una medida independiente de la amplitud de la señal digital muestreada y que varíe de una manera controlada con cambios de resolución.

La invención también comprende la generalización de la medida de singularidad como se ha descrito en el apartado

10 de antecedentes anterior, al proponer una nueva medida de singularidad basada en el conjunto o variedad de puntos impredecibles, UPM por sus siglas en inglés; en otras palabras, se parte de considerar la agrupación de todos los puntos impredecibles, en oposición a los demás puntos que son predecibles.

DESCRIPCIÓN DETALLADA DE LA INVENCIÓN

15 El objeto de esta invención, como se ha indicado en la sección anterior, comprende el cálculo de una medida de reconstructibilidad para calcular de manera precisa los exponentes de singularidad de señales digitales de fenómenos físicos, y que dichos exponentes permitan obtener reconstrucciones de gran calidad. Los requisitos básicos para definir una medida de singularidad μ basada en el conjunto de puntos impredecibles UPM son los que

20 se indican a continuación:

i) La medida μ debe hacer referencia al comportamiento singular local de las funciones. ii) La medida μ debe conducir a una variedad más singular MSM tan cercana al conjunto de puntos impredecibles UPM como sea posible.

25 Las medidas del conjunto de puntos impredecibles son medidas de singularidad que también tienen en cuenta el grado de predecibilidad de los puntos según la ecuación

(5)

30 donde F es la UPM y el superíndice "c" representa el conjunto complementario, en otras palabras, los puntos predecibles; esta ecuación (5) es una consecuencia de la ecuación (4) como se describe en [12]. Por lo tanto, la ecuación (5) muestra que la divergencia del gradiente tomado sólo sobre los puntos predecibles se anula.

El inventor propone en este documento que la mejor manera de continuar trabajando sobre singularidades es definir

35 las medidas basadas en el conjunto de puntos impredecibles como proyecciones de ondículas de medidas de gradiente estándar. De este modo, la medida del conjunto de puntos impredecibles es una proyección de ondícula de la medida de gradiente expresamente diseñada para que se penalice la impredecibilidad. Esto conlleva generalizar el concepto de proyección de ondícula, a fin de producir proyecciones de ondícula con valores vectoriales. El uso de proyecciones de ondícula con valores vectoriales es bien conocido desde hace algún tiempo y

40 no introduce especiales complejidades en la forma de abordar el problema.

Otra diferencia principal en relación con el análisis de singularidad estándar detallado en los antecedentes que se han mencionado anteriormente es que en la propuesta de esta invención no se realizan proyecciones de ondícula de la medida de singularidad a diversas escalas r para extraer los exponentes de singularidad mediante una regresión

45 logarítmica aplicada a la ecuación (3).

Ha de destacarse que proyectar la medida sobre una ondícula a diversas escalas es costoso en tiempo de

50 computación y sólo sirve para mejorar la resolución de las estructuras menos singulares a expensas de empeorar la de las más singulares (véase la argumentación al respecto en la referencia [17]). Dado que el objetivo fundamental en relación a esta invención es extraer las estructuras más singulares, resulta perjudicial realizar las proyecciones a lo largo de múltiples escalas; en lugar de ello, se propone utilizar un estimador puntual (véanse las referencias [17],

[9]) de los exponentes de singularidad, a saber:

Ψμ o Ψμ fo

= +

1

fo

(7)

fo

5 donde {Tl μ(·,r0)} es la media de la proyección de ondícula a lo largo de toda la señal y sirve para disminuir la amplitud relativa de la corrección o(1/log r0). Al aplicar la ecuación (7), se precisa que r0 sea lo suficientemente pequeño para despreciar esta corrección. La escala r0 se define como la menor accesible, es decir, la escala de un píxel. Convencionalmente se asigna una medida de Lebesgue de 1 a la totalidad del dominio espacial, con lo que, para el caso de una imagen de N x M píxeles, el valor de r0 se fijaría en:

fo = (8)

NM

así que, en general, se precisa que las imágenes sean lo suficientemente grandes para hacer del primer término de la parte derecha de la ecuación (7) una buena aproximación al exponente de singularidad. Esto implica típicamente

15 tener imágenes de al menos 100 píxeles en una de las direcciones.

Un aspecto importante de la presente invención reside en el diseño de ondículas digitales con el fin de implementar medidas de singularidad basadas en variedades de puntos impredecibles. Más adelante se presentan dos implementaciones de medidas basadas en variedades de puntos impredecibles de este tipo, que proporcionan un

20 buen resultado en aplicaciones prácticas. El diseño está orientado, en su conjunto, al procesamiento de señales digitales y, en consecuencia, las ondículas se definen (implícitamente) mediante unos pesos numéricos, aunque la presentación está basada en una teoría y es fácil de generalizar a un esquema continuo.

Otro elemento importante del procedimiento propuesto reside en la manera de definir y/o establecer estimaciones

25 numéricas estables del gradiente ∇s para que el núcleo de reconstrucción sea estable numéricamente. Para ello se proponen dos opciones posibles: diferencias de un píxel o punto hacia la derecha y diferencias de medio píxel, en otras palabras, las diferencias de valor al desplazarse una posición a la derecha en el primer caso, o la interpolación equivalente a la diferencia que se obtendría al desplazarse media posición a la derecha y a la izquierda del punto, en el segundo caso. Ambas están definidas por núcleos de derivación descritos en el espacio de Fourier.

30 En otras palabras, la derivada estable de la etapa a1) del procedimiento, que se ha mencionado anteriormente, se obtiene por derivada de incrementos hacia la derecha de un punto o centrados de medio punto.

En las fórmulas que se indican a continuación se caracterizará ∂x, aunque la caracterización de ∂y es análoga. Este

35 operador actúa sobre una señal digital multiplicando simplemente la transformada de Fourier de la señal por el núcleo de derivación, y después anti-transformando el resultado. Se asume que hay Nx píxeles o puntos en la dirección x y Ny en la dirección y.

Diferencia de un píxel/punto hacia la derecha: 40

n

2π.

N

e

−

1

∂

(9)

()=

n

Diferencia de medio píxel/punto:

n

.sen π ; n <N

2

()n = N

∂

(10)

N −n

−.sen π ; n ≥N

N

Otro aspecto básico del procedimiento propuesto consiste en la introducción del nuevo concepto de transformada de Fourier en cruz. Para estimar el grado de predecibilidad de un punto dado, se aplica una fórmula de reconstrucción, 5 que se expresa por la siguiente ecuación

(11)

para el menor número de vecinos posible de un punto, concretamente sus primeros vecinos que constituyen el

10 conjunto C. Debe apreciarse que C no es la MSM F∞; se incluyen puntos en C si hay vecinos del punto x, ignorando el valor de los exponentes de singularidad en ellos. En 2D (siendo d la dimensión de la señal; por lo tanto, d = 2 en este caso) esto consiste en 4 puntos vecinos, que forman, junto con el punto original, una cruz. Para cualquier cantidad p(x) se representan los vecinos de cualquier punto x0 mediante un vector de 5 componentes que comprende dicho punto y sus cuatro vecinos más cercanos, siguiendo la convención de indexación señalada en el

15 siguiente dibujo que ilustra la representación esquemática de la indexación de los puntos en la cruz en 2D. De este modo se asignará al punto central el índice 0, al punto de su derecha el índice 1, al punto de su izquierda el índice 2, al punto de encima suyo el índice 3 y al punto de debajo suyo el índice 4. Por lo tanto, el entorno de primeros vecinos del punto en estudio se convierte en el vector (p0, p1, p2, p3, p4).

En otras palabras, respecto al centro de la cruz, la posición del resto de puntos se corresponde con los desplazamientos de ±1 (en unidades de un punto) ya sea en la dirección x o en la dirección y. Para definir una transformada de Fourier especializada o adaptada a esta configuración en cruz, se ha de tener en cuenta que la

25 frecuencia básica de Nyquist en cada dirección es de 2π/3. Para simplificar la notación se introduce el elemento complejo básico ζ:

2π.

ζ=ζ +ζ=e

= (2π

)+.sen(2π

) (12)

.

30 De acuerdo con la invención, se define la transformada directa de Fourier en cruz de cualquier vector de 5

componentes P =(PPP PP4 ) como el vector complejo de 5 componentes

o 12

P=(P PP PP4 )

o12

35 obtenido a partir de la siguiente fórmula:

P = FP (13)

en la que F es la siguiente matriz compleja de 5 x 5:

1
1
1
1
1

1: ζ ζ *
1
1

5: = 1F 1 ζ * ζ 1 1 (14)

1
1
1: ζ ζ *

1
1
1: ζ * ζ

Esta matriz representa la combinación lineal de los armónicos asociados a los desplazamientos en la cruz y está diseñada para representar con la máxima fidelidad la composición en el centro de la cruz, en base a los puntos más cercanos. La inversa de esta matriz se puede calcular fácilmente,

− 11 1 1 1 1 ζ * ζ oo

F−1 = 1 ζζ *o o (15)

1oo ζ * ζ 1oo ζζ *

y esta última matriz es necesaria para realizar la anti-transformada de Fourier en cruz.

15 Es necesario definir implementaciones del gradiente y la fórmula de reconstrucción restringida a la cruz, a fin de evaluar rápidamente el grado de predecibilidad del punto central en función de los vecinos. Por este motivo, esta invención propone unas implementaciones adecuadas del gradiente y de la fórmula de reconstrucción de dicho gradiente, sobre la base de la transformada de Fourier en cruz.

20 Una primera implementación es la del operador de gradiente en cruz en funciones de operador de gradiente local, que es el operador

(∂ ∂)= F−1 ⋅(∂ ∂)⋅ F .

En el espacio de Fourier, dicho operador actúa simplemente multiplicando cualquier función por las funciones 25

∂

y

30 ∂

para obtener las coordenadas x e y, respectivamente. La función

∂

35 se define para entornos en cruz como se indica a continuación:

(16)

y, de manera análoga, se tiene:

(17)

que se definen de tal manera que representan diferencias de medio píxel; de hecho

) .

= 2sen(π

Una segunda implementación es la del operador de reconstrucción en cruz, que es una de las inversas del operador

10 de gradiente en cruz. Como el operador de gradiente elimina cualquier constante sumada a cada componente del vector de 5 componentes que representa los vecinos, la reconstrucción está completamente definida excepto por un cambio de esa constante; nuestra implementación del operador de reconstrucción en cruz que se propone es tal que

el vector de 5 componentes resultante tiene una media de cero, � P = o. Por este motivo, las señales deben

..=o

tener la media sustraída antes de aplicarse estos dos operadores. Para ello se aplican los elementos de matriz de la

15 primera línea de la transformada de Fourier en cruz inversa para no introducir armónicos cuando se aplique el operador de reconstrucción. Como la suma de los elementos de esa primera línea es (2 x d) -1 (cuyo resultado es 3 en el caso de señales 2D), para sustraer la media se suman todos los valores del vector entorno (primeros vecinos) y se divide por ((2 x d) -1), y al resultado se le suma la primera componente del vector de entorno y se le restan las otras componentes.

20 = F −1 ⋅⋅ F

La reconstrucción en cruz es el operador . En el espacio de Fourier tiene dos componentes funcionales,

=( );

25 el operador actúa como la suma del producto de cada componente con la componente correspondiente (x e y) del

gradiente sobre el que opera. La componente se define para un entorno de cruz como se indica a continuación:

(18)

30 y, de manera análoga para ,

(19)

De este modo la medida de singularidad definida en la etapa b1) del procedimiento de esta invención puede 35 describirse mediante los siguientes pasos para una señal genérica definida en un espacio de dimensión d arbitraria:

-: se extrae el entorno de los (2 x d) primeros vecinos de un punto de base x, obteniendo los primeros vecinos del punto x, modificando consecutivamente cada uno y sólo uno de los índices de coordenadas de dicho punto x, primero sumando -1 y luego sumando +1, formando un vector de (2 x d)+1 componentes, cuya primera componente

40 es el valor de la señal en el punto x, la segunda el valor de la señal en el punto obtenido al modificar la primera coordenada sumándole -1, la tercera el valor de la señal en el punto obtenido al modificar la primera coordenada sumándole +1, la cuarta el valor de la señal en el punto obtenido al modificar la segunda coordenada sumándole -1 y así sucesivamente;

-: se extrae la tendencia de este vector, que se define como la suma de sus valores dividido por ((2 x d)-1) y esta 45 tendencia se aplica al vector, añadiéndola a la componente referida al punto de base x y sustrayéndola de las otras componentes, de manera que de este modo el nuevo vector obtenido tiene media nula;

-: se aplica un operador de gradiente local al vector de media nula que se ha mencionado anteriormente, lo que devuelve (2 x d)+1 vectores de gradiente, cada uno de ellos de d componentes, que definen el gradiente local:

-: se anulan las componentes de dicho gradiente local asociadas al punto x;

5 - se aplica al gradiente local, con las componentes anuladas, un operador de reconstrucción asociado unívocamente al citado operador de gradiente local obteniendo un vector de (2 x d)+1 componentes, que se denomina señal estimada;

-: se aplica una vez más el operador de gradiente local sobre dicho vector de (2 x d)+1 componentes o señal estimada y se obtienen (2 x d)+1 vectores de d componentes, que definen el gradiente local estimado;

10 - se obtienen (2 x d)+1 vectores de d componentes que expresan la diferencia de gradientes entre dicho gradiente local y dicho gradiente local estimado, y a partir de estos (2 x d) +1 vectores de diferencia de gradientes se obtiene la medida de la singularidad asociada al punto x.

Los operadores de gradiente en cruz y de reconstrucción en cruz son dos de los procedimientos incluidos en esta

15 invención susceptibles de implementarse mediante algoritmos básicos materializados de manera ventajosa en programas de ordenador ejecutables en un entorno computacional, para el diseño y cálculo de las medidas de singularidad basadas en variedades de puntos impredecibles. En particular, tales programas o partes de los mismos pueden incluirse en rutinas almacenadas en microprocesadores o microchips. Estos operadores pueden simplificarse a una forma matricial de ((2 x d)+1) x ((2 x d)+1), para una implementación numérica más rápida.

20: A continuación,invención: se describen dos medidas de singularidad diseñadas de acuerdo con los principios de esta

25: - medida de singularidad de correlación local (lcsm); y - medida de singularidad de correlación global (gcsm).

Ambas medidas pueden implementarse mediante algoritmos específicos materializados de manera ventajosa en programas de ordenador ejecutables en un entorno computacional. En particular tales programas o partes de los 30 mismos pueden incluirse en rutinas almacenadas en microprocesadores o microchips.

La medida de singularidad de correlación local se ha concebido para medir la impredecibilidad de un punto dado, simplemente calculando la diferencia entre el valor real de la señal sin media (es decir, una vez suprimida la media) en un punto dado y el valor inferido a partir de sus cuatro vecinos (cuando d = 2). Esta medida tiene como objetivo

35 evaluar μ( ) en un punto dado x0 y en el caso d = 2 comprende los siguientes pasos:

Ψoo

lcsm

1. Se convierten los vecinos de x0 en un vector de 5 componentes s =(sss ss4) según el esquema de

o 12

indexación en cruz del gráfico que se ha mostrado anteriormente.

2. Se rectifica el vector convenientemente: en primer lugar se obtiene S = s , y el vector rectificado,

..=1

40 P =( P PP PP4), se define como:

o 12

Po = so + S

P. = so − S. = 1...4 (20)

45 3. Se aplica el operador de gradiente en cruz a P para obtener los vectores g y g .

4.: Se conserva el valor de las componentes asociadas al índice 0 de dichos vectores para su uso posterior, A = g o A = g o.

5.: Se ajustan dichas dos componentes a cero, g = g = o.

oo

6.: Se aplica el operador de reconstrucción en cruz a los vectores resultantes g y g , para obtener la señal

reconstruida .

7.: Se aplica de nuevo el operador de gradiente en cruz a para obtener los gradientes estimados ρ y ρ .

8.: Se define la medida de singularidad de correlación local como el módulo de la diferencia de los gradientes en cruz

en el centro de dicha cruz, a saber: 5

μ( )= (A −ρ )2 +(A −ρ )2 (21)

Ψ oo oo

lcsm

De hecho, este último paso significa conservar el módulo de una proyección de ondícula con valores vectoriales, pero para simplificar la notación se deja tal cual está. 10

9. A continuación, se obtiene el exponente de singularidad h(x0) aplicando la ecuación (7).

En otras palabras, la medida de la singularidad asociada al punto x comprende:

15 - retener de los (2 x d) + 1 vectores de diferencia de gradiente local obtenidos las d componentes asociadas al punto x, y - obtener la medida de singularidad como la raíz cuadrada de la suma de los cuadrados de estas d componentes

con lo que se obtiene una medida de singularidad de correlación local apta para medir la impredecibilidad de un 20 punto dado.

La medida de singularidad de correlación global mejora la de correlación local teniendo en cuenta no sólo el tamaño de las desviaciones entre las señales estimadas y las reales, sino también la diferencia entre las direcciones de los gradientes obtenidos. Por este motivo, los datos iniciales no son solamente la señal s , sino también el gradiente

25 ∇s . Es muy importante proporcionar una caracterización estable de ∇s ; para ello se han utilizado los dos núcleos que se han mencionado anteriormente: un núcleo de diferencias de un píxel hacia delante y un núcleo de incrementos de medio píxel.

La medida de singularidad de correlación global tiene una estructura más compleja; sin embargo, el inventor ha

30 comprobado que es la más eficaz para evaluar singularidades y a la vez garantizar una gran calidad de reconstrucción. La obtención de esta medida se realiza en dos etapas: en primer lugar, se obtiene una diferencia de gradiente para todos los puntos; a continuación, se construye la medida en cada punto x0 combinando las

diferencias de gradiente asociadas a dicho punto y el gradiente ∇s en cada grupo de vecinos de dicho punto.

Esta medida tiene como objetivo evaluar gcsm μ( ) en un punto dado x0 y comprende los siguientes pasos:

Ψoo

35 Primera etapa: obtener las diferencias de gradiente en cada punto x0.

1. Se convierten los vecinos de x0 en un vector de 5- componentes s =(s ss ss4 ) según el esquema de

o 12

indexación en cruz del gráfico que se ha mostrado anteriormente.

40 2. Se rectifica el vector convenientemente: en primer lugar se obtiene S = s , y el vector rectificado,

..=1

P =( PPPPP4 ), se define como:

o 12

Po = so + S

P. = so − S. = 1... 4 (22) 45

3.: Se aplica el operador de gradiente en cruz a P para obtener los vectores g y g .

4.: Se conserva el valor de las componentes asociadas al índice 0 de dichos vectores para su uso posterior,

A =g o A =g o .

5. Se ajustan dichas dos componentes a cero, g =g =o.

oo

6. Se aplica el operador de reconstrucción en cruz a los vectores resultantes g y g , para obtener la señal

reconstruida .

5 7. Se aplica de nuevo el operador de gradiente en cruz a para obtener los vectores de gradiente estimado ρy

ρ .

8. Se genera la diferencia de gradiente asociada al punto central, (ε ε)=(ρ−A ρ−A ).

Segunda etapa: se evalúa la medida de singularidad de correlación global combinando las diferencias de gradiente y 10 los gradientes del grupo de vecinos de cada punto conforme a los siguientes pasos:

1. Para cada punto x0, se considera la ventana de 3 x 3 centrada a su alrededor. En esta ventana, cada punto tiene unas coordenadas +(dd ), donde dd pueden tomar los valores -1, 0, 1.

o

2. Se calcula la autoproyección de las diferencias de gradientes en esta ventana, S():

o

15 S()=ε() ε(+(dd ))+ε() ε(+(dd )) (23)

ooo oo dd =−1o1 dd =−1o1

Es decir, la medida de la singularidad asociada al punto x comprende:

20 - tomar el hipercubo d-dimensional que rodea a un punto x dado, formado por los puntos obtenidos cuando se suma -1, 0 ó +1 a los índices de coordenadas, lo que proporciona 3d puntos;

-: retener para cada punto de dicho hipercubo el vector d dimensional la diferencia de gradiente local de ese punto, y

-: sumar estos 3d vectores de d dimensiones, y calcular el producto escalar del vector resultante con el vector d

dimensional de diferencia de gradiente asociado al punto x, lo que proporciona un índice de alineación de 25 diferencias de gradiente local.

El índice de alineación de diferencias de gradiente permite deducir la existencia de una coherencia espacial entre los errores cometidos al prescindir del punto central cuando la señal es reconstruida, lo cual permite diferenciar entre ruido (de orientación aleatoria) y señal coherente.

3. Se obtiene la energía del gradiente asociada a esta ventana, E():

o

E()= (∂s(+(dd ))2 +∂s(+(dd ))2 ) (24)

oo o dd =−1o1

35 4. Se obtiene la energía de diferencia de gradiente del punto x0, e():

o

e()=ε()2 +ε()2 (25)

oo o

40 5. Finalmente se define la medida de singularidad de correlación global como:

S()

o

gcsm μ(o )= e o

E() (26)

Ψ o () o

En este caso, la definición es mucho más compleja y la linealidad se ha perdido por completo, incluso si se consideran las proyecciones de ondícula con valores vectoriales.

Se observa que se ha obtenido la energía de gradiente de dicho hipercubo como la suma de los módulos al 5 cuadrado de los gradientes de cada punto del hipercubo.

Por otro lado, se observa que la medida de singularidad global es el producto de la medida de singularidad de correlación local obtenida como se ha explicado anteriormente, por la raíz cuadrada del valor absoluto del índice de alineación de diferencias de gradiente local dividido por la energía de gradiente del hipercubo.

6. A continuación, se obtiene el exponente de singularidad h(x0) aplicando la ecuación (7) .

La invención descrita hasta este punto puede ponerse en la práctica mediante técnicas de computación que se ejecutaran en unidades operativas o de cálculo. La implementación del procedimiento comprende un sistema para 15 análisis de singularidades en señales digitales caracterizado que comprende en una versión básica:

-: medios para obtener para cada punto de la señal un entorno local que comprende los primeros vecinos; y

-: medios para calcular para cada punto x de la señal una medida de reconstructibilidad, o medida de singularidad, en base al correspondiente entorno asociado a cada punto, construida a partir de la inferencia del valor de cada punto

20 en base al valor de los puntos de dicho entorno utilizando para el cálculo la siguiente función o fórmula de reconstrucción

25 en la que

-: s es una señal dada,

-: C es el entorno local de x,

es el gradiente esencial de s sobre ese entorno, 30 - g es el núcleo de reconstrucción universal que se ha explicado anteriormente, y

-: el símbolo · significa producto escalar de convolución,

conteniendo dicha medida de singularidad la diferencia entre el valor medido y el valor inferido para cada punto.

35 De acuerdo con una realización mejorada, el sistema comprenderá además medios para realizar al menos una transformación logarítmica sobre dicha medida de reconstructibilidad destinada a suprimir la dependencia del número de puntos de la señal y proporcionar un exponente de singularidad para cada punto de la señal y en general:

medios para obtener una función derivada estable de una señal digital, muestreada a intervalos regulares;

40 medios para obtener para cada punto de dicha señal digital muestreada una medida de singularidad de la función en ese punto, ponderando las contribuciones de un entorno local del punto y el valor de la derivada en todos los puntos de dicho entorno local; y medios para realizar al menos una transformación logarítmica sobre dicha medida de singularidad con el fin de obtener una medida independiente de la amplitud de la señal digital muestreada y que varíe de una manera

45 controlada bajo cambios de resolución.

Dichos medios comprenderán en general unidades de cálculo o de procesamiento de datos integradas en el sistema

o materializadas en forma de circuito integrado o unidad de procesado dedicada. Las instrucciones para la ejecución

de las etapas del procedimiento estarán grabadas en programas cargables en las unidades operativas o integradas 50 en circuitos electrónicos.

Se describen a continuación unos ejemplos, que han de ser tomados a título no limitativo de aplicación del procedimiento de acuerdo con la invención en diferentes campos. Todas las señales digitales procesadas en los ejemplos siguientes han sido obtenidas de bases de datos públicas, menos las correspondientes a las figuras 5 y 10.

Todas las señales han sido procesadas utilizando un programa escrito por el inventor en lenguaje C y ejecutado en un ordenador personal con sistema operativo Linux. Los exponentes de singularidad obtenidos han sido convertidos en imágenes digitales por el mismo programa.

5 Ejemplo 1: Detección de estructuras y reconocimiento de patrones

Los exponentes de singularidad permiten reconocer estructuras muy sutiles que son difíciles de detectar a simple vista. Esto es así porque los exponentes miden el grado de transición (es decir, su brusquedad) de la señal en cada punto independientemente de su amplitud real. Esto se puede utilizar para detectar pequeñas modificaciones en un

10 medio y para evidenciar la existencia de nuevas estructuras en imágenes. El rango de aplicaciones cubre todo tipo de imágenes, desde imaginería médica hasta teledetección, así como la detección de fotografías manipuladas.

La Figura 1 de los dibujos indica la detección de ondas internas oceánicas a partir de una imagen de satélite MeteoSat (véase la referencia [14]). La imagen de la izquierda muestra una porción de una imagen del canal visible 15 del satélite MeteoSat V adquirida el 27 de diciembre de 2004 sobre la cresta submarina de Mascareña (área al Nordeste de Madagascar); la imagen tiene una resolución de 2,5 kilómetros x 2,5 kilómetros aproximadamente, y comprende 500 x 500 píxeles (lo que corresponde a un área de 1250 km x 1250 km). Las nubes aparecen como áreas blancas, borrosas, mientras que el mar es el fondo oscuro. La imagen de la derecha muestra los exponentes de singularidad asociados, representados con una paleta que asigna a los valores menores los colores más 20 brillantes. La obtención de dichos exponentes de singularidad llevó alrededor de 10 segundos en un ordenador portátil con dos procesadores Centrino a 1,8 Mhz (los tiempos de los otros ejemplos están referidos al mismo ordenador). Además de una estructura más rica asociada a las nubes y a los flujos atmosféricos, se reveló la existencia de frentes concéntricos oceánicos de hasta 500 km de largo, probablemente ondas internas, en el centro de la imagen. Hoy en día se sabe que tener un buen conocimiento de las ondas internas oceánicas es clave para

25 conocer los procesos de disipación de energía y de mezcla (de nutrientes, dispersión de contaminantes, etc.) en los océanos; a pesar de ello, existe muy poca información, y poco sistemática, sobre las zonas del planeta afectadas por estas ondas. Por ejemplo, las citadas en la imagen de la derecha, a pesar de su enorme extensión (diversos frentes de hasta 500 km de longitud separados hasta 300 km) no se habían publicado hasta la fecha.

30 En la Figura 2, la parte superior muestra una imagen de una proliferación de algas en el lago Mendota (Suiza) en falso color, como una combinación de varios canales para aumentar el contraste de dicha proliferación de algas. En la parte inferior de la figura se muestran los exponentes de singularidad, que se obtuvieron tras apenas 3 segundos de cálculo, representados usando una paleta de tonos de gris en la que los valores menores son más brillantes.

35 En la Figura 3 en la parte superior se muestra una vista de la bahía de los Alfacs (NE España, en el Delta del Ebro) registrada por la banda 8 del LandSat, en fecha indefinida. La resolución de esta imagen es de 2,5 metros, y la zona representada cubre 500 x 500 píxeles. La parte inferior de la figura muestra los exponentes de singularidad (tiempo de cálculo: 10 segundos). Algunas embarcaciones apenas evidentes en la imagen de arriba aparecen con contornos bien cerrados en la imagen de abajo; se observan también diversos frentes de olas.

40 La Figura 4 a la izquierda muestra una mamografía de formato digital con una resolución de 1976 x 4312 píxeles, extraída de un archivo público de mamografías en formato digital (USF Digital Mammography Home Page, http://marathon.csee.usf.edu/Mammography/Database.html) al que se accedió el 10 de octubre de 2008. En la parte derecha de la misma figura se muestran los exponentes de singularidad asociados (tiempo de cálculo:

45 aproximadamente 4 minutos). El análisis revela la estructura de los diferentes tejidos que conforman la mama. Este análisis puede permitir una detección más precoz de lesiones. Al mismo tiempo, la capacidad de detección de las líneas de singularidad independientemente del contraste hace posible reducir la exposición radiológica necesaria para la detección de los patrones.

50 En la Figura 5, la parte superior muestra una imagen de 200 x 200 píxeles del núcleo de una célula de cebolla en interfase, obtenida por microscopía óptica (imagen cortesía de Elisenda Gendra y Mónica Pons, Instituto de Biología Molecular de Barcelona, Consejo Superior de Investigaciones Científicas). Esta imagen fue adquirida a partir de un microscopio confocal SP1 de Leica, en modo transmisión (Nomarski), con iluminación por laser de argón a 488 nm de longitud de onda. En la parte inferior de la misma figura se muestran los exponentes de singularidad asociados

55 (tiempo de cálculo: aproximadamente dos segundos). El análisis de singularidades revela la existencia de líneas coherentes dentro del núcleo y en su periferia, posiblemente asociadas a estructuras relacionadas con elementos del núcleo tales como la cromatina y la membrana nuclear, siendo dichas estructuras difíciles o imposibles de resolver o evidenciar por medios ópticos, en particular en ausencia de cualquier tipo de tinción o marcaje. Los exponentes de singularidad parecen evidenciar por ejemplo la estructura de doble membrana del núcleo periféricamente, y también estructuras asociadas con las fibras de cromatina que llenan el núcleo.

Ejemplo 2: Compresión de imágenes y eliminación de ruido

5 Debido a la fórmula de reconstrucción, es posible regenerar una imagen a partir del conjunto de los puntos más singulares con gran calidad. Dicho conjunto suele estar bastante disperso, constituyendo el 20-30% de los puntos totales. Para completar la descripción, se debe registrar y almacenar el gradiente sobre dichos puntos, y se debe codificar de manera compacta. Se ha constatado que los gradientes cambian con suavidad por las líneas de variedad más singular MSM y se estima factible poder codificarlos de manera compacta. Por lo tanto, la

10 reconstrucción de imágenes a partir de la variedad más singular MSM ha evidenciado tener el potencial de proporcionar códigos de compresión de alta calidad para imágenes.

En la Figura 6, la parte superior muestra una imagen de van Hateren identificada en la referencia [18] como imk01020.imc. Esta imagen se ha obtenido con una cámara CCD con distancia focal de 28 mm y está definida por

15 una matriz de 1536 x 1024 píxeles; los datos están codificados como niveles de gris en 12 bits nominales. La obtención de los exponentes de singularidad llevó unos 50 segundos. En la parte central de la figura se muestra el 30% de puntos más singulares. En la parte inferior se ha reconstruido la imagen a partir de los gradientes sobre la MSM mostrada en la parte central, obteniéndose una calidad medida usando la Relación de Pico Señal a Ruido (PSNR) de 37 dB, lo que indica alta calidad.

20 La Figura 7 muestra cómo la reconstrucción a través de la MSM permite reducir el ruido presente en la señal. En la parte de arriba de la figura se muestra la imagen original (imagen de Lena, estándar IEEE en procesamiento de imagen) de una resolución de 200 x 200 píxeles; en la parte inferior, la reconstrucción a partir de la MSM asociada. Los contornos y bordes contenidos en la MSM se conservan en la reconstrucción, pero las transiciones asociadas al

25 ruido, que no forman frentes coherentes, son eliminadas en su mayor parte, lo cual es particularmente perceptible en algunas áreas de la cara, en la imagen reconstruida.

Ejemplo 3: Determinación de líneas de flujo en imágenes geofísicas y de otros tipos

30 Debido a las raíces teóricas en las que se fundamenta la definición de exponentes de singularidad, los últimos son particularmente útiles cuando se utilizan para analizar imágenes de variables escalares en flujos turbulentos. La teoría predice que las singularidades son advectadas (en otras palabras, arrastradas por el fluido), lo que puede usarse para trazar las líneas de corriente. En esencia, se puede seguir el camino de las corrientes analizando simplemente las imágenes asociadas a la temperatura, la concentración de clorofila y otros indicadores análogos.

35 Los resultados de las singularidades derivadas de la temperatura de la superficie del mar estimados por sensores de microondas (MW SST) embarcados en los satélites Modis Acqua y TRMM se comparan con los mapas de altimetría.

Los datos de altimetría son difíciles de producir y tienen una resolución espacial muy mala, que se ha de filtrar mediante un filtro de paso bajo. Además, para producir mapas de altimetría de calidad, se deben combinar diversos 40 altímetros activos, pero desde 2003 sólo siguen en funcionamiento dos satélites, y pronto sólo uno de ellos estará activo, o incluso ninguno. Por el contrario, la MW SST es mucho más económica, se obtiene sinópticamente en grandes zonas y es fácil de procesar. Tal como muestra la comparación, las singularidades delinean bastante bien los patrones de circulación, demostrando que están canalizadas por el flujo. Por lo tanto, determinar corrientes utilizando el análisis de singularidades emerge con fuerza como una alternativa interesante para sistemas

45 oceanográficos operacionales para la gestión del riesgo medioambiental.

La Figura 8 muestra en la parte inferior los exponentes de singularidad obtenidos a partir de una imagen de temperatura de la superficie del mar obtenida mediante microondas (MW SST) - AMSR-E-TMI mostrada en la parte superior, correspondiente al 1 de febrero del 2003 (imagen descargada de Remote Sensing Systems, 50 http://www.ssmi.com/; tiempo de cálculo del análisis de singularidades: aproximadamente 5 segundos). La zona mostrada corresponde a la corriente del Golfo de México. El mapa de temperaturas está dado en una malla en proyección cilíndrica con una resolución angular constante de ¼ de grado. La parte superior de la Figura 9 muestra el campo de corrientes geostrófico obtenido por la interpolación de cuatro satélites altimétricos, para la misma fecha del 1 de febrero de 2003; la parte inferior de la figura muestra la superposición de los dos campos (exponentes de

55 singularidad de la temperatura de la figura anterior y campo de velocidades geostrófico del panel superior de esta figura).

Ejemplo 4: Análisis dinámico de secuencias de variables en fluidos turbulentos.

Hoy en día la simulación numérica de fluidos es una herramienta indispensable en tareas como la predicción meteorológica y oceanográfica, el prototipado aerodinámico de modelos o diversos problemas de análisis de reacciones químicas y combustión, de interés industrial. Sin embargo, dado el carácter caótico de los fluidos en régimen de turbulencia no es posible hacer una descripción de precisión acotada con un número finito de grados de

5 libertad como los que imponen las mallas discretizadas que se usan en simulación numérica. Este problema es consecuencia de que cuando se describe un fluido con un paso de malla de tamaño determinado no se pueden resolver los movimientos que tienen lugar a escalas más pequeñas, que por el carácter caótico del fluido no se pueden predecir.

10 La estrategia habitual para tratar estas escalas no resueltas es la introducción de unos coeficientes de viscosidad empírica (para el campo de velocidad) y difusividad empírica (para cada variable considerada en la simulación), también llamados coeficientes de viscosidad turbulenta y difusividad turbulenta. Estos coeficientes representan la dispersión más o menos aleatoria y homogénea de la variable considerada en estas escalas sin resolver. Dichos coeficientes sirven para modelizar el efecto de las escalas no resueltas sobre las escalas resueltas por la simulación

15 bajo ciertas circunstancias; por ejemplo, si la turbulencia está completamente desarrollada o si el tiempo de integración de la simulación es suficientemente grande comparado con el tiempo típico de dispersión de las escalas no resueltas.

La determinación de los coeficientes de viscosidad y difusividad turbulentas en fluidos es fundamental para describir

20 la evolución de los fluidos turbulentos mediante modelos numéricos con suficiente calidad. Una buena determinación de estos coeficientes es crucial para tener mayor precisión sobre las variables predichas con la simulación numérica, así como para extender el horizonte temporal de validez de estas predicciones. Sin embargo, en la mayoría de los modelos numéricos que se utilizan hoy en día se toman estos coeficientes como una constante en todo el dominio del fluido. Esta constante se estima de manera heurística para cada realización numérica, aunque se suele

25 comparar su valor con un valor experimental evaluado del análisis dinámico de secuencias, de acuerdo a la fórmula siguiente:

30 donde es el coeficiente de difusividad empírico global, es la variable analizada, el subíndice 0 hace mención a la escala a la cual se está procesando y los paréntesis triangulares significan promedio sobre todo el

dominio espacial del fluido. Si en lugar de

se toma la función de corriente lo que se evalúa es la viscosidad, en vez de la difusividad.

35 En realidad, se puede pasar de la estimación global de la difusividad que se ha discutido anteriormente a una estimación local de este coeficiente para cada punto del dominio. Para ello, los promedios de derivadas temporales y gradientes de la expresión anterior se sustituyen por promedios ponderados con un peso decreciente con la distancia al punto de evaluación. Sin embargo, si la fórmula de evaluación global que se ha presentado anteriormente ya tiene cierta inestabilidad, la fórmula local es extremadamente inestable, dando lugar a valores de la

40 difusividad local negativos en ciertos puntos, lo que es físicamente inaceptable. La aplicación del procedimiento de la presente invención permite obtener una variable más estable (la densidad de la MSM), a partir de la cual se obtienen evaluaciones de la difusividad global muy estables y evaluaciones de la difusividad local que no son negativas en ningún punto.

45 Como muestra de ejemplo, se han evaluado las difusividades locales a partir de imágenes de un colorante dispersado en un experimento de laboratorio. La fila superior de la figura 10 muestra un colorante dispersado en un medio en turbulencia bidimensional en un laboratorio en dos instantes de tiempo, t = 0 s para la columna izquierda y t = 10 s para la columna derecha (Imágenes de la fila superior, cortesía de Patrick Tabeling, École Normale Supérieure, París). En la fila central de la misma figura se muestra la evaluación de la difusividad turbulenta local

50 para todos los puntos en los mismos instantes de las imágenes de la fila superior, usando como variable la concentración de colorante, la cual se estima a partir del nivel de gris de las figuras de dicha fila superior. Los valores de difusividad local obtenidos de esta manera se han representado usando una paleta de dos colores extremos (rojo para los negativos, azul para los negativos), con un color intermedio (blanco) para los valores próximos a cero. Para facilitar la comprensión de la figura se ha superpuesto una trama de trazo grueso y oblicuo a

las zonas con valores negativos de mayor magnitud, y una trama horizontal y más fina a las zonas de mayores valores positivos. Como muestran las imágenes de esta fila central, la estimación de la difusividad en base a la concentración da lugar a amplias zonas con valores negativos; además, cuando se procesa esta secuencia se puede observar que la determinación es muy inestable, ya que los valores estimados de difusividad local de algunas 5 zonas cambian de manera brusca en determinados instantes de tiempo. Por último, la fila inferior presenta las evaluaciones de difusividad local para los dos mismos instantes de tiempo obtenidas en base a la función de densidad de la MSM, la cual se calcula a partir de los exponentes de singularidad evaluados usando esta invención. Como muestra la figura, en estas evaluaciones de la difusividad local no se presentan regiones con valores negativos; además, la observación de la secuencia entera muestra una evolución suave y continua de los valores de

10 difusividad local en todos los puntos.

A continuación, se incluye una serie de referencias a publicaciones científicas del estado de la técnica que reflejan aspectos explicados en esta invención.

15 Referencias

[1]: A. Arneodo. Wavelet analysis of fractals: from the mathematical concepts to experimental reality. In G. Erlebacher,

M.: Yousuff Hussaini y L.M. Jameson, editors, Wavelets. Theory and applications, página 349. Oxford University Press. ICASE/LaRC Series in Computational Science and Engineering, Oxford, 1996.

20 [2] A. Arneodo, F. Argoul, E. Bacry, J. Elezgaray y J. F. Muzy. Ondelettes, multifractales et turbulence. Diderot Editeur, Paris, Francia, 1995.

[3] U. Frisch. Turbulence. Cambridge Univ. Press, Cambridge MA, 1995.

[4] S. Mallat. A theory for multiresolution signal decomposition: the wavelet representation. IEEE Transaction on Pattern Analysis and Machine Intelligence, 11: 67-93, 1989.

25 [5] S. Mallat y W. L. Huang. Singularity detection and processing with wavelets. IEEE Trans. in Inf. Th., 38: 617-643, 1992.

[6] S. Mallat y S. Zhong. Wavelet transform maxima and multiscale edges. In Ruskai M. B. y col., editor, Wavelets and their applications. Jones and Bartlett, Boston, 1991.

[7] D. Marr. Vision. Freeman and Co. Nueva York, 1982.

30 [8] G. Parisi y U. Frisch. On the singularity structure of fully developed turbulence. In M. Ghil, R. Benzi y G. Parisi, editors, Turbulence and Predictability in Geophysical Fluid Dynamics. Proc. Intl. School of Physics E. Fermi, pages 84-87, Ámsterdam, 1985. Holanda del norte.

[9] O. Pont, A. Turiel y C. Perez- Vicente. Application of the microcanonical multifractal fonnalism to monofractal systems. Physical Review E, 74: 061110, 2006.

35 [10] Z. R. Struzik. Determining local singularity strengths and their spectra with the wavelet transfonn. Fractals, 8(2): 163-179, junio de 2000.

[11] A. Turiel. Relevance of multifractal textures in static images. Electronic Letters on Computer Vision and Image Analysis, 1(1): 35-49, 2003.

[12] A. Turiel y A. del Pozo. Reconstructing images from their most singular fractal manifold. IEEE Trans. Im. Proc., 40 11: 345-350, 2002.

[13] A. Turiel, J. Isern-Fontanet, E. Garcia-Ladona y J. Young. Detection of wave fronts in the Indian Ocean from geostationary sunglint satellite imagery. Próxima aparición en el International Journal of Remote Sensing, 2007.

[14] A. Turiel y N. Parga. The multi-fractal structure of contrast changes in natural images: from sharp edges to textures. Neural Computation, 12: 763-793, 2000.

45 [15] A. Turiel y C. Perez- Vicente. Multifractal geometry in stock market time series. Physica A, 322: 629-649, mayo de 2003.

[16] A. Turiel y C. Perez- Vicente. Role of multifractal sources in the analysis of stock market time series. Physica A,

355: 475-496, septiembre de 2005.

[17] A. Turiel, C. Perez- Vicente y J. Grazzini. Numerical methods for the estimation of multifractal singularity spectra 50 on sampled data: a comparative study. Journal of Computational Physics, 216(1): 362-390, julio de 2006.

[18] J. H. van Hateren y A. van der Schaaf. Independent component filters of natural images compared with simple cells in primary visual cortex. Proc. R. Soc. Lond., B265: 359-366, 1998.

Claims

REIVINDICACIONES

1. Procedimiento para el análisis de singularidades en señales digitales de fenómenos físicos, caracterizado porque comprende las siguientes etapas:

a) determinar para cada punto x de la señal un entorno local C que comprende los primeros vecinos de dicho punto x; b) calcular para cada punto x de la señal una medida de singularidad que proporciona una medida de reconstructibilidad de señal en base al entorno local asociado, construida a partir de la inferencia del valor de la

10 señal en dicho punto en base al valor de los puntos de dicho entorno local utilizando la siguiente fórmula de reconstrucción

15 en la que:

-

s es una señal dada,

es el gradiente esencial de s sobre dicho entorno local, es decir, el gradiente restringido al entorno local, cuyos componente son ∂x,Cs(x) = ∂xs(x) para x ∈ C, ∂x,Cs(x) = 0 para x ∉ C, ∂y,cs(x) = ∂yx(x) para x ∈ C, ∂y,cs(x) = 0 20 para x ∉ C, ∂

-

g es el nucleo de reconstrucción universal, es decir, el nucleo vectorial de reconstrucción universal, que en el espacio de Fourier se da por la expresión:

25 y

, donde * se refiere al producto de convolución de funciones,

30 en el que dicha medida de singularidad se basa en la diferencia entre el valor de la señal medida en dicho punto x y el valor de la señal inferido para dicho punto x de dicha fórmula de reconstrucción, estando adaptada la implementación de dicha fórmula de reconstrucción a dicho entorno local limitándola a dicho entorno local.
2. Procedimiento de acuerdo con la reivindicación 1, caracterizado porque también incluye una tercera

35 etapa c) que comprende realizar al menos una transformación logarítmica sobre dicha medida de singularidad que suprime la dependencia de la medida del número de puntos de la señal, obteniendo un exponente de singularidad para cada punto de la señal.
3. Procedimiento de acuerdo con la reivindicación 2, caracterizado porque comprende, antes de

40 realizar la etapa a), obtener una función derivada estable de dicha señal digital que está muestreada a intervalos regulares, y porque en la etapa b) comprende obtener para cada punto de dicha señal digital muestreada una medida de singularidad de la función en ese punto, ponderando las contribuciones de un entorno local del punto y el valor de la derivada en todos los puntos de dicho entorno local.

45 4. Procedimiento según la reivindicación 3, caracterizado porque la derivada estable de la señal digital que precede a la etapa a) se obtiene por derivada de incrementos de un punto hacia la derecha o incrementos centrados de medio punto, estando en ambos casos la derivada definida en el espacio de Fourier por la multiplicación de la señal por los correspondientes núcleos de derivación, donde asumiendo que hay Nx puntos en una dirección coordenada x en la que se quiere derivar, los núcleos de derivación se expresan como se indica a

continuación: Diferencia de un punto hacia la derecha:

n

2π.

N

()∂ = e −1

n

5 Diferencia centrada de medio punto:

.sen π n ; n <N

2

()∂ n = N N −n

−.sen π ; n ≥N

N

comprendiendo las siguientes etapas: 10

-

la aplicación de la transformada de Fourier a dicha señal;

-

la multiplicación de una copia de la transformada de Fourier de la señal por el núcleo asociado a cada una de las d componentes; y

-

la aplicación de la antitransformada a estas d componentes. 15
5. Procedimiento de acuerdo con la reivindicación 2, caracterizado porque la transformación logarítmica de la etapa c) que es al menos una, se realiza como se indica a continuación:

-

para cada punto de la señal digital muestreada se toma la medida obtenida en la etapa b) y se divide por la media 20 de las medidas de todos los puntos; y

-

el logaritmo del resultado se divide por el logaritmo de la escala mínima de la señal digital muestreada, el cual se define como la raíz d-ésima del número total de puntos de la señal, donde d es la dimensión o número de variables propias de la señal.

25 6. Procedimiento de acuerdo con la reivindicación 5, caracterizado porque la medida de singularidad definida en la etapa b) se calcula mediante los siguientes pasos:

-

calcular el vector de entorno de los (2 x d) primeros vecinos de un punto de base x, obteniendo los primeros vecinos del punto x sumando consecutivamente a cada uno y sólo uno de los índices de coordenadas de dicho

30 punto x, primero -1 y después +1, formando el vector de entorno de (2 x d)+1 componentes, cuya primera componente es el valor de la señal en el punto x, la segunda el valor de la señal en el punto obtenido al sumar -1 a la primera coordenada de x, la tercera el valor de la señal en el punto obtenido al sumar +1 a la primera coordenada de x, la cuarta el valor de la señal en el punto obtenido al sumar -1 a la segunda coordenada de x, y así sucesivamente;

35 - extraer la tendencia de este vector, que se define como la suma de sus componentes dividida por ((2 x d)-1) y aplicar esta tendencia al vector de entorno, añadiéndola a la componente referida al punto de base x y sustrayéndola de las otras componentes, de manera que de este modo el nuevo vector de entorno obtenido tiene media nula;

-

aplicar un operador de gradiente local sobre el citado vector de media nula, lo cual devuelve (2 x d)+1 vectores de 40 gradiente cada uno de ellos de d componentes, que definen el gradiente local:

-

anular las componentes de dicho gradiente local asociadas al punto x;

-

aplicar a dicho gradiente local, con las componentes anuladas, un operador de reconstrucción local asociado unívocamente al citado operador de gradiente local obteniendo un vector de (2 x d)+1 componentes, que se denomina señal estimada;

45 - aplicar una vez más el operador de gradiente local a dicho vector de (2 x d)+1 componentes o señal estimada y se obtienen (2 x d)+1 vectores, uno por cada punto del entorno local, de d componentes cada uno, que definen el gradiente local estimado para ese entorno;

-

obtener (2 x d)+1 vectores de d componentes que expresan la diferencia de gradientes entre dicho gradiente local y dicho gradiente local estimado, y - obtener la medida de singularidad asociada al punto x usando estos (2 x d) +1 vectores de diferencia de gradientes.
7. Procedimiento de acuerdo con la reivindicación 6, caracterizado porque dicho operador de gradiente local aplicado al vector de (2 x d)+1 componentes comprende, para cada entorno de un punto x definido por el vector de (2 x d)+1 componentes que incluye el valor de la señal en el punto x y en sus (2x d) vecinos, ejecutar una transformada de Fourier local, que sólo tiene en cuenta este entorno.
8. Procedimiento de acuerdo con la reivindicación 7, caracterizado porque dicha transformada de Fourier local está construida como una matriz de ((2 x d)+1)) x ((2 x d)+1)), cuyos elementos son todos de valor 1 excepto los de la diagonal principal y de las diagonales adyacentes, donde todos los elementos de la diagonal principal excepto el primero por la izquierda, que vale 1, valen la exponencial compleja 2 x π x i/3, donde i es la raíz

15 cuadrada de -1 y los elementos de las diagonales adyacentes valen consecutivamente 1, exponencial de -2 x π x i/3, 1, y así sucesivamente, comenzando desde arriba a la izquierda hacia abajo a la derecha.
9. Procedimiento de acuerdo con la reivindicación 8, caracterizado porque se calcula la transformada

de Fourier local aplicando matricialmente la matriz descrita de ((2 x d)+1)) x ((2 x d)+1)) a un vector de entorno de (2 20 x d)+1 componentes.
10. Procedimiento de acuerdo con la reivindicación 9, caracterizado porque se calcula la antitransformada de Fourier local aplicando la matriz inversa de la descrita en la reivindicación 8, la cual existe siempre.
11. Procedimiento de acuerdo con la reivindicación 10 caracterizado porque para cada entorno de un punto definido por un vector P de (2 x d)+1 componentes, la aplicación a P de dicho operador de gradiente local da un resultado expresado por el vector de gradiente local para el punto x y los puntos de su entorno y

comprende las siguientes etapas: 30

-

aplicar la transformada de Fourier local a P ;

-

construir la derivada a lo largo de una dirección coordenada dada multiplicando por .

la componente del vector

transformada de Fourier de P asociada al punto que se obtiene cuando se modifican las coordenadas del punto x

al sumar -1 al índice de dicha dirección coordenada y multiplicando por − .

la componente de dicho mismo 35 vector obtenida al modificar las coordenadas del punto x al sumar +1 a dicho índice coordenado, y anulando las restantes componentes, obteniendo así d vectores de derivada, uno por cada coordenada;

-

aplicar la antitransformada de Fourier local a estos d vectores de (2 x d)+1 componentes, representando así cada vector la derivada a lo largo de cada una de las d direcciones coordenadas en todos los puntos del entorno local; y

-

reordenar las componentes de estos d vectores, agrupando para cada uno de los puntos del entorno local las d

40 derivadas asociadas a ese punto, obteniendo (2 x d)+1 vectores de gradiente local, de d componentes cada uno, que reproducen el gradiente en cada punto del entorno local.
12. Procedimiento de acuerdo con la reivindicación 6, caracterizado porque dicho operador de

reconstrucción aplicado al gradiente local se define como el inverso del operador de gradiente local, y comprende las 45 siguientes etapas:

-

aplicar la transformada de Fourier local a los d vectores de derivada a lo largo de cada dirección coordenada, cada uno de (2 x d)+1 componentes;

-

construir el vector de reconstrucción a lo largo de una dirección coordenada dada dividiendo por .

la

50 componente del vector transformada de Fourier local de P asociada al punto que se obtiene cuando se modifican

las coordenadas del punto x al sumar -1 al índice de dicha dirección coordenada y dividiendo por − .

la componente de dicho mismo vector obtenida al modificar las coordenadas del punto x al sumar +1 a dicho índice coordenado, y anulando las restantes componentes, obteniendo así d vectores de reconstrucción a lo largo de una dirección, uno por cada coordenada;

-

sumar estos d vectores de reconstrucción; y

-

aplicar la antitransformada de Fourier local al vector de (2 x d)+1 componentes resultante del paso anterior.
13. Procedimiento de acuerdo con la reivindicación 6, caracterizado porque el paso final de la etapa b) 5 por el que se obtiene la medida de singularidad asociada al punto x comprende:

-

retener de los (2 x d) + 1 vectores de diferencia de gradientes obtenidos las d componentes asociadas al punto x, y - obtener la medida de singularidad como la raíz cuadrada de la suma de los cuadrados de estas d componentes

10 con lo que se obtiene una medida de singularidad de correlación local apta para medir la impredecibilidad de un punto dado.
14. Procedimiento de acuerdo con la reivindicación 6, caracterizado porque el paso final de la etapa b)

por el que se obtiene la medida de singularidad asociada al punto x comprende: 15

-

tomar un hipercubo d-dimensional que rodea un punto x dado, formado por los puntos obtenidos al sumar -1, 0 ó +1 a cada índice de coordenadas de x lo que proporciona 3d puntos;

-

retener para cada punto de dicho hipercubo el vector d-dimensional asociado a la componente central (asociada al punto de base) de la diferencia de gradiente, y

20 - sumar estos 3d vectores, y calcular el producto escalar del vector resultante con el vector d-dimensional de diferencia de gradiente asociada al punto x,

con lo que se obtiene un índice de alineación de diferencias de gradiente que permite deducir la existencia de una coherencia espacial entre los errores cometidos al prescindir del punto central cuando la señal es reconstruida, lo 25 cual permite diferenciar entre ruido (de orientación aleatoria) y señal coherente.
15. Procedimiento de acuerdo con la reivindicación 14, caracterizado porque comprende además:

-

antes de realizar la etapa a), obtener una función derivada estable de dicha señal digital que está muestreada a

30 intervalos regulares, y obtener en la etapa b) para cada punto de dicha señal digital muestreada una medida de singularidad de la función en ese punto, ponderando como se indica a continuación las contribuciones de un entorno local del punto y el valor de la derivada en todos los puntos de dicho entorno local; - obtener la energía de gradiente del hipercubo que se ha mencionado anteriormente sumando los módulos al cuadrado de los gradientes de cada punto del hipercubo;

35 - obtener una medida de singularidad de correlación local en el punto x mediante las operaciones siguientes:

-

retener de los (2 x d) + 1 vectores de diferencia de gradientes obtenidos las d componentes asociadas al punto x, y - obtener la medida de singularidad como la raíz cuadrada de la suma de los cuadrados de estas d componentes; y

40 - obtener la medida de singularidad de correlación global como el producto de la medida de singularidad de correlación local por la raíz cuadrada del valor absoluto del índice de alineación de diferencias de gradiente local, dividido este último por la energía de gradiente del hipercubo.
16. Procedimiento según la reivindicación 15, caracterizado porque la derivada estable de la señal digital

45 que precede a la etapa a) se obtiene por derivada de incrementos de un punto hacia la derecha o incrementos centrados de medio punto, estando en ambos casos definida la derivada en el espacio de Fourier por la multiplicación de la señal por los correspondientes núcleos de derivación, donde asumiendo que hay Nx puntos en una dirección coordenada x en la que se quiere derivar dichos núcleos de derivación, los núcleos de derivación se expresan como se indica a continuación:

50 Diferencia de un punto hacia la derecha:

n

2π.

N

()∂ = e −1

n

Diferencia centrada de medio punto:

.sen π n ; n <N

2

()n = N

∂ N −n

−.sen π ; n ≥N

N

5 comprendiendo las siguientes etapas:

-

la aplicación de la transformada de Fourier a dicha señal;

-

la multiplicación de una copia de la transformada de Fourier de la señal por el núcleo asociado a cada una de las d 10 componentes; y

-

la aplicación de la antitransformada a estas d componentes.
17. Procedimiento según la reivindicación 6, caracterizado porque dicha señal digital muestreada es una de:

15 series temporales, transectos de variables físicas seleccionadas entre un grupo que comprende temperatura, concentración de especies químicas, intensidades eléctricas, fuerza, presión, y densidad, en el caso de d = 1; imágenes de entornos reales, tales como imágenes fotográficas, imágenes biomédicas, tales como ecografías, radiografías e imágenes de radiodiagnóstico, e imágenes de medicina nuclear en general, tales como

20 gammagrafías, TAC, PET, RMN, e imágenes obtenidas mediante cualquier otra técnica, imágenes de microscopía de tipo óptica, electrónica y de cualquier otro tipo, imágenes geofísicas, imágenes obtenidas desde satélites y medios aerotransportados, terrestres, sumergidos o de otro tipo, variables distribuidas bidimensionalmente captadas por sensores en tierra, mar, aire, satélite y otros medios en el caso de d = 2; secuencias temporales de imágenes y variables bidimensionales de los casos anteriores en volúmenes

25 tridimensionales en el caso de d = 3; secuencias temporales de variables en volumen en el caso de d = 4.
18. Procedimiento según la reivindicación 6, caracterizado porque la señal digital muestreada en cuestión es representativa de una variable en un fluido turbulento, y porque comprende la obtención de los

30 exponentes de singularidad para la estabilización de dicha variable a fin de realizar un análisis dinámico del fluido, obteniéndose nuevas magnitudes, tales como la difusividad turbulenta de la variable, la viscosidad turbulenta del fluido y otras magnitudes representativas de las escalas no resueltas del fluido.
19. Sistema para el análisis de singularidades en señales digitales de fenómenos físicos, caracterizado 35 porque comprende:

-

medios para obtener para cada punto x de la señal un entorno local que comprende los primeros vecinos de dicho punto x; y

-

medios para calcular para cada punto x de la señal una medida de singularidad proporcionando una medida de

40 reconstructibilidad de señal en base al entorno local asociado, construida a partir de la inferencia del valor de la señal en dicho punto en base al valor de los puntos de dicho entorno local utilizando la siguiente fórmula de reconstrucción

en la que:

-

s es una señal dada,

es el gradiente esencial de s sobre dicho entorno local, es decir, el gradiente restringido al entorno local,

cuyos componente son ∂x,Cs(x) = ∂xs(x) para x ∈ C, ∂x,Cs(x) = 0 para x ∉ C, ∂y,cs(x) = ∂yx(x) para x ∈ C, ∂y,cs(x) = 0

para x ∉ C.

-

g es el núcleo de reconstrucción universal, es decir, el núcleo vectorial de reconstrucción universal, que en el 5 espacio de Fourier se da por la expresión:

y

-

el símbolo · se refiere al producto escalar de convolución, es decir,

, donde * se refiere al producto de convolución de funciones,

en el que dicha medida de singularidad se basa en la diferencia entre el valor de la señal medida en dicho punto x y el valor de la señal inferido para dicho punto x de dicha fórmula de reconstrucción, estando adaptada la

15 implementación de dicha fórmula de reconstrucción a dicho entorno local limitándola a dicho entorno local.
20. Sistema de acuerdo con la reivindicación 19, caracterizado porque incluye además medios para realizar al menos una transformación logarítmica sobre dicha medida de reconstructibilidad que suprime la dependencia del número de puntos de la señal, lo que proporciona un exponente de singularidad para cada punto de

20 la señal.
21. Sistema de acuerdo con la reivindicación 20, caracterizado porque comprende adicionalmente:

medios para obtener una función derivada estable de la señal digital, muestreada a intervalos regulares; y

25 medios para obtener para cada punto de dicha señal digital muestreada una medida de singularidad de la señal en ese punto, ponderando las contribuciones de dicho entorno local del punto y el valor de dicha derivada estable en todos los puntos de dicho entorno local.