Skip to Content

Wavelet Packet Fractal Analysis -
Software Operating Instructions

Author : Cameron L. Jones

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 currently 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 3.2 or 3.3). 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 scanner input of image data. Also available on the Swinburne Server - Download from here!

A conference paper which describes fully the application of this software is available online.

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 for the s.12 wavelet filter. 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. Download the source code to handle the Image Size you require from the table below. Source code is provided in three popular formats. (1) Microsoft Word ver 6.x for Windows format e.g. *.DOC (2) Plain Text file e.g. *.TXT (3) Adobe Acrobat Portable Document Format e.g. *.PDF. Information regarding reading and using files in *.PDF format can be found at the Adobe site. Users of the software should read the Copyright Notice and Permission Notice at the end of this document for information regarding software re-distribution and source code modification.

    Please fill in the FREE online registration form for this software so that I can alert you of software upgrages!

    2-D WPA Fractal Programs
    S-Plus
    Filename    		*.DOC*.TXT*.PDFImage Size
    (pixels)
    wp.frac512wpfr512a.docwpfr512a.txtwpfr512b.pdf512x512
    wpfrac384wpfr384a.docwpfr384a.txtwpfr384b.pdf384x384
    wp.frac256wpfr256a.docwpfr256a.txtwpfr256b.pdf256x256
    wp.frac192wpfr192a.docwpfr192a.txtwpfr192b.pdf192x192
    wp.frac128wpfr128a.docwpfr128a.txtwpfr128b.pdf128x128
    wp.frac3200wpf3200a.docwpf3200a.txtwpf3200b.pdf300x200
    wp.frac64wpfr64a.docwpfr64a.txtwpfr64b.pdf64x64

  3. The program prompts the user to "Enter the path and file name for the *.RAW image file: ". For example if the path is on the C: drive in the C:\photos directory and the file name is *.raw, type c:\photos\*.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. In the example below, the image is called rfa3.raw and is 512x512 = 262144 pixels.

  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. s12
  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. Coefficient energies which are <500 are listed by coefficient energy number and highlighted in red on the sample graph

    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 the 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 find 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). Note that for small image arrays such as 64x64, one can repeatedly test many different wavelet filters by following the onscreen instructions. 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: n (i.e. "No") 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.

    Program File Installation:

    To add the programs to your S-Plus operating environment follow this procedure. The following instructions work in S-Plus for Windows runnning under the MS Windows operating environment. It is also assumed that you have available a text editor such as "Notepad", and another editor or equivalent from a word processing package installed on your system. If you have downloaded *.DOC or *.PDF files, then these have to be converted into *.TXT format in order to be installed within the S-Plus uperating system. For this example, "Notepad" and "Ms Word for Windows" was used. S-Plus prompts are indicated by ">", and "enter" means press the Enter key after typing the command.

    1. Start S-Plus and S+Wavelets
    2. Start MS Word for Windows. Open the program file of your choice. e.g. wpfr512a.txt
    3. In S-Plus, make a new function called "func512". [In the section titled: "Command Line Instructions for Applying 2-D WPA", this function was called: "wp.frac512"]. 
      	>fix(func512)  enter
    4. This loads the default text editor which in this case is "Notepad". You should see the following screen: 
      	function()	
	{	
	}
    5. From Edit choose Select All and then delete with the "delete key" on the keyboard
    6. In MS Word for Windows, choose Select All and Copy. The program text should now be highlighted in black.
    7. In Notepad, choose Edit then Paste. The program text has now been copied across.
    8. In Notepad, choose File then Save. This saves the program to the S-Plus operating system as a new function called: func512
    9. Close Notepad.
    10. Confirm that the program is loaded into S-Plus by typing: 
      	>func512 enter
      You should see the program scroll past in the Command window.
    11. To run the program type: 
      	>func512() enter
    12. Follow the onscreen commands.

    All other programs to handle alternative image file sizes should be installed into S-Plus as appropriate. Note that in the Command Line Instruction section, the program file names are those given in the 2-D WPA Fractal Program table. Therefore wp.frac512() is equivalent to func512() described above.

    Sample Analyses:

    The following links show the actual results which were calculated using the wp.frac512 program. All sample images were 512x512 pixel size. The tables present the Coefficient Index position number, slope, and R squared (r^2) values which were returned. The tables also show which value was used as the cut-off to estimate the fractal dimension. It is clear that the higher order wavelet filters (i.e. s8, s10, and s12) are more efficient at estimating the Fractal Dimension. These filters return an estimated value which is closest to theory. It is stressed that there are no hard and fast rules to follow to ensure that an optimal wavelet filter is used for a particular sample image. The mathematical guidelines indicate that some wavelet filters will be better suited for particular signals than others. For example, if the image is highly symmetrical, then a symmetrical wavelet filter would be more sensible to use, whereas an image with very fine detail might be better approximated with a wavelet filter with several vanishing moments. It is noteable that for the images used here at size: 512x512, the s12 wavelet returned Fractal Dimension estimates which were closest to theory. It is suggested that all images should be analyzed with multiple wavelet filters in order to gain an understanding of how the Fractal Dimension is dependent on wavelet filter architecture, and how this may correlate with specific image features. In this way, it was found that for the sample images described here, the s12 wavelet was the optimal filter in this general situation.

    Table showing results for the d4, d6, and d8 wavelet filters.

    Table showing the results for the s8, s10, and s12 wavelet filters.

    Recent work has focussed on applying 2-D WPA to the analysis of small image arrays (i.e. 64x64 pixels). Contact the author for up-to-date tables and graphs which guide the user in the analysis of small images.

    Current Applications:

    The full set of programs has been tested on selected binary images. Binary images have only 2 grey-scale levels, indicating a filled pixel (black) or an unfilled pixel (white), hence the name binary to indicate the 2 possible states, 1 or 0. Current work in progress is applying 2-D WPA to: (i) Self-affine and multifractal 256 - grey scale elevation or surface images. (ii) Sierpinski carpet binary and grey scale images. (iii) Minimal Spanning Tree graphs. Contact the author for up to date information regarding 2-D WPA software applications to the above image types. These pages will be updated periodically to reflect these changes.

    Sample Images:

    Sample images of 512x512 pixel size in *.RAW image format are available for downloading here. Each file is zipped (*.ZIP) for easy retrieval. In this way, users can replicate the results obtained using the program: wp.frac512 and compare their results with mine (see tables above), before attempting to analyze new images. The following images from left to right are: Fungal colony (Ritz and Crawford, 1990), Neurite (Caserta et al., 1992), Vicsek Snowflake (Vicsek, 1992), Snowflake (Nittman and Stanley, 1987), Plaster dissolution pattern (Lemaire and Van Damme, 1989), Quadratic Kock curve (Mandelbrot, 1982), and a Real Fractal Fungal Colony (Jones et al., 1993).

    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 ]
    • Jones, C.L. (1996) 2-D wavelet packet analysis of structural self-organization and morphogenic regulation in filamentous fungal colonies. Complex Systems Conference - "From Local Interactions to Global Phenomena" July 14-17, 1996 Charles Sturt University; Albury - Wodonga, N.S.W. Australia. Also in: Complex Systems 96 - From Local Interactions to Global Phenomena. (Stocker, R., Jelinek, H., Durnota, B. and Bossomaier, T. eds.) IOS Press, Amsterdam, 1996. pp. 12-23. [Experimental application of the S-Plus software code described here to estimate the global Fractal Dimension of branching fungal colonies; a preprint of this conference paper is available from: http://www.swin.edu.au/chem/bio/cs96/camjones.htm]
    • 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 paper should be referenced as follows:

    TITLE: "Wavelet Packet Fractal Analysis - 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: 11 August 1997 (ver. 4 - updates the previous version ver. 3 [21 June 1996 & 19 April 1996])
    URL REFERENCE: http://www.swin.edu.au/chem/bio/s+code/wpafrac1.htm


    All Web material and S-Plus software code is under copyright © All Rights Reserved. Permission to use the S-Plus code is granted provided that the following information is attached or referenced as appropriate with all distributions and uses thereof:

    Copyright Notice
    Copyright © 1995-1997 by Cameron L. Jones, Program algorithm developed and compiled in S-Plus for Windows (ver. 3.2 & 3.3) by Cameron L. Jones; Centre for Applied Colloid and BioColloid Science,School of Chemical Sciences,Swinburne University of Technology - Hawthorn, Victoria, 3122, AUSTRALIA
    Email: CJONES@swin.edu.au Tel: +613 9214 8187 Fax:+613 9819 0834

    Permission Notice
    Permission to use, copy, and distribute this software and its documentation for any purpose and without fee is hereby granted, provided that the above copyright notice appears in all copies and that both the copyright notice and this permission notice appear in all supporting documentation. Permission to modify the software algorithm or code is not allowed without the prior written consent of the author.

    The author welcomes questions, comments, suggestions, bug reports, etc.

    Document Last Modified: August 11 1997 - Overwrites April 19 1997 version
    [Authors Bibliography]
    Web page production used a combination of HotDog 32-bit v.1.0 and Ant software:
    Web-Counter says you are visitor number