*=======================================================================
* erodil - find the fractal dimension by a sequence of erosions and dilations
* by Andy Allinger, 2012, released to the public domain
* This program may be used by any person for any purpose.
*
* refer to:
* Petros Maragos and Fang-Kuo Sun
* "Measuring the Fractal Dimension of Signals: Morphological Covers
* and Iterative Optimization"
* IEEE Transactions on Signal Processing v.41 n.1. pp.108-121
* January, 1993
* kindly made available on the first author's web page
*
*=======================================================================
* I/O LIST
* Name Type In/Out Description
* F[n] REAL in sequence of data points
* n INTEGER in # of points
* h REAL out the fractal dimension
*=======================================================================
SUBROUTINE ERODIL(F, n, h)
IMPLICIT NONE
INTEGER n
REAL F, h
DIMENSION F(n)
* local variables
INTEGER i, j, jold, jnew ! indices
REAL
& A(n/2), ! area under curve at each step of dilation
& D(n, 0:1), ! dilated data at prior, current step
& E(n, 0:1), ! eroded data at prior, current step
& L(n/2), ! length scale at each step of dilation
& sum
DOUBLE PRECISION
& x, y, sx, sx2, sy, sxy ! least-squares sums
*=======================================================================
* begin. This version is QUADRATIC complexity in N. To obtain
* the advertised linear time performance -- actually N @ LOG(N) --
* revise this program so that instead of considering each length L
* from 1...N/2, a shortened schedule of lengths L is used, so that L[i]
* is double or triple (ternary logic would be entirely appropriate here)
* the length L[i-1]. The linear time claim seems to assume that i_max
* is a small constant, when it usually is a significant fraction of N.
* Observe the cascade of maxima and minima values. In binary logic,
* comparing D[i-k] and D[i+k], the half-length k (s = 2@k+1) may be doubled
* at each step. In ternary logic, comparing D[i-k], D[i], and D[i+k], the
* length scale may be tripled at each step.
* This speedup is implicit in the discussion in:
* Benoit Dubuc
* "On Estimating Fractal Dimension"
* McGill University Master's Thesis, 9 March 1988.
*
* Dubuc's adviser, Claude Tricot, deserves mention for his work
* on covering algorithms.
*
*=======================================================================
* first step, copy data from F into E and D
E(1,0) = MIN(F(1), F(2)) ! special treatment for endpoints
D(1,0) = MAX(F(1), F(2))
DO i = 2, n-1
E(i,0) = MIN(F(i-1), F(i), F(i+1))
D(i,0) = MAX(F(i-1), F(i), F(i+1))
END DO
E(n,0) = MIN(F(n-1), F(n))
D(n,0) = MAX(F(n-1), F(n))
* length and area at first step
L(1) = 2. / n
sum = 0.
DO i = 1, n
sum = sum + D(i,0) - E(i,0)
END DO
A(1) = sum
* compound the erosions and dilations. alternate between D(i,0) and D(i,1)
DO j = 2, n/2
IF (BTEST(j, 0)) THEN
jold = 1 ; jnew = 0
ELSE
jold = 0 ; jnew = 1
END IF
E(1,jnew) = MIN(E(1,jold), E(2,jold))
D(1,jnew) = MAX(D(1,jold), D(2,jold))
DO i = 2, n-1
E(i,jnew) = MIN(E(i-1,jold), E(i+1,jold))
D(i,jnew) = MAX(D(i-1,jold), D(i+1,jold))
END DO
E(n,jnew) = MIN(E(n-1,jold), E(n,jold))
D(n,jnew) = MAX(D(n-1,jold), D(n,jold))
* length and area
L(j) = 2 * FLOAT(j) / n
sum = 0.
DO i = 1, n
sum = sum + D(i,jnew) - E(i,jnew)
END DO
A(j) = sum
END DO ! next j
* fit a line through the points:
* X[j] = LOG(1 / L[j])
* Y[j] = LOG(A[j] / L[j]^2)
sx = 0. ; sx2 = 0. ; sy = 0. ; sxy = 0.
DO j = 1, n/2
x = LOG(1. / L(j))
y = LOG(A(j) / L(j)**2)
sx = sx + x
sx2 = sx2 + x * x
sy = sy + y
sxy = sxy + x * y
END DO
x = (n/2) * sx2 - sx * sx
y = (n/2) * sxy - sx * sy
h = SNGL(y / x) ! slope of the least-squares line
END ! of ERODIL