Skip to Content

The latest version of this document is available by following this link

Preprint of Wavelet Packet Analysis Fractal-Software Operating Instructions

Step-by-Step Instructions:

Introduction:

The fractal analysis programs which are described here, use two dimensional Wavelet Packet Analysis (2-D WPA) to quantify the global Fractal Dimension, D. Multifractal analysis is not possible using this software. A two dimensional image is analyzed for its frequency content, and an estimate of the scaling correlation between pixels at multiple resolutions and spatial translations is returned. The programs operate within the S-Plus object-oriented programming language. This computer language is available for MS Windows and Unix computer platforms (although the fractal software has been written for the Windows version). The fractal programs also require the S+Wavelets module. This is an add-on to S-Plus but allows 2-D images to be analysed. S+Wavelets requires the 2-D digital image files to be in an appropriate file format in order to be read by the operating language, and in turn be analyzed by the 2-D WPA fractal program. All images must be supplied as a vector of x coordinates, a vector of y coordinates, and a matrix of z values (usually the grey-scale elevation). Standard image file formats such as *.TIFF, *.GIF, *.PCX, *.PICT are not acceptable. Headerless *.RAW files however can be input with ease. Unfortunately many common graphics utility programs such as Adobe Photoshop fail to create a *.RAW file (following conversion from other formats) which is acceptable image input to S-Plus. The following graphics manipulation program solves this problem, and converts between standard image file formats. It is available as shareware from the SimTel Software Repository, located under 'Graphics' Windows software (wlabv30.zip - includes directory information), which includes TWAIN Support for desktop scanners.

Why Wavelets Work:

Wavelets are similar to, but an extension of Fourier analysis and the Wavelet Transform is computationally similar in principle to the Fast Fourier Transform (FFT). The FFT uses cosines, sines and exponentials to represent a signal, and is most useful for representing linear functions. Since many 1- and 2-D signals display non-linear, chaotic or fractal behaviour, Fourier analysis is less suitable for analyzing such signals. Wavelets offer a new method to analyze complex signals. The aim of the wavelet transform is to 'express' an input signal into a series of coefficients of specified energy. The discrete numbers associated with each coefficient contain all the information needed to completely describe the signal provided one knows which analyzing wavelet was used for the decomposition. The wavelet transform partitions a signal with respect to spatial frequency. This is achieved by filtering the signal with a pair of dyadic orthogonal filters called a quadrature mirror filter (QMF). A QMF comes in a pair: termed a 'father' wavelet and a 'mother' wavelet. The father wavelet provides an approximate or blurred version of the signal at successive resolutions, while the mother wavelet captures the detail at each resolution. The wavelet transform of a 2-D signal returns a coefficient matrix which maps all the spatial relationships at multiple scales in the horizontal, vertical and diagonal directions. The coefficient matrix is however a redundant representation, since certain coefficients are large while others are negligible, and practically zero. Wavelet Packets are a subset of the generalised wavelet transform, and offer greater flexibility for the detection of oscillatory or periodic behaviour. This is important since for many 1- and 2-D signals, one may be interested in 'fractal scaling'. In turn, it has been shown that many stochastic (real-world) fractals display either repetitive (persistent) or alternate (antipersistent) scaling trends. This behaviour equates to oscillatory and periodic motion. A statistical technique called the best basis is used to find the optimal subset of coefficients which represent the signal. This subset contains a listing of wavelet packet energies which can be ordered (from highest to lowest) and plotted against index position. A least squares linear regression is plotted through the points to estimate fractal scaling.

The relative accuracy of 2-D WPA has been compared against the Sandbox technique - which is an established technique for measuring the mass-density of 2-D fractal objects. Test object images were 512x512 pixel size. Both methods have been applied to selected test structures of known (theoretical) dimension. The Sandbox method was applied 5 times to each object to ensure good statistics, whereas 2-D WPA was applied once per object. The Sandbox method takes approximately 5 minutes per analysis to collect raw data, with additional graph processing required to arrive at the fractal dimension. In contrast, 2-D WPA takes approximately 2 minutes to perform the transform and another 1 minute to compute the fractal dimension and generate graph output. A table showing the comparative performance of both methods has been prepared. The Sandbox method showed a standard deviation between the theoretical and predicted exponent of ±0.0162. 2-D WPA showed a standard deviation of ±0.0239. Although the Sandbox method showed a lower standard deviation, indicating greater accuracy, the computational time was considerably longer. Also Vicsek (1990) has recently suggested that Sandbox analysis must be performed using multiple origin points to reduce the problem of analyzing an asymmetric structure with symmetric 'circles' or 'square' boxes of known radii. This would add even more time to the computation. In contrast, 2-D WPA takes into consideration all statistical fluctuations in all directions, by analyzing the object in parallel.

Command Line Instructions for Applying 2-D WPA:

It is assumed that S-Plus is resident on the host computer and has been started. The wavelet module (S+Wavelets) needs to be initialised in order for the subroutines to be called from the fractal program. The following steps illustrate the type-written commands which must be entered (manually type in usually followed by ENTER) on the screen at the appropriate 'prompts' to successfully 'run' the described program.
  1. Type: module(wavelets) ENTER
  2. Type: wp.frac512( ) ENTER
    Note that wp.frac512 operates on images of 512x512 pixel size. Programs which work on alternative image sizes are also available, for example:
    2-D WPA Fractal Programs
    NameImage Size (pixels)
    wp.frac512512x512
    wp.frac384384x384
    wp.frac256256x256
    wp.frac192192x192
    wp.frac128128x128
    wp.frac3200300x200


  3. The program prompts the user to "Enter the path and file name for your *.RAW image file: ". For example if the path is on the C: drive in the C:\temp directory and the file name is *****512.raw, type c:\temp\*****512.raw. Press ENTER to read the file. The screen addresses the command, and provides details of the file such as the number of pixels read. If the input file is corrupted or is an incorrect format, an error will be written to the screen.
  4. The user is then prompted to select a Wavelet Transform Filter to use for 2-D WPA from the list. Twenty-two different filters can be individually chosen from a list including: the Haar wavelet, Daublets, Symmlets, and Coiflets. The Haar wavelet is a square wave, Daublets are continuous orthogonal wavelets with compact support and are asymmetric. Symmlets have compact support and are nearly symmetric, while the Coiflets are nearly symmetrical but have vanishing moments for both 'mother' and 'father' wavelet functions - which is not true for the Symmlets.
  5. When prompted for your choice, type in your selection using the same combination of lower case font and number followed by RETURN. e.g. d8
  6. The fractal program now evaluates the input image and transforms it using the chosen wavelet packet filter. This takes approximately 2 minutes depending on processor speed.
  7. At this point a sample graph is plotted showing Wavelet Packet Coefficient Index Position (x-axis) versus Wavelet Packet Energy (y-axis). This provides a rough overview of the scaling behaviour. Underneath this Graphic Screen in the Command Window are shown the individual 'Wavelet Packet Coefficient Energies'. The program prints the following prompt to the screen: "Take particular note of the Coefficient Index Position number where the energy falls below 1." This is a guide since usually values below 1 can be rejected because the end-point coefficient to include is close to but >1. For binary images of 512x512 pixel size [or any other image dimension], the Fmin (minimum resolution of scaling in the image - i.e. usually the thickness of the discrete units which make up the structure) to include in the least squares regression is the one that returns the highest r^2 value (i.e. closest to 1).
  8. Therefore one proceeds systematically and chooses for example 4-5 Coefficient Index Positions to test to find the highest r^2 value. Type in the number at the prompt to indicate which Coefficient Index Position you choose.
  9. A final graph is now plotted in the graphics window. This shows 'Sorted Coefficient Index Position' along th x-axis, versus 'Wavelet Packet Energy Magnitude' along the y-axis. This graph can be cut-and-pasted into the Clipboard, or placed directly into popular wordprocessing packages such as MS Word for Windows.
    In the command window, the least squares regression slope is displayed in a table showing the intercept, and directly below this the slope estimate is printed [log10(x2)]. The standard error is also shown, and the r^2 value is called the 'Multiple R-Squared' value. These values should be written down (especially the slope and the r^2), and another coefficient index position should be tested to detect the highest r^2 value.
  10. The program prompts the user to do this by displaying - "Do you want to re-compute the regression using fewer data points? Type YES=y , NO=n : "
    Just type y for "yes" and n for "no". The list of wavelet packet coefficient energies is again displayed and another number is selected. Again the 'final' graph is re-computed and all statistical data is re-calculated and displayed in the command window.
  11. At this point, the program prompts the user - "Do you want to decompose the image using another wavelet filter? Type YES=y , NO=n : "
    For images of size >384x384 ALWAYS select "no" since the available RAM is usually insufficient to hold ALL wavelet packet coefficient energy tables for more than 1 decomposition (even when onboard RAM is >=16MB). Also, by typing "n", one gets to continue testing for the Fmin value. The program once again prompts the user to re-compute the regression using fewer data points. Exit from the fractal program by typing NO (i.e. "n") at the prompts to perform various functions. Screen and memory usage clean-up is performed by typing q() for QUIT.
  12. When you have found the highest r^2 value, the slope value is used to compute the global Hurst Exponent as follows. The Hurst Exponent, H is estimated from: (slope)+1=H, and the global Fractal Dimension, D is determined from: 2-H=D for 2-dimensional images. This calculation is performed manually.

    Useful references to orient the reader in the experimental application of wavelet analysis:

    • Argoul, F., Arneodo, A., Elezgaray, J., Grasseau, G. and Murenzi, R. (1989) Wavelet transform of fractal aggregates. Phys. Lett. A. 135: 327-336. [2-D signal analysis]
    • Arneodo, A., Bacry, E. and Muzy, J.F. (1995) The thermodynamics of fractals revisited with wavelets. Physica A. 213: 232-275. [Application of the Wavelet Transform which is a predecessor of the Wavelet Packet Transform - 1-D and 2-D signals]
    • Arneodo, A., Argoul, F., Muzy, J.F. and Tabard, M. (1992) Structural five-fold symmetry in the fractal morphology of diffusion-limited aggregates. Physica A. 188: 217-242. [2-D signal analysis]
    • Caserta, F., Hausman, R.E., Eldred, W.D., Kimmel, C. and Stanley, H.E. (1992) Effect of viscosity on neurite outgrowth and fractal dimension. Neurosci. Lett. 139: 198-202. [Test structure]
    • Coifman, R.R. and Wickerhauser, M.V. (1993) Wavelets and adapted waveform analysis - a toolkit for signal processing and numerical analysis (computer program and book) (Massachusetts: A.K. Peters). [Best basis; Application to 1-D signals]
    • Hubbard, B.B. (1996) The world according to wavelets - the story of a mathematical technique in the making. (Massachsetts: A.K. Peters). [General introduction to wavelets]
    • Jones, C.L., Lonergan, G.T. and Mainwaring, D.E. (1993) A rapid method for the fractal analysis of fungal colony growth using image processing. Binary. 5: 171-180. [Test structure]
    • Jones, C.L., Lonergan, G.T. and Mainwaring, D.E. (1996) Wavelet packet computation of the Hurst exponent. J. Phys. A: Math. Gen. 29: 2509-2527. [Application of WPA to 1-D signals, but the mathematics for 2-D images is a simple extension. Full details discussed at length ]
    • Lemaire, E. and Van Damme, H. (1989) Mise en evidence d'une transition vers la fracturation en digitation viscoelastique (From flow to fracture in viscoelastic fingering) Comptes Rendus de L'academie des Sciencs. 309(II): 859-863. [Test structure]
    • Mandelbrot, B.B. (1982) The fractal geometry of nature. (New York: W.H. Freeman and Company) [Seminal work on fractal geometry; Test structure]
    • Nittman, J. and Stanley, H.E. (1987) Non-deterministic approach to anisotropic growth patterns with continuously tunable morphology: the fractal properties of some real snowflakes. J. Phys A: Math. Gen. 20: L1185-1191. [Test structure]
    • Ritz, K. and Crawford, J. (1990) Quantification of the fractal nature of colonies of Trichoderma viride. Mycol. Res. 94(8): 1138-1152. [Test structure]
    • Vicsek, T. (1990) Mass multifractals. Physica A. 168: 490-497. [Problems with Sandbox method]
    • Vicsek, T. (1992) Fractal growth phenomena (2nd. ed.) Singapore: World Scientific Publishing Co. [Test structure]
    • Wickerhauser, M.V. (1994 ) Adapted wavelet analysis from theory to software. (Massachusetts: Peters). [Best basis; Application to both 1-D and 2-D signals]

    This preprint should be referenced as follows:

    TITLE: "Preprint of Wavelet Packet Analysis Fractal-Software Operating Instructions"
    AUTHOR: Cameron L. Jones
    AFFILIATION: Centre for Applied Colloid and BioColloid Science, Swinburne University of Technology, School of Chemical Sciences, GPO Box 218, Hawthorn, Victoria, 3122, Australia
    EMAIL: CJONES@swin.edu.au
    DATE: 28 May 1996 (ver. 2 - updates the previous version)
    URL REFERENCE: http://www.swin.edu.au/chem/bio/splus/splusfaq.htm