function(image)
{
# Set default page/screen width in characters
options(width = 80) # Program credit acknowledgement
cat("WP.FRACTAL program Copyright 1995/1996 Cameron L. Jones\n email: \nCJONES@swin.edu.au \nVersion 1.0 Release 1 for S-Plus (MS Windows \n3.1x)\nWavelet Packet Computation of 2-D objects or image data\nto estimate the \nFractal Dimension from computation of the Hurst exponent.\n")
# Image size requirement
cat("This program operates on images of size 64x64 (nrow,height x \nncol,width)\n\n")
# Wavelet Packet Computation to determine the fractal dimension of 2-D objects.
cat("Enter the path and file name for your *.RAW image file:\n")
# Read input file and select the wavelet filter to use for wavelet packet transform of
# image.
fil <- readline()
fila <- read.image(nrow = 64, ncol = 64, file = fil)
myimage <- fila
cat("Input Wavelet Transform Filter to use for analysis:\n\nChoose from the \nfollowing list:\nDAUBLETS: haar(d2) d4 d6 d8 d10 d12 d14 d16 d18 \nd20\nSYMMLETS: s8 s10 s12 s14 s16 s18 s20\nCOIFLETS: c6 c12 c18 c24 \nc30\n\nYour choice is: ")
ab <- readline()
wavelet.options(wavelet = ab)
# Perform wavelet packet transform and sort wavelet packet coefficient energies
# from highest to lowest.
myimage.wpc <- wp.costs.2d(myimage)
myimage.wpbb <- best.basis(myimage.wpc, data = myimage)
wpbb.sort <- rev(sort(unclass(myimage.wpbb)^2))
# Select subset of sorted coefficient index positions and retain corresponding
# coefficient energies.
wpbb.sorta <- wpbb.sort[c(1, 2, 3, 4, 5, 6, 8, 10, 13,
16, 20, 25, 31, 39, 49, 61, 76, 95, 119, 149,
186, 233, 291, 364, 455, 569, 711, 889, 1111,
1389, 1736, 2170, 2713, 3391, 4096)]
win.graph()
# Examine coefficient energy values and select cut-off for preliminary graph plotting.
y1 <- wpbb.sorta
x1 <- c(1, 2, 3, 4, 5, 6, 8, 10, 13, 16, 20, 25, 31, 39,
49, 61, 76, 95, 119, 149, 186, 233, 291, 364,
455, 569, 711, 889, 1111, 1389, 1736, 2170,
2713, 3391, 4096)
print(y1)
cat("Examine the list of Wavelet Packet Coefficient Energies.\nTake particular \nnote of the Coefficient Index Position number where\nthe energy falls below 1.\nFor \nbinary images of 64x64 pixel size,the 'Fmin'size to\ninclude in the least squares \nregression is the one that returns the highest r^2 value (i.e. closest to 1)\nCoefficient \nEnergies <500 are highlighted in RED.\n")
par(pty = "s")
plot(x1, y1, log = "xy")
label <- y1 < 1
text(x1, y1, label, col = 2, adj = 0)
cat("Enter the number of the last value\nyou want included in the regression: ")
# Read user input.
v <- readline()
mode(v) <- "numeric"
y2 <- y1[1:v]
x2 <- x1[1:v]
# Prepare final version of graphical output.
plot(x2, y2, log = "xy", xlab =
"Sorted Coefficient Index Position", ylab =
"Wavelet Packet Energy Magnitude", type = "n")
points(x2, y2, pch = 18, mkh = 2)
mtext((ab), line = 0.3, adj = 0)
mtext("Wavelet Packet Transform", line = 0.3, adj =
0.345)
# Perform least squares linear regression through data and compute slope exponent.
fit <- lm(log10(y2) ~ log10(x2))
abline(fit, col = 4)
print(summary(fit))
cat("The Hurst Exponent,H is estimated from: (slope)+1=H, and the\nFractal \nDimension,D is determined from: 2-H=D.\nThe slope is the log10(x2) value.\n\n")
repeat {
cat("Do you want to re-compute the regression using fewer data \npoints?\n Type YES=y , NO=n :")
r <- readline()
if(r == "n")
break
else y1 <- wpbb.sorta
x1 <- c(1, 2, 3, 4, 5, 6, 8, 10, 13, 16, 20, 25,
31, 39, 49, 61, 76, 95, 119, 149, 186,
233, 291, 364, 455, 569, 711, 889, 1111,
1389, 1736, 2170, 2713, 3391, 4096)
print(y1)
cat("Examine the list of Wavelet Packet Coefficient Energies.\nTake \nparticular note of the Coefficient Index Position number where\nthe energy falls below \n1.\nFor binary images of 64x64 pixel size, the 'Fmin'size to\ninclude in the least \nsquares regression is the one that returns the highest r^2 value (i.e. closest to \n1 ; see caveat for 64x64 images)\nCoefficient Energies <500 are highlighted in RED.\n")
par(pty = "s")
plot(x1, y1, log = "xy")
label <- y1 < 1
text(x1, y1, label, col = 2, adj = 0)
cat("Enter the number of the last value\nyou want included in the \nregression:")
# Read user input.
v <- readline()
mode(v) <- "numeric"
y2 <- y1[1:v]
x2 <- x1[1:v]
# Prepare final version of graphical output.
plot(x2, y2, log = "xy", xlab =
"Sorted Coefficient Index Position",
ylab = "Wavelet Packet Energy Magnitude",
type = "n")
points(x2, y2, pch = 18, mkh = 2)
mtext((ab), line = 0.3, adj = 0)
mtext("Wavelet Packet Transform", line = 0.3, adj = 0.345)
# Perform least squares linear regression through data and compute slope exponent.
fit <- lm(log10(y2) ~ log10(x2))
abline(fit, col = 4)
print(summary(fit))
cat("The Hurst Exponent,H is estimated from: (slope)+1=H, and \nthe\nFractal Dimension,D is determined from: 2-H=D.\nThe slope is the log10(x2) \nvalue.\n\n")
repeat {
cat("Do you want to decompose the image using another \nwavelet \t\t\t\tfilter?\n Type YES=y , NO=n :")
m <- readline()
if(m == "n")
break
else cat("Input Wavelet Transform Filter to use for analysis: ")
ab <- readline()
wavelet.options(wavelet = ab)
# Perform wavelet packet transform and sort wavelet packet coefficient energies
# from highest to lowest.
myimage.wpc <- wp.costs.2d(myimage)
myimage.wpbb <- best.basis(myimage.wpc,data = myimage)
wpbb.sort <- rev(sort(unclass(myimage.wpbb)^2))
# Select subset of sorted coefficient index positions and retain corresponding
# coefficient energies.
wpbb.sorta <- wpbb.sort[c(1, 2, 3, 4, 5,
6, 8, 10, 13, 16, 20, 25, 31,
39, 49, 61, 76, 95, 119, 149,
186, 233, 291, 364, 455, 569,
711, 889, 1111, 1389, 1736,
2170, 2713, 3391, 4096)]
win.graph()
# Examine coefficient energy values and select cut-off for preliminary graph plotting.
y1 <- wpbb.sorta
x1 <- c(1, 2, 3, 4, 5, 6, 8, 10, 13, 16,
20, 25, 31, 39, 49, 61, 76, 95,
119, 149, 186, 233, 291, 364,
455, 569, 711, 889, 1111, 1389,
1736, 2170, 2713, 3391, 4096)
print(y1)
cat("Examine the list of Wavelet Packet Coefficient \nEnergies.\nTake particular note of the Coefficient Index Position number where\nthe \nenergy falls below 1.\nFor binary images of 64x64 pixel size, the 'Fmin'size to\ninclude \nin the least squares regression is the one that returns the highest r^2 value (i.e. closest \nto 1; see caveat for 64x64 images)\nCoefficient Energies <500 are highlighted in RED.\n")
par(pty = "s")
plot(x1, y1, log = "xy")
label <- y1 < 1
text(x1, y1, label, col = 2, adj = 0)
cat("Enter the number of the last value\nyou want included in \nthe regression: ")
# Read user input.
v <- readline()
mode(v) <- "numeric"
y2 <- y1[1:v]
x2 <- x1[1:v]
# Prepare final version of graphical output.
plot(x2, y2, log = "xy", xlab =
"Sorted Coefficient Index Position",
ylab = "Wavelet Packet Energy Magnitude",type = "n")
points(x2, y2, pch = 18, mkh = 2)
mtext((ab), line = 0.3, adj = 0)
mtext("Wavelet Packet Transform", line= 0.3, adj = 0.345)
# Perform least squares linear regression through data and compute slope exponent.
fit <- lm(log10(y2) ~ log10(x2))
abline(fit, col = 4)
print(summary(fit))
cat("The Hurst Exponent,H is estimated from: (slope)+1=H, \nand the\nFractal Dimension,D is determined from: 2-H=D.\nThe slope is the log10(x2) \nvalue.\n\n")
#Note that when you load this program into S-Plus the formula: (slope)+1=H will
#be re-written with (slope)+1=*H*. This is equivalent to saying take the absolute part
#of the Hurst exponent, and disregard the negative sign.
}
}
}