Patent application title: METHOD AND SYSTEM FOR THE ANALYSIS OF PECULIARITIES IN FINGERPRINTS

Abstract:

The proposed method comprises determining for each point of the signal an
environment comprising the first neighbours and calculating for each
point of the signal a reconstructibility measurement, or singularity
measurement, weighting the contributions of said local environment of the
point. The value of the signal at each point is inferred from the values
of the signal at the points of the local environment using a
reconstruction formula. The singularity measurement includes the
difference between the value of the signal at the point and the value
estimated by its local environment. A logarithmic transformation is
carried out on said singularity measurement with a view to obtaining an
independent measurement of the amplitude of the sampled digital signal
that varies in a controlled manner under changes in resolution.
A system is proposed with means for obtaining said singularity measurement
for each point and for carrying out the abovementioned logarithmic
transformation and other calculations.

Claims:

1. Method for the analysis of singularities in digital signals, wherein it
comprises the following stages:a) determining for each point x of the
signal an environment of first neighbours or local environment; andb)
calculating for each point x of the signal a reconstructibility
measurement, or singularity measurement, based on the associated local
environment, constructed from inference of the value of the signal at
said point on the basis of the value of the points of said local
environment using the following reconstruction formulas(x)=({right arrow
over (g)})(x)where:s is a given signal, is the most singular manifold or
MSM of said signal s, is the essential gradient of s,{right arrow over
(g)} is the universal reconstruction kernel, andsymbol • means
convolution scalar productsaid reconstruction formula being adapted to
the abovementioned environment, in such a way that the singularity
measurement contains the difference between the measured value and the
value inferred from the reconstruction formula.

2. Method according to claim 1, wherein it also includes a third stage c)
which comprises carrying out at least one logarithmic transformation on
said singularity measurement, which suppresses the dependency of the
measurement on the number of points of the signal, obtaining a
singularity exponent for each point of the signal.

3. Method according to claim 2, wherein it comprises, before carrying out
stage a), obtaining a stable derivative function of said digital signal
that is sampled at regular intervals, and in that in stage b) it
comprises obtaining for each point of said sampled digital signal a
singularity measurement of the function at that point, weighting the
contributions of a local environment of the point and the value of the
derivative at all points of said local environment.

4. Method according to claim 3, wherein the stable derivative of the
digital signal which precedes stage a) is obtained by derivative of one
point increases to the right or centred half point increases, the
derivative being defined in both cases in the Fourier space by
multiplication of the signal by the corresponding derivative cores, where
assuming that there are Nx points in a coordinate direction x to be
derived, the derivative cores are expressed as follows:Difference of one
point to the right: ( ∂ x ) n = ( 2 π
n N x - 1 ) ##EQU00015## Difference of half point
centred: ( ∂ x ) n = { i sin (
π n N x ) ; n < N x / 2 - i sin
( π N x - n N x ) ; n ≧ N x / 2
##EQU00016## comprising the following stages:application of the Fourier
transform to said signal;multiplication of a copy of the Fourier
transform of the signal by the core associated to each one of the d
components; andapplication of the anti-transform to these d components.

5. Method according to claim 2, wherein the logarithmic transformation of
stage c) which is at least one, is carried out as follows:for each point
of the sampled digital signal the measurement obtained in stage b) is
taken and divided by the mean of the measurements of all points; andthe
logarithm of the result is divided by the logarithm of the minimum scale
of the sampled digital signal, which is defined as the d-nth root of the
total number of points of the signal, where d is the dimension or number
of variables inherent to the signal.

6. Method according to claim 5, wherein the singularity measurement
defined in stage b) is calculated by means of the following
steps:calculating the environment vector of the (2.times.d) first
neighbours of a base point x, obtaining the first neighbours of point x
by adding consecutively to each one and only one of the coordinate
indices of said point x, first -1 and then +1, forming the environment
vector of (2.times.d)+1 components, whose first component is the value of
the signal at point x, the second the value of the signal at the point
obtained by adding -1 to the first coordinate of x, the third the value
of the signal at the point obtained by adding +1 to the first coordinate
of x, the fourth the value of the signal at the point obtained by adding
-1 to the second coordinate of x, and so on successively;extracting the
tendency of this vector, which is defined as the sum of its components
divided by ((2.times.d)-1) and applying this tendency to the environment
vector, adding it to the component referring to base point x and
subtracting it from the other components, in such a way that this new
obtained environment vector has a zero mean;applying a local gradient
operator to said zero mean vector, which returns (2.times.d)+1 gradient
vectors each one of them of d components, which define the local
gradient:cancelling out the components of said local gradient associated
to point x;applying to said local gradient, with the cancelled out
components, a local reconstruction operator unequivocally associated to
said local gradient operator obtaining a vector of (2.times.d)+1
components, which is referred to as an estimated signal;applying once
again the local gradient operator to said vector of (2.times.d)+1
components or estimated signal and obtaining (2.times.d)+1 vectors, one
for each point of the local environment, of d components each one, which
define the estimated local gradient for that environment;obtaining
(2.times.d)+1 vectors of d components which express the difference in
gradients between said local gradient and said local estimated local
gradient, andobtaining using these (2.times.d)+1 vectors of difference in
gradients the singularity measurement associated to point x.

7. Method according to claim 6 wherein said local gradient operator
applied to the vector of (2.times.d)+1 components comprises, for each
environment of a point x defined by the vector of (2.times.d)+1
components which includes the value of the signal at point x and at its
(2.times.d) neighbours, executing a local Fourier transform, which only
takes this environment into account.

8. Method according to claim 7, wherein said local Fourier transform is
constructed as a matrix of ((2.times.d)+1))×((2.times.d)+1)), whose
elements all have a value of 1 except those of the main diagonal and of
the adjacent diagonals, where all the elements of the main diagonal
excluding the first on the left, which is worth 1, are worth the complex
exponential 2.times.π×i/3 where i is the square root of -1 and
the elements of the adjacent diagonals are worth consecutively 1,
exponential of -2.times.π×i/3, 1, and so on successively,
starting from top left towards bottom right.

9. Method according to claim 8, wherein the local Fourier transform is
calculated by matricially applying the described matrix of
((2.times.d)+1))×((2.times.d)+1)) to an environment vector of
(2.times.d)+1 components.

10. Method according to claim 9, wherein the local Fourier anti-transform
is calculated by applying the inverse matrix of the one described in
claim 8, which always exists.

11. Method according to claim 10 wherein for each environment of a point x
defined by a vector {right arrow over (p)} of (2.times.d)+1 components,
application to {right arrow over (p)} of said local gradient operator
gives a result expressed by the local gradient vector for point x and the
points of its environment and comprises the following stages:applying the
local Fourier transform to {right arrow over (p)};constructing the
derivative throughout a given coordinate direction by multiplying by if i
{square root over (3)} the component of the Fourier transform vector of
{right arrow over (p)} associated to the point that is obtained when the
coordinates of point x are modified by adding -1 to the index of said
coordinate direction and by multiplying by -i {square root over (3)} the
component of said same vector obtained by modifying the coordinates of
point x through adding +1 to said coordinate index, and cancelling out
the remaining components, thereby obtaining d derivative vectors, one for
each coordinate;applying the local Fourier anti-transform to these d
vectors of (2.times.d)+1 components, each vector thereby representing the
derivative throughout each one of the d coordinate directions at all
points of the local environment; andreordering the components of these d
vectors, by grouping together for each one of the points of the local
environment the d derivatives associated to that point, obtaining
(2.times.d)+1 local gradient vectors, of d components each one, which
reproduce the gradient at each point of the local environment.

12. Method according to claim 6 wherein said reconstruction operator
applied to the local gradient is defined as the inverse of the local
gradient operator, and comprises the following stages:applying the local
Fourier transform to the d derivative vectors throughout each coordinate
direction, each one of (2.times.d)+1 components;constructing the
reconstruction vector throughout a given coordinate direction by dividing
by i {square root over (3)} the component of the local Fourier transform
vector of {right arrow over (p)} associated to the point that is obtained
when the coordinates of point x are modified through adding -1 to the
index of said coordinate direction and dividing by -i {square root over
(3)} the component of said same vector obtained when the coordinates of
point x are modified through adding +1 to said coordinate index, and
cancelling out the remaining components, thereby obtaining d
reconstruction vectors throughout a direction, one for each
coordinate;adding these d reconstruction vectors; andapplying the local
Fourier anti-transform to the vector of (2.times.d)+1 components
resulting from the previous step.

13. Method according to claim 6, wherein the final step of stage b)
whereby the singularity measurement associated to point x is obtained
comprises:withholding from the obtained (2.times.d)+1 vectors of
difference in gradient the d components associated to point x,
andobtaining the singularity measurement as the square root of the sum of
the squares of these d componentswhereby a local correlation singularity
measurement is obtained suitable for measuring the unpredictability of a
given point.

14. Method according to claim 6, wherein final step of stage b) whereby
the singularity measurement associated to point x is obtained
comprises:taking a d-dimensional hypercube that surrounds a given point
x, made up of the points obtained by adding -1, 0 or +1 to each
coordinate index of x which provides 3.sup.d points;withholding for each
point of said hypercube the d-dimensional vector associated to the
central component (associated to the base point) of difference in
gradient, andadding these 3.sup.d vectors, and calculating the scalar
product of the resulting vector with the d-dimensional vector of
difference in gradient associated to point x,whereby an alignment index
of differences in gradient is obtained which allows deduction of the
existence of a spatial coherence between the errors committed by omitting
the central point when the signal is reconstructed, allowing distinction
between noise (of random orientation) and coherent signal.

15. Method according to claim 14, wherein it comprises additionally:before
carrying out stage a), obtaining a stable derivative function of said
digital signal that is sampled at regular intervals, and obtaining in
stage b) for each point of said sampled digital signal a singularity
measurement of the function at that point, weighting as follows the
contributions of a local environment of the point and the value of the
derivative at all points of said local environment;obtaining the gradient
energy of the abovementioned hypercube by adding the squared modules of
the gradients of each point of the hypercube;obtaining a local
correlation singularity measurement at point x by means of the following
operations:withholding from the obtained (2.times.d)+1 vectors of
difference in gradients the d components associated to point x,
andobtaining the singularity measurement as the square root of the sum of
the squares of these d components; andobtaining the global correlation
singularity measurement as the product of the local correlation
singularity measurement through the square root of the absolute value of
the alignment index of local differences in gradient, the latter divided
by the gradient energy of the hypercube.

16. Method according to claim 15, wherein the stable derivative of the
digital signal that precedes stage a) is obtained by derivative of one
point increases to the right or centred half-point increases, in both
cases the derivative being defined in the Fourier space by multiplication
of the signal by the corresponding derivative cores, where assuming that
there are Nx points in a coordinate direction x where said
derivative cores are to be derived, the derivative cores are expressed as
follows:Difference of one point to the right: ( ∂ x )
n = ( 2 π n N x - 1 ) ##EQU00017##
Difference of half-point centred: ( ∂ x ) n = {
i sin ( π n N x ) ; n < N x / 2
- i sin ( π N x - n N x ) ; n
≧ N x / 2 ##EQU00018## comprising the following
stages:application of the Fourier transform to said signal;multiplication
of a copy of the Fourier transform of the signal by the core associated
to each one of the d components; andapplication of the anti-transform to
these d components

17. Method according to claim 6, wherein said sampled digital signal is
representative of:time series, transects of physical variables selected
from a group comprising, by way of illustration and not limitation,
temperature, concentration of chemical species, electrical intensities,
strength, pressure, and density, in the case of d=1;images of real
environments, which comprises by way of illustration and not limitation,
photographic images, biomedical images (such as ultrasound scans, X-rays
and radiodiagnosis and nuclear medicine images in general, such as gamma
graphs, CAT, PET, MRI, and images obtained by means of any other
technique), microscopic images (optical, electronic and of any other
type), geophysical images, images obtained from satellites or craft
conveyed by air, land, submerged or of any other type, distributed
variables captured in two dimensions by sensors in land, sea, air,
satellite and other media in the case of d=2;image time sequences and
two-dimensional variables of the preceding cases in three-dimensional
volumes in the case of d=3;time sequences of variables in volume in the
case of d=4; orin any number of dimensions, results of numerical
simulations and synthesised signals.

18. Method according to claim 6, wherein the sampled digital signal in
question is representative of a variable in a turbulent fluid, and in
that it comprises obtaining the singularity exponents for the
stabilisation of said variable in order to carry out a dynamic analysis
of the fluid, obtaining new magnitudes such as the eddy diffusivity of
the variable, the eddy viscosity of the fluid and other representative
magnitudes of the fluid's unresolved scales.

19. System for the analysis of singularities in digital signals, wherein
it comprises:means for obtaining for each point of the signal a local
environment that comprises the first neighbours; andmeans for calculating
for each point x of the signal a reconstructibility measurement, or
singularity measurement based on the associated local environment,
constructed from inference of the value of the signal at said point on
the basis of the value of the points of said local environment using the
following reconstruction formula)s(x)=({right arrow over (g)})(x)where:s
is a given signal, is the most singular manifold or MSM of said signal s,
is the essential gradient of s,{right arrow over (g)} is the universal
reconstruction kernel, andsymbol • means convolution scalar
product,said singularity measurement containing the difference between
the measured value and the value inferred for each point.

20. System according to claim 19 wherein it includes also means for
carrying out at least one logarithmic transformation on said
reconstructibility measurement that suppresses the dependency on the
number of points of the signal, which provides a singularity exponent for
each point of the signal.

21. System according to claim 20, wherein it comprises additionally:means
for obtaining a stable derivative function of the digital signal, sampled
at regular intervals; andmeans for obtaining for each point of said
sampled digital signal a singularity measurement of the signal at that
point, weighting the contributions of said local environment of the point
and the value of said stable derivative at all points of said local
environment.

Description:

FIELD OF THE INVENTION

[0001]This invention relates to the analysis of digital signals, in other
words, signals that are sampled at regular intervals, applying wavelet
analysis for implementation of the proposed method, which makes it
possible to identify various subsets of particular points in space to a
scale and position that is more informative about the signal than others.

[0002]The invention likewise relates to a system for the implementation of
the proposed method.

[0003]Throughout this description a digital signal shall be understood as
any structured collection of data uniformly sampled that may be
represented by means of a multidimensional matrix whose positions are
referred to as signal points.

[0004]The invention provides useful techniques and tools for processing,
reconstructing and compressing digital signals based on partial
information regarding their gradient and in particular operating on the
basis of gradient-based measurements obtained through finite increases.
Said techniques and tools can be implemented by means of automatic
algorithms materialised advantageously in computer programs that can be
run in computational environments.

[0005]The invention, given the great efficiency that it provides mainly in
digital signal reconstruction tasks, in particular those representing
images, is applicable in numerous fields, among which as specific
applications one could mention digital signal compression (including
image compression) and the evaluation of flow lines in fluid-related
signals (including determination of current lines in images representing
physical phenomena); and as more general applications, the detection of
structures and recognition of patterns in images of real environments,
such as photographic, geophysical, biomedical and other types of images.

[0006]The invention concerns signals defined in any number of dimensions,
although once the method has been described for a certain number of
dimensions (for example, two), it will be fairly obvious for an expert in
the art to generalise it to signals defined in any number of dimensions.
For this reason, and for the purposes of simplicity, many of the
equations and derivatives presented throughout this description have been
written for 2D, in other words two-dimensional signals, likely to
constitute elements such as images. However, useful results have also
been obtained in other numbers of dimensions and in particular in the
processing of 1D signals, such as the stock market time series (see
references [16], [17]).

[0008]U.S. Pat. No. 6,434,261 describes a method for the detection and
segmentation of digital images in order to locate targets in said images
based on determination of an adaptive threshold to carry out a wavelet
analysis of the digital images that are decomposed into different scale
channels.

[0009]U.S. Pat. No. 7,181,056 concerns a method for the automatic
detection of regions of interest in a digital image representing at least
one portion of biological tissue, wherein a wavelet-based representation
is generated of the regions to be explored.

[0010]U.S. Pat. No. 7,062,085 relates to a method for detecting aspects in
regions of colour images where reference is made to texture
characteristics materialised through coefficients derived from a wavelet
transform based on a multiresolution analysis of the colour digital
image.

[0011]Patent application US-A-2005/0259889 relates to a method for
denoising an X-ray image comprising the application of a complex wavelet
transform to the image bearing the motif, operating with wavelet
coefficients in order to reduce the noise.

[0012]Patent application WO-A-2004/068410 concerns a method for detecting
points of interest in a digital image which implements a wavelet
transform, by associating a subsampled image with an original image.

[0013]The analysis of singularity (see reference [14]) which implements
the concept of describing the local behaviour of a function f(x)
estimated in Rm and defined over Rd around each of its domain
points x in accordance with the so-called Holder singularity exponent, or
Hurst exponent, denoted by h(x), is very useful for many signal
processing tasks, and in particular, is very relevant for compression
purposes and as a pattern recognition tool, and, depending on the
context, can also be used to reveal information regarding the evolution
and dynamic of complex signals.

[0014]U.S. Pat. No. 6,745,129 concerns a method based on wavelets for the
analysis of singularities in seismic data based on the processing of a
representative time series of a record of the phenomenon. The object of
this patent is to calculate the Holder exponent over seismic records
through a continuous wavelet transform. Using this method, when the
signal is analysed (as shown in said patent's FIG. 2b), instabilities
occur that affect both spatial resolution as well as the quality of
determination of the Holder exponent for each point (see discussion of
this issue in reference [11]). This problem in fact makes it impossible
to use the method of U.S. Pat. No. 6,745,129 for tasks of digital signal
reconstruction, unlike the proposals of the method of this invention.
This invention provides a more accurate determination of the singularity
exponents, both in terms of their position, as well as their value. The
difference in accuracy between this invention and U.S. Pat. No. 6,745,129
is due to the use of gradient measurements (which eliminates the
undesirable fluctuations associated to complex wavelets, see reference
[17]) and also because said measurement incorporates an indicator of the
degree of reconstructibility. In accordance with the above, this
invention also makes it possible to reconstruct a high quality signal
based on partial information, unlike the method of U.S. Pat. No.
6,745,129 (see reference [11]).

[0015]Within the field of wavelet based signal analysis, applied in
particular to digital signal processing one of the most commonly used
methods is the so-called Wavelet Transform Modulus Maxima (known as WTMM)
which is determined by the local maxima of the wavelet projections.
Mallat and Zhong (see references [4], [5] and [6]) surmised that this set
can be used to reconstruct the signal completely. Subsequently it has
been verified that the set leads to an attenuated signal and that various
empirical coefficients have to be introduced in order to be able to
reproduce the correct signal amplitudes. Ever since publication of the
document by Mallat and Zhong, there have been several attempts to obtain
high quality reconstruction based on WTMM. In any case, what is most
interesting about the WTMM method is that in the case of images the
highest quantity of lines is concentrated around the borders and
contours; and since it has been known for years (see reference [7]) that
borders and contours contain most of the information about a visual
scene, WTMM has proven itself to be a good candidate for extracting
perceptual information (borders) using an automatic canonical algorithm
based on said method.

[0016]Another branch of research also focused on the use of WTMM was
initiated by Arneodo and collaborators (see references [1] and [2]) who
recognised the capacity of WTMM to deal with multiscale signals,
centering their studies on systems known to present scale invariance
properties, such as turbulent flows and chaotic systems.

[0017]The main inconvenience of all WTMM-based proposals is the
impossibility to systematically extract the maxima when these accumulate
(topologically), a situation that occurs whenever dealing with real
signals, as discussed in [17]. In [10] the problem was studied and
partial solutions were proposed.

[0018]Aside from the use of WTMM, this field of the technique is aware of
recent research by this inventor, A. Turiel, regarding the analysis of
singularities based on gradient measurements. In these works, a gradient
measurement μ associated to a signal s(x) over an arbitrary set A is
defined as the integral of the gradient module over that set:

μ()=dx|∇s|(x) (1)

[0019]This research by A. Turiel shows that when working with real data,
discretised and with noise, it is necessary to operate with wavelet
transforms of the measurements in the following manner: given a wavelet
Ψ, the wavelet transform of gradient measurement μ at a scale r
and at a point x is given by the expression:

[0020]The wavelet transform of measure μ makes it possible to determine
the local singularity exponent, since the dominant term when r is small
depends on r as a power law (see reference [14])

T.sub.Ψμ(x,r)=α.sub.Ψ(x)rh(x)+o(rh(x)) (3)

[0021]By introducing gradient measurements it is possible to improve the
spatial resolution of the singularity exponents (see references [11] and
[14]). In this way, instead of having a singularity exponent every ten
pixels or requiring control of oscillations due to the wavelet, it is
possible to assign a singularity exponent to every pixel with a minimum
dispersion of the point. It is appreciated that there are differences in
the resolution capacities of the different wavelets, meaning that the
need to find an optimised wavelet suitable for dealing with discretised
data has been recognised.

[0022]An important element in the construction of wavelets with a capacity
for optimised resolution is the concept of reconstructing signals using
partial information regarding their gradient. The theoretical approaches
and practical implementations of a gradient reconstruction algorithm
appear in reference [12], which introduces a discussion regarding the
structure of multifractal signals (see references [8], [3] on the
multifractal structure of turbulent fluids). For signals with a
multifractal structure, the set associated to the upper vertex of the
hierarchy is well known, at least from a theoretical point of view, and
is known as the Most Singular Manifold (MSM), which is the set comprising
the points with the most singular (in other words, most negative) values
of h(x).

[0023]Reference [12] mentioned above maintains the theory that MSM
contains sufficient information to reconstruct the signal completely,
analysing the reconstruction of images although the formulae are valid
for any number of dimensions. The reconstruction formula obtained in [12]
for infinitely large domains is as follows:

s(x)=({right arrow over (g)})(x) (4)

where: [0024]s is a given signal, [0025] is the most singular manifold
or MSM of said signal s [0026](x)=(,)(x)) is the essential gradient
vector of s, in other words, the gradient restricted to the MSM, whose
components are (x)=∂xs(x) si x.di-elect cons., (x)=0 si
x(x)=∂ys(x) si x.di-elect cons., (x)=0 si x,
[0027]{right arrow over (g)}(x)=gx(x),gy(x)) is the vectorial
kernel of universal reconstruction, which in the Fourier space is given
by the expression:

[0027] g _ ^ ( k ) = i k k 2 ##EQU00002##

and [0028]notation • means the convolution scalar product, in
other words, ({right arrow over (g)})(x)=(gx*)(x)+(gy*)(x),
where * means the ordinary convolution product of functions.

[0029]From the studies carried out by this inventor (see reference [12])
taking gradient measurements into consideration, it has been concluded
that, if one exists, there is only one possible algorithm for
reconstructing signals using the gradient based on MSM in infinitely
large domains. Likewise, it is known that MSM leads to very good
reconstructions with this algorithm (see references [11] and [12] for
images, and [15], [16] for time series).

BRIEF DESCRIPTION OF THE INVENTION

[0030]This invention proposes a method for the analysis of singularities
in digital signals, which comprises

[0031]a) determining a first-neighbour environment for each point of the
signal; and

[0032]b) calculating for each point of the signal x a reconstructibility
or reconstruction capacity measurement provided by the environment (which
we will refer to, indiscriminately, as "singularity measurement") based
on the abovementioned environment, constructed from inference of the
signal's value at x according to the value of the signal at the points of
said environment of x using the reconstruction function explained in
reference [12] and indicated in formula (4) above, but adapted to the
environment, in such a way that a singularity measurement is obtained
that contains the difference between the measured value of the point and
the value inferred from its environment.

[0033]The proposed method also comprises advantageously a third stage c)
wherein at least one logarithmic transformation is carried out on said
reconstructibility measurement which suppresses the measurement's
dependence on the total number of points of the signal, thus obtaining a
singularity exponent for each point of the signal.

[0034]In an improved embodiment, the proposed method comprises the
following stages:

[0035]a1) obtaining a stable derivative function of a digital signal,
which is sampled at regular intervals;

[0036]b1) obtaining for each point of said digital signal, a singularity
measurement of the function at that point, weighting the contributions of
a local environment of the point (using a reconstruction function as
indicated previously) and the value of the derivative at all points of
said local environment, and

[0037]c1) carrying out at least one logarithmic transformation on said
singularity measurement with a view to obtaining an independent
measurement of the amplitude of the sampled digital signal and that
varies in a controlled manner under changes in resolution.

[0038]The invention also comprises the generalisation of the singularity
measurement as described under the preceding background heading, by
proposing a new singularity measurement based on the set or manifold of
unpredictable points, UPM in its English acronym; in other words, one
starts out by considering the grouping of all unpredictable points, as
against the other points that are predictable.

DETAILED DESCRIPTION OF THE INVENTION

[0039]The object of this invention, as indicated in the previous section,
comprises the calculation of a reconstructibility measurement in order to
calculate accurately the singularity exponents of a digital signal, and
for said exponents to enable obtaining high quality reconstructions. The
basic requirements for defining a singularity measurement μ based on
the unpredictable points manifold UPM are as follows: [0040]i)
Measurement μ must refer to the local singular behaviour of the
functions. [0041]ii) Measurement μ must lead to a more singular
manifold MSM as close to the unpredictable points manifold UPM as
possible.

[0042]The measurements of the unpredictable points manifold are
singularity measurements that also take into account the level of
predictability of the points according to the equation

div()=0 (5)

where F is the UPM and the superscript "c" denotes the complementary set,
in other words, the predictable points; this equation (5) is a
consequence of equation (4) as described in [12]. Therefore, equation (5)
shows that the divergence of the gradient taken only over the predictable
points is cancelled out.

[0043]The inventor proposes herein that the best way of continuing to work
with singularities is to define the measurements based on the set of
unpredictable points as wavelet projections of standard gradient
measurements. In this way, the measurement of the set of unpredictable
points is a wavelet projection of the gradient measurement expressly
designed to penalise unpredictability. This entails generalising the
concept of wavelet projection, with a view to producing wavelet
projections with vectorial values. The use of wavelet projections with
vectorial values has been well known for some time and does not introduce
special complexities into the way of approaching the problem.

[0044]Another main difference in relation to the standard singularity
analysis described in the abovementioned background is that the proposal
of this invention does not carry out wavelet projections of the
singularity measurement on various scales r in order to extract the
singularity exponents by means of a logarithmic regression applied to the
equation (3).

T.sub.Ψμ(x,r)=α.sub.Ψ(x)rh(x)+o(rh(x)) (6)

[0045]One must highlight that projecting the measurement on a wavelet at
various scales is costly in computation time and only serves to improve
the resolution of less singular structures at the expense of worsening
the more singular ones (see the arguments in this respect in reference
[17]). Given that the fundamental objective in relation to this invention
is to extract the most singular structures, it is detrimental to carry
out the projections across multiple scales; instead, the proposal is to
use a point estimator (see references [17], [9] of the singularity
exponents, namely:

where {T.sub.Ψμ(•,(r0)} is the wavelet projection
measurement across the entire signal and serves to decrease the relative
amplitude of the correction o(1/log r0). In applying the equation
(7), r0 must be sufficiently small to disregard this correction. The
scale r0 is defined as the least accessible, in other words, the
scale of one pixel. Conventionally a Lebesgue measurement of 1 is applied
to the entire spatial domain, meaning that, in the case of an image of
N×M pixels, the value of r0 would be set at:

r 0 = 1 NM ( 8 ) ##EQU00004##

thus, in general, sufficiently large images are required to make the first
term on the right hand side of the equation (7) a good approximation to
the singularity exponent. This typically implies having images of at
least 100 pixels in one of the directions.

[0046]An important aspect of this invention lies in the design of digital
wavelets in order to implement singularity measurements based on
unpredictable points manifold. Below two implementations of measurements
based on unpredictable points manifold of this type are presented, which
provide a good result in practical applications. The design is oriented,
overall, at the processing of digital signals and, consequently, the
wavelets are defined (implicitly) by means of numerical weights, although
the presentation is based on a theory and is easy to generalise to a
continuous scheme.

[0047]Another important element of the proposed method lies in the way of
defining and/or establishing numerical estimations of the gradient
∇s so that the reconstruction is numerically stable. To do this,
two possible options are suggested: differences of one pixel or point to
the right and differences of half a pixel, in other words, the
differences in value when moving one position to the right in the first
instance, or the interpolation equivalent to the difference that would be
obtained by moving half a position to the right and to the left of the
point in the second instance. Both are defined by the derivative cores
described in the Fourier space.

[0048]In other words, the stable derivative of stage a1) of the method,
mentioned above, is obtained derived from increases to the right of one
point or centred of half a point.

[0049]In the formulae that follow, ∂x will be
characterised, although the characterisation of ∂y is
analogous. This operator acts on a digital signal by simply multiplying
the signal's Fourier transform by the derivative cores, and then
anti-transforming the result. It is assumed that there are Nx pixels
or points in direction x and Ny in direction y.

[0050]Another basic aspect of the proposed method consists of introducing
the new concept of the cross Fourier transform. In order to estimate the
level of predictability of a given point, a reconstruction formula is
applied, which is expressed by the following equation

s(x)=({right arrow over (g)})(x) (11)

for the least possible number of neighbours of a point, specifically its
first neighbours. In 2D (where d is the signal dimension; therefore, d=2
in this case) it consists of 4 neighbouring points, which, together with
the original point, form a cross. For any quantity p(x) the neighbours of
any point x0 are represented by means of a 5-component vector
comprising said point and its four closest neighbours, following the
indexation convention indicated in the following drawing, which
illustrates in schematic outline the indexation of the points of the
cross in 2D. In this way, the central point will be assigned index 0, the
point to its right index 1, the point to its left index 2, the point
above index 3 and the point below index 4. Thus, the first neighbour
environment of the point under study becomes vector (p0, p1,
p2, p3, p4).

##STR00001##

In other words, in respect of the centre of the cross, the position of all
other points corresponds to shifts of ±1 (in units of one point)
either in direction x or in direction y. In order to define a Fourier
transform specialised or adapted to this cross configuration, one must
taken into account that the basic Nyquist frequency in each direction is
2π/3. In order to simplify the notation, the basic complex element
ζ is introduced:

ζ=ζR+iζI=e2πi/3=cos(2π/3)+i
sin(2π/3) (12)

[0051]According to the invention, the direct cross Fourier transform of
any 5-component vector {right arrow over (p)}=(p0, p1, p2,
p3, p4) is defined as the complex T-component vector {right
arrow over ({circumflex over (p)}=({circumflex over (p)}0,
{circumflex over (p)}1, {circumflex over (p)}2, {circumflex
over (p)}3, {circumflex over (p)}4) obtained from the following
formula:

[0052]This matrix represents the linear combination of the harmonics
associated to the shifts in the cross and is designed to represent as
faithfully as possible the composition in the centre of the cross, based
on the nearest points. The inverse of this matrix can be easily
calculated,

and this last matrix is necessary in order to carry out the anti-cross
Fourier transform.

[0053]It is necessary to define implementations of the gradient and the
reconstruction formula restricted to the cross, in order to rapidly
evaluate the level of predictability of the central point according to
its neighbours. For this reason, this invention proposes suitable
implementations of the gradient and of the reconstruction formula of said
gradient, on the basis of the cross Fourier transform.

[0054]A first implementation is that of the cross gradient operator in
local gradient operator functions, which is operator
(∂x,∂y)=F-1({circumflex over
(∂)}x,{circumflex over (∂)}y)F. In
the Fourier space, said operator acts simply by multiplying any function
by functions {circumflex over (∂)}x and {circumflex
over (∂)}y in order to obtain coordinates x and y,
respectively. Function {circumflex over (∂)}x is
defined for cross environments as follows:

which are defined in such a way as to represent differences of half a
pixel; in fact {square root over (3)}=2 sin(π/3).

[0055]A second implementation is that of the cross reconstruction
operator, which is one of the inverses of the cross gradient operator.
Since the gradient operator eliminates any constant added to each
component of the 5-component vector representing the neighbours, the
reconstruction is fully defined except by a change in that constant; our
proposed implementation of the cross reconstruction operator is such that
the resulting 5-component vector has a zero mean,
Σi=04pi=0. For this reason, the signals must have
the mean subtracted before applying these two operators. To this effect
the elements of the matrix of the first line of the Fourier transform are
applied in inverse cross so as not to introduce harmonics when the
reconstruction operator is applied. Since the sum of the elements of that
first line is (2×d)-1 (whose result is 3 in the case of 2D
signals), in order to subtract the mean all the values of the environment
vector (first neighbours) are added up and divided by ((2×d)-1),
adding to the result the first component of the environment vector and
subtracting the other components.

[0056]The cross reconstruction is the operator R=F-1{circumflex over
(R)}F. In the Fourier space {circumflex over (R)} has two functional
components, {circumflex over (R)}=({circumflex over
(R)}x,{circumflex over (R)}y); the operator acts as the sum of
the product of each component with the corresponding component (x and y)
of the gradient on which it operates. The component {circumflex over
(R)}x is defined for a cross environment as follows:

[0057]In this way, the singularity measurement defined in stage b1) of the
method of this invention can be described by means of the following steps
for a generic signal defined in a space of arbitrary dimension d:
[0058]the environment of the (2×d) first neighbours of a base point
x is extracted, obtaining the first neighbours of point x, consecutively
modifying each one and only one of the coordinate indices of said point
x, first by adding -1 and then by adding +1, forming a vector of
(2×d)+1 components, whose first component is the value of the
signal at point x, the second the value of the signal at the point
obtained by modifying the first coordinate through adding -1, the third
the value of the signal at the point obtained by modifying the first
coordinate through adding +1, the fourth the value of the signal at the
point obtained by modifying the second coordinate through adding -1 and
so on successively; [0059]the tendency of this vector is extracted, which
is defined as the sum of its values divided by ((2×d)-1) and this
tendency is applied to the vector, adding it to the component referring
to base point x and subtracting it from the other components, so that in
this way the new vector obtained has zero mean; [0060]a local gradient
operator is applied to the abovementioned zero mean vector, which returns
(2×d)+1 gradient vectors each one of them of d components, which
define the local gradient: [0061]the components of said local gradient
associated to point x are cancelled out; [0062]the local gradient, with
the cancelled out components, is applied a reconstruction operator
associated unequivocally to the abovementioned local gradient operator
obtaining a vector of (2×d)+1 components, referred to as estimated
signal; [0063]once more the local gradient operator is applied to said
vector of (2×d)+1 components or estimated signal and (2×d)+1
vectors of d components are obtained, which define the estimated local
gradient; -(2×d)+1 vectors of d components are obtained, which
express the difference in gradients between said local gradient and said
estimated local gradient, and from these (2×d)+1 vectors of
difference in gradient the singularity measurement associated to point x
is obtained.

[0064]The cross gradient and cross reconstruction operators are two of the
procedures included in this invention capable of implementation by means
of basic algorithms materialised advantageously in computer programs that
can be run in a computational environment, for the design and calculation
of singularity measurements based on unpredictable point manifolds. In
particular such programs or parts thereof may be included in routines
stored in microprocessors or microchips. These operators can be
simplified to a matrix form of ((2×d)+1)×((2×d)+1), for
faster numerical implementation.

[0065]Below two singularity measurements are described designed in
accordance with the principles of this invention: [0066]local
correlation singularity measurement (lcsm); and [0067]global correlation
singularity measurement (gcsm)

[0068]Both measurements can be implemented by means of specific algorithms
materialised advantageously in computer programs that can be run in a
computational environment. In particular such programs or parts thereof
may be included in routines stored in microprocessors or microchips.

[0069]The local correlation singularity measurement has been conceived in
order to measure the unpredictability of a given point, simply by
calculating the difference between the real value of the signal without
mean (in other words, once the mean has been eliminated) at a given point
and the value inferred on the basis of its four neighbours (when d=2).
The purpose of this measurement is to evaluate
T.sub.Ψlcsmμ(x0,r0) at a given point x0 and in
the case d=2 comprises the following steps: [0070]1. The neighbours of
x0 are converted into a 5-component vector {right arrow over
(s)}=(s0, s1, s2, s3, s4) following the cross
indexation scheme of the drawing shown above. [0071]2. The vector is
appropriately rectified: in the first instance obtaining

pi=s0- S, i=1, . . . , 4 (20) [0072]3. The cross gradient
operator is applied to {right arrow over (p)} in order to obtain vectors
{right arrow over (g)}x and {right arrow over (g)}y. [0073]4.
The value of the components associated to index 0 of said vectors is
preserved for subsequent use, Ax=gx,0, Ay=gy,0.
[0074]5. Said two components are adjusted to zero, gx,0=gy,0=0.
[0075]6. The cross reconstruction operator is applied to the resulting
vectors {right arrow over (g)}x and {right arrow over (g)}y, in
order to obtain the reconstructed signal {right arrow over (r)}. [0076]7.
The cross gradient operator is applied once more to {right arrow over
(r)} in order to obtain the estimated gradients, {right arrow over
(ρ)}x and {right arrow over (ρ)}y. [0077]8. The local
correlation singularity measurement is defined as the module of the
difference of the cross gradients in the centre of said cross, namely:

[0077]T.sub.Ψlcsmμ(x0,r0)= {square root over
((Ax-ρx,0)2+(Ay-ρy,0)2)}{square
root over ((Ax-ρx,0)2+(Ay-ρy,0)2)}
(21) [0078]In fact, this last step means preserving the module of a
wavelet projection with vectorial values, but in order to simplify the
notation it is left as is. [0079]9. Next the singularity exponent
h(x0) is obtained by applying equation (7).

[0080]In other words, the singularity measurement associated to point x
comprises: [0081]withholding from the obtained (2×d)+1 vectors of
local difference in gradient the d components associated to point x, and
[0082]obtaining the singularity measurement as the square root of the sum
of the squares of these d componentswhereupon a local correlation
singularity measurement is obtained suitable for measuring the
unpredictability of a given point.

[0083]The global correlation singularity measurement improves that of
local correlation by taking into account not only the size of deviations
between the estimated and real signals, but also the difference between
the directions of the obtained gradients. For this reason, the initial
data is not only the signal s({right arrow over (x)}), but also the
gradient ∇s({right arrow over (x)}). It is very important to
provide a stable characterisation of ∇s({right arrow over (x)});
to this effect, the two cores mentioned above have been used: a core of
differences of one pixel forward and a core of half pixel increases.

[0084]The global correlation singularity measurement has a more complex
structure; however, the inventor has verified that it is the most
effective for evaluating singularities and at the same time guaranteeing
a high quality of reconstruction. Obtaining this measurement is carried
out in two stages: in the first place, a gradient difference is obtained
for all points; next, the measurement at each point x0 is
constructed by combining the differences in gradient associated to said
point and the gradient ∇s in each group of neighbours of said
point.

[0085]The purpose of this measurement is to evaluate
T.sub.Ψgcsmμ(x0,r0) at a given point x0 and
comprises the following steps:

[0086]First stage: obtaining the differences in gradient at each point
x0. [0087]1. The neighbours of x0 are converted into a
5-component vector {right arrow over (s)}=(s0, s1, s2,
s3, s4) following the cross indexation scheme of the drawing
shown above. [0088]2. The vector is appropriately rectified: in the first
instance obtaining

pi=s0- S, i=1, . . . , 4 (22) [0089]3. The cross gradient
operator is applied to {right arrow over (p)} in order to obtain vectors
{right arrow over (g)}x and {right arrow over (g)}y. [0090]4.
The value of the components associated to index 0 of said vectors is
preserved for subsequent use, Ax=gx,0, Ay=gy,0.
[0091]5. Said two components are adjusted to zero, gx,0=gy,0=0.
[0092]6. The cross reconstruction operator is applied to the resulting
vectors {right arrow over (g)}x and {right arrow over (g)}y, in
order to obtain the reconstructed signal {right arrow over (r)}. [0093]7.
The cross gradient operator is applied once more to {right arrow over
(r)} in order to obtain the estimated gradient vectors {right arrow over
(ρ)}x and {right arrow over (ρ)}y. [0094]8. The
difference in gradient associated to the central point is generated,
(.di-elect cons.x,.di-elect cons.y)=(ρx-Ax,
ρy-Ay).

[0095]Second stage: the global correlation singularity measurement is
evaluated by combining the differences in gradient and the gradients of
the group of neighbours of each point in accordance with the following
steps: [0096]1. For each point x0, the window of 3×3 centred
around it is taken into consideration. In this window, each point has
coordinates x0+(dx,dy), where dx,dy can have the
values -1, 0, 1. [0097]2. The autoprojection of the differences in
gradients in this window is calculated, S(x0):

[0098]In other words, the singularity measurement associated to point x
comprises: [0099]taking the d-dimensional hypercube that surrounds a
given point x, formed by the points obtained when -1, 0 or +1 is added to
the coordinate indices, which provides 3d points; [0100]withholding
for each point of said hypercube the d dimensional vector, the local
difference in gradient of that point, and [0101]adding these 3d
vectors of d dimensions, and calculating the scalar product of the
resulting vector with the d dimensional vector of difference in gradient
associated to point x, which provides an alignment index of local
differences in gradient. [0102]The alignment index of differences in
gradient allows us to deduce the existence of a spatial coherence between
the errors committed by omitting the central point when the signal is
reconstructed, which allows us to differentiate between noise (of random
orientation) and coherent signal. [0103]3. The energy of the gradient
associated to this window is obtained, E(x0):

[0105] T Ψ gcsm μ ( x 0 , r 0 ) = (
x 0 ) S ( x 0 ) E ( x 0 ) ( 26 )
##EQU00013## [0106]In this case, the definition is much more complex
and the linearity has been totally lost, even if wavelet projections with
vectorial values are considered. [0107]One can see that the gradient
energy of said hypercube has been obtained as the sum of the squared
modules of the gradients of each point of the hypercube. [0108]On a
separate note, one can see that the global singularity measurement is the
product of the local correlation singularity measurement obtained as
explained above, from the squared root of the absolute value of the
alignment index of local differences in gradient divided by the gradient
energy of the hypercube. [0109]6. Next the singularity exponent
h(x0) is obtained by applying equation (7).

[0110]The invention described hereto can be implemented by means of
computation techniques run on operating or calculation units. The
implementation of the method comprises a system for the analysis of
singularities in digital signals characterised in that it comprises in a
basic version: [0111]means for obtaining for each point of the signal a
local environment which comprises the first neighbours; and [0112]means
for calculating for each point x of the signal a reconstructibility
measurement, or singularity measurement, based on the corresponding
environment associated to each point, constructed from inference of the
value of each point based on the value of the points of said environment
using the following function or reconstruction formula for the
calculation

[0112]s(x)=({right arrow over (g)})(x)

where [0113]s is a given signal, [0114] is the most singular manifold or
MSM of said signal s, [0115] is the essential gradient of s, [0116]{right
arrow over (g)} is the universal reconstruction kernel, and [0117]symbol
• means the convolution scalar product, said singularity
measurement containing the difference between the measured value and the
inferred value for each point.

[0118]According to an improved embodiment, the system will additionally
comprise means for carrying out at least one logarithmic transformation
on said reconstructibility measurement designed to suppress the
dependency on the number of points of the signal and providing a
singularity exponent for each point of the signal and in general:

[0119]means for obtaining a stable derivative function of a digital
signal, sampled at regular intervals;

[0120]means for obtaining for each point of said sampled digital signal a
singularity measurement of the function at that point, weighting the
contributions of a local environment of the point and the value of the
derivative at all points of said local environment; and

[0121]means for carrying out at least one logarithmic transformation of
said singularity measurement with a view to obtaining an independent
measurement of the amplitude of the sampled digital signal and that
varies in a controlled manner under changes in resolution.

[0122]Said means will comprise in general calculation or data processing
units integrated in the system or materialised in the form of an
integrated circuit or dedicated processing unit. The instructions for
executing the stages of the method will be recorded in programs that can
be loaded on the operating units or integrated in electronic circuits.

[0123]Below several examples are described, to be considered
non-limitative of the application of the method in accordance with the
invention in different fields. All of the digital signals processed in
the following examples have been obtained from public databases, except
those corresponding to FIGS. 5 and 10. All signals have been processed
using a programme written by the inventor in C language and run on a
personal computer with a Linux operating system. The singularity
exponents obtained have been converted into digital images by the same
program.

Example 1

Detection of Structures and Recognition of Patterns

[0124]Singularity exponents make it possible to recognise very subtle
structures which are difficult to detect at first sight. This is so
because the exponents measure the degree of transition of the signal (in
other words, their blur) at each point irrespective of their real
amplitude. This can be used to detect small modifications in a medium and
to prove the existence of new structures in images. The range of
applications covers all manner of images, from medical imagery through to
remote detection, as well as the detection of manipulated photographs.

[0125]FIG. 1 of the drawings indicates the detection of internal oceanic
waves from a MeteoSat satellite image (see reference [14]). The image on
the left shows a portion of the visible channel of the MeteoSat V
satellite acquired on 27 Dec. 2004 over the submarine Mascarene Ridge
(Northeast area of Madagascar); the image has a resolution of 2.5
kilometres×2.5 kilometres approximately, and comprises
500×500 pixels (which corresponds to an area of 1250 km×1250
km). The clouds appear as blurry, white areas, while the sea is the dark
background. The image on the right shows the associated singularity
exponents, represented with a palette that assigns the brightest colours
to the lowest values. It took approximately 10 seconds to obtain said
singularity exponents on a laptop computer with two Centrino processors
at 1.8 Mhz (the times of the other examples refer to the same computer).
In addition to a richer structure associated to the clouds and
atmospheric flows, the existence of concentric oceanic fronts of up to
500 km long was revealed, probably internal waves, in the centre of the
image. We now know that having good knowledge of internal oceanic waves
is key to understanding processes of dissipation of energy and mixture
(of nutrients, spread of contaminants, etc) in oceans; despite this,
there is very little, and not very systematic information regarding the
areas of the planet affected by these waves. For example, those mentioned
in the image on the right, despite their enormous extension (various
fronts of up to 500 km long separated by up to 300 km) had not been
published to date.

[0126]In FIG. 2 the top part shows an image of algae proliferation in lake
Mendota (Switzerland) in false colour, as a combination of various
channels in order to increase the contrast of said algae proliferation.
The bottom part of the figure shows the singularity exponents, which were
obtained following barely 3 seconds of calculation, represented using a
palette of shades of grey in which the lower values are brightest.

[0127]In FIG. 3 the top part shows a view of Alfacs Bay (NE Spain, River
Ebro Delta) registered by band 8 of LandSat, on an unspecified date. The
resolution of this image is 2.5 metres, and the represented zone covers
500×500 pixels. The bottom part of the figure shows the singularity
exponents (calculation time: 10 seconds). Several boats that are barely
visible in the top image appear with well-defined contours in the bottom
image; various wave fronts can also be observed.

[0128]FIG. 4 on the left shows a mammography in digital format with a
resolution of 1976×4312 pixels, extracted from a public archive of
digital format breast scans (USF Digital Mammography Home Page,
http://marathon.csee.usf.edu/Mammography/Database.html) accessed on 10
Oct. 2008. The right hand side of the same figure shows the associated
singularity exponents (calculation time: approximately 4 minutes). The
analysis reveals the structures of the various tissues forming the
breast. This analysis could allow an earlier detection of damage. At the
same time, the capacity to detect singularity lines irrespective of the
contrast makes it possible to reduce the exposure to X-rays required for
the detection of patterns.

[0129]In FIG. 5 the top part shows an image of 200×200 pixels of the
nucleus of an onion cell in interphase, obtained through optical
microscopy (image courtesy of Elisenda Gendra and Monica Pons, Molecular
Biology Institute of Barcelona, Higher Council of Scientific Research).
This image was acquired from a Leica SP1 confocal microscope, in
transmission mode (Nomarski), with argon laser lighting at a 488 nm
wavelength. The bottom part of the same figure shows the associated
singularity exponents (calculation time: approximately two seconds). The
analysis of singularities reveals the existence of coherent lines inside
the nucleus and on its periphery, possibly associated to structures
related to elements of the nucleus such as chromatin and the nuclear
membrane, such structures being difficult or impossible to resolve or
reveal using optical media, in particular in the absence of any form of
staining or marking. The singularity exponents appear to reveal for
example the double membrane structure of the nucleus peripherally, and
also structures associated to the chromatin fibres that fill the nucleus.

Example 2

Image Compression and Denoising

[0130]Due to the reconstruction formula, it is possible to regenerate an
image based on the set of most singular points with enormous quality.
Said set tends to be fairly disperse, constituting 20-30% of the total
points. In order to complete the description, the gradient over said
points must be recorded and stored, and coded in a compact manner. It has
been observed that the gradients change smoothly over the lines of most
singular manifold MSM and it is considered viable to be able to encode
them in a compact manner. Therefore, reconstruction of images based on
the most singular manifold MSM has been proven to have the potential to
provide compression codes for high quality images.

[0131]In FIG. 6, the top part shows an image of van Hateren identified in
reference [18] as imk01020.imc. This image has been obtained using a CCD
camera with a focal distance of 28 mm and is defined by a matrix of
1536×1024 pixels; the data is encoded as shades of grey in 12
nominal bits. It took about 50 seconds to obtain the singularity
exponents. The middle part of the figure shows 30% of the most singular
points. In the bottom part, the image has been reconstructed on the basis
of the gradients over the MSM shown in the middle part, obtaining a
quality measured using the Peak Signal-to-Noise Ratio (PSNR) of 37 dB,
which indicates high quality.

[0132]FIG. 7 shows how reconstruction through MSM makes it possible to
reduce the noise present in the signal. The top part of the drawing shows
the original image (image of Lena, IEEE standard in image processing)
with a resolution of 200×200 pixels; the bottom part, the
reconstruction based on the associated MSM. The contours and borders
contained in the MSM are preserved in the reconstruction, but the
transitions associated to noise, which do not form coherent fronts, are
mostly eliminated, which is particularly noticeable in some areas of the
face, in the reconstructed image.

Example 3

Determination of Flow Lines in Geophysical and Other Types of Images

[0133]Given the theoretical roots on which the definition of singularity
exponents is founded, the latter are very useful when used to analyse
images of scalar variables in turbulent flows. The theory predicts that
singularities are advected (in other words, carried by the fluid), which
can be used in order to trace current lines. In essence, the path of
currents can be followed simply by analysing the images associated to
temperature, chlorophyll concentration and other similar indicators. The
results of the singularities derived from sea surface temperature
detected by microwave sensors (MW SST) on board satellites Modis Acqua
and TRMM are compared with altimeter maps.

[0134]Altimeter data is very difficult to produce and has a very poor
spatial resolution, requiring filtration using a low step filter.
Additionally, to produce quality altimeter maps, various active
altimeters need to be combined, but since 2003 only two satellites remain
in operation, and soon only one of them will be active, or even none.
However, MW SST is much cheaper, can be synoptically obtained over large
zones and is easy to process. As shown in the comparison, the
singularities outline circulation patterns fairly well, demonstrating
that they are channelled by the flow. Therefore, determining currents
using singularities analysis emerges strongly as an interesting
alternative for operational oceanographic systems for the management of
environmental risk.

[0135]FIG. 8 shows on the bottom part the singularity exponents derived
from an image of sea surface temperature obtained from microwaves (MW
SST)-AMSR-E-TMI shown in the upper part, corresponding to 1 Feb. 2003
(image downloaded from Remote Sensing Systems, http://www.ssmi.com/;
calculation time for the analysis of singularities: about 5 seconds). The
zone shown corresponds to the current in the Gulf of Mexico. The
temperature map is given on a cylindrical projection grid with a constant
angular resolution of 1/4 degree. The top of FIG. 9 shows the field of
geostrophic currents obtained by interpolation of four altimeter
satellites, for the same date of 1 Feb. 2003; the bottom part of the
figure shows the overlap of the two fields (the singularity exponents of
temperature of the previous figure and the geostrophic speed field of the
upper panel of this figure).

Example 4

Dynamic Analysis of Variable Sequences in Turbulent Flows

[0136]Nowadays, the numerical simulation of fluids is an essential tool in
tasks such as weather and oceanographic forecasting, the aerodynamic
prototyping of models or various problems of analysis of chemical and
combustion reactions, of industrial interest. However, given the chaotic
nature of fluids in a turbulent regime it is not possible to make an
accurate description limited by a finite number of degrees of freedom
such as those imposed by the discretised grids used in numerical
simulation. This problem is a consequence of the fact that when a fluid
is described using a specific size of grid step the movements which take
place on smaller scales cannot be resolved, which given the chaotic
nature of the fluid cannot be predicted.

[0137]The usual strategy for dealing with these unresolved scales is to
introduce empirical viscosity coefficients (for the speed field) and
empirical diffusivity (for each variable considered in the simulation),
also known as eddy viscosity and eddy diffusivity. These coefficients
represent the more or less random and homogeneous dispersion of the
variable considered in these unresolved scales. Said coefficients serve
to model the effect of unresolved scales on resolved scales through
simulation under certain circumstances; for example, if turbulence is
fully developed or if the integration time of the simulation is
sufficiently large compared to the typical dispersion time of unresolved
scales.

[0138]Determination of the eddy viscosity and eddy diffusivity
coefficients in fluids is fundamental for describing the evolution of
turbulent fluids using numerical models with sufficient quality. A good
determination of these coefficients is crucial in order to obtain greater
accuracy for the abovementioned variables with the numerical simulation,
as well as to extend the time horizon of the validity of these forecasts.
However, in most numerical models used nowadays, these coefficients are
taken as a constant throughout the fluid's domain. This constant is
estimated in a heuristic manner for each numerical execution, although
its value is usually compared with an evaluated experimental value of the
dynamic sequence analysis, according to the following formula:

κ 0 = - 1 2 ∂ t θ 0 2
∇ θ 0 2 ##EQU00014##

where κ0 is the global empirical diffusivity coefficient,
θ0 is the analysed variable, subscript 0 refers to the scale
at which processing is taking place and the triangular parentheses mean
average throughout the spatial domain of the fluid. If instead of
θ0 on takes the current function, what is evaluated is
viscosity, instead of diffusivity.

[0139]In reality, it is possible to pass from a global estimation of
diffusivity as discussed above to a local estimation of this coefficient
for each point of the domain. To do this, the averages of time
derivatives and gradients of the expression above are replaced with
weighted averages with a decreasing weight with the distance from the
point of evaluation. However, if the global evaluation formula presented
above is already somewhat unstable, the local formula is extremely
unstable, giving rise to negative local diffusivity values at certain
points, which is physically unacceptable. Application of the method of
this invention makes it possible to obtain a more stable variable (the
density of the MSM), on the basis of which very stable evaluations of
global diffusivity are obtained and evaluations of local diffusivity that
are not negative at any point.

[0140]As an example, local diffusivities have been evaluated using images
of a colouring agent dispersed in a laboratory experiment. The top row of
FIG. 10 shows a colouring agent dispersed in a 2D turbulent medium in a
laboratory at two moments in time, t=0 s for the column on the left and
t=10 s for the column on the right. (Images of the top row, courtesy of
Patrick Tabeling, Ecole Normale Superieure, Paris). In the middle row of
the same figure, the evaluation of local eddy diffusivity is shown for
all points at the same moments as the images of the top row, using as
variable θ0 the concentration of colouring agent, which is
estimated on the basis of the shade of grey of the figures of said top
row. The values of local diffusivity thereby obtained have been
represented using a palette of two extreme colours (red for negatives,
blue for negatives), with an intermediate colour (white) for values close
to zero. In order to facilitate comprehension of the figure a shaded area
of thick slanted lines has been overlaid on the zones with the largest
size negative values, and a shaded area in a thinner horizontal line on
the zones with highest positive values. As these middle row images show,
the estimation of diffusivity based on concentration gives rise to
extensive areas with negative values; plus, when this sequence is
processed one can see that the determination is very unstable, since the
estimated values of local diffusivity in certain areas change suddenly at
specific moments in time. Finally, the bottom row presents the local
diffusivity evaluations for the same two moments in time obtained based
on the density function of the MSM, which is calculated on the basis of
the singularity exponents evaluated using this invention. As the figure
shows, these evaluations of local diffusivity do not present regions with
negative values; also, observation of the entire sequence shows a gentle
and continuous evolution of the local diffusivity values across all
points.

[0141]Next, a series of references to scientific publications on the state
of the art are included, which reflect aspects explained in this
invention.