Image Anal Stereol 2011;30:63-76
Review article
GENERALIZED FRACTAL TRANSFORMS AND SELF-SIMILARITY: RECENT RESULTS AND APPLICATIONS
Davide La Torre^ and Edward R. Vrscay^
^Department of Economics, Business and Statistics, University of Milan, Italy; ^Department of Applied Mathematics, University of Waterloo, Waterloo, Ontario, Canada N2L 3G1 e-mail: davide.latorre@unimi.it, ervrscay@uwaterloo.ca (Accepted June 7, 2011)
ABSTRACT
Most practical as well as theoretical works in image processing and mathematical imaging consider images as real-valued functions, u : X ^ R^, where X denotes the base space or pixel space over which the images are defined and r^ C m is a suitable greyscale space. a variety of function spaces .^{X) may be considered depending on the apphcation. Fractal image coding seeks to approximate an image function as a union of spatially-contracted and grey scale-modified copies of subsets of itself, i.e., u« Tu, where T is the so-called Generalized Fractal Transform (GFT) operator The aim of this paper is to show some recent developments of the theoiy of generahzed fractal transforms and how they can be used for the purpose of image analysis (compression, denoising). This includes the formulation of fractal transforms over various spaces of multifunctions, i.e., set-valued and measure-valued functions. The latter may be useful in nonlocal image processing.
Keywords: fractal transforms, iterated function systems, measure-valued functions, multifunctions, nonlocal image processing, self-similarity.
INTRODUCTION
In his classic work. The Fractal Geometry oT Nature, Mandelbrot (1977) presented the first description, along with an extensive catalog, of selT-slmllar sets, namely, sets that may be expressed as unions of contracted copies of themselves. He called these sets "fractals," because their (fractional) Hausdorff-Besicovitch dimensions exceeded their (integer-valued) topological dimensions. The ternary Cantor set and the von Koch "snowflake curve" are two of the most famous examples of such sets.
Hutchinson (1981) and, shortly thereafter, Bamsley and Demko (1985) showed how systems of contractive maps with associated probabilities, referred to as "iterated fiinction systems" (IPS) by the latter, can be used to construct fractal, self-similar sets and measures. These sets and measures are attractive fixed points of Tractal transTorm operators. (We shall briefly review IPS in the next section.) But Bamsley and Demko were the first to see the potential of using IPS for the purpose of approximation-. Given a "target" self-similar set (or measure), say 5, find an IPS fractal transform operator T with fixed point 5 that is as close as possible to 5. More on this below.
The formulation of IPS-type methods over various complete metric spaces has been an ongoing research programme. It involves the construction of appropriate IPS-type operators, or generalized
Tractal transTorms (GPT), over these spaces, including various fianction spaces and distributions (Cabrelli ez-aZ, 1992; Porte and Vrscay, 1998a;b), vector-valued measures (Mendivil and Vrscay, 2002), integral transforms (Porte e? a/., 1999), wavelet transforms (Mendivil and Vrscay, 1997; Vrscay, 1998). More recently, we have formulated GPTs over set-valued fianctions and measures, i.e., multifianctions, e.g., Kunze e?a/. (2007; 2008); La Torre and Mendivil (2008); La Torre e? a/. (2009a;b); La Torre and Mendivil (2009).
The action of a generalized fractal transform T : on an element u of the complete metric space {ß^{X),d) can be summarized in the following steps:
1.	It first produces a set of N spatially-contracted copies of u.
2.	It then modifies the values of these copies by means of a suitable range-mapping.
3.	Pinally, it recombines these altered copies by means of an operator appropriate to the space
to produce the element ve	i.e., v =
Tu.
Under conditions appropriate for each space, the generalized fractal transform T is a contraction mapping which, by Banach's fixed point theorem, guarantees the existence of a unique fixed point u = TU.
Most practical as well as theoretical works in image processing and mathematical imaging consider images as real-valued functions. There are, however, situations in which it is useful to consider the greyscale value of an image u at a point x as a random variable that can assume a range of values R^ C M. This is an example of a multifunction representation of image functions. But it is often not enough to simply know the greyscale values that may be assumed by an image u at a point x. one must also have an idea of the probabilities (or frequencies) of these values. As such, it may be more useful to represent images by measure-valued functions, for example, ß ^{Rg), where ^{Rg) is the set of probability measures supported on Rg (La Torre et af, 2009b). This is another example of multifiinction representation of an image. Later in this paper, we outline this formulation along with an appropriate class of fractal transforms acting on this space.
The IFS-based inverse problem, which has become important in a number of applications, may then be phrased as follows:
Given a "target" element v e -^{X), find a (point-to-point) contraction mapping T : .^{X) ^ .^{X) with fixed point ü such that c/(v; u) is as small as possible.
From a practical perspective, however, it is difiicult to construct solutions to this problem so one relies on the following simple consequence of Banach's fixed point theorem, known in the fractal coding literature as the
Collage Theorem (BsmsXcy et af, 1985):
Theorem 1 For any v G .^{X), where c is the contractivity factor ofT.
(1)
Instead of trying to minimize the error d{v,u), one looks for a contraction mapping T that minimizes the so-called collage error d{ v, Tv). As we shall describe below, this is the essence of fractal image coding (Fisher, 1995; Lu, 2003). However, this method of collage coding may be applied in other situations where contractive mappings are encountered. We have shown this to be the case for inverse problems involving differential equations. In the simplest case of ordinary differential equations, the contractive mapping is the Picard integral operators associated with the initial value problem (Kunze and Vrscay, 1999).
At this point, it should be mentioned that in collage coding, the contractive (fractal) transform T is
generally defined in terms of a finite set of parameters. In fractal image coding, this set is often referred to as the fractal code associated with the image. Solving the inverse problem using collage coding is based on the following continuity property of fixed points of contractive mappings (Centore and Vrscay, 1994):
Theorem 2 Let {J^{X),d) be a complete metric space and Con{^ {X)) a set of contraction mappings T : .^{X) ^ .^{X). Let_Ti,T2 G_Con{.^{X)) with respective ßxed points, u\ and U2 and contraction factors ci and C2. Deßne the distance between T\ and Ti as follows,
dconiX){TuT2)= sup d{71 U, T2u) .

Then
(l{ui,U2) < —--dcon{S^{X)){T\, T2) ,
^ Qnin
where c^m = min(ci ,02).
ITERATED FUNCTION SYSTEMS (IPS)
IFS: Here we briefly review the IFS formalism of Hutchinson (1981) and Bamsley and Demko (1985). In what follows, {X, d) denotes a compact metric "base space" (or "pixel space"), typically [0,1]". Let w = {ivi, • • • , wn} be a set of 1 -1 contraction maps Wj-.X^ X, to be referred to as an A'-map IFS. Let q G [0,1) denote the contraction factors of the Wi and define c = maxi<KA^Q. Note that cG [0,1).
Now let ^{X) denote the set of nonempty compact subsets of X and dt the Hausdorff metric. Then c//,) is a complete metric space (Hutchinson, 1981). Associated with the IFS maps Wi is a set-valued mapping w : Jif{X) Jif{X) the action of which is defined to be
n

(2)
i=i
where Wi{S) := {wi{x),xe 5} is the image of 5under wj,i= 1,2, •••,7V.
Theorem 3 (Hutchinson, 1981) w is a contraction mapping on {Jif{X), c//,) .•
dt{^{A)MB))<cdt{A,B), A,BeJf{X). (3)
Corollary 1 There exists a unique set A G Jf(X), such thatw{A) = A, the so-called attractor of the IFS w. This implies that
A=(]wlÄ). i=i
In other words, the attractor A is self-similar since it maybe expressed as a union of copies of itself
Moreover, A is globally attractive; For any S G .^{X), dt{w{S),A) ^ 0 35/3 ^
Simple examples:
1.	[0,1] and Af = 2, with wi(x) = ^x, W2{x) =
Then the attractor A is simply [0,1],
2.	X= [0,1] and Af = 2, with wi(x) = ^x, W2{x) =
+ I. Then the attractor A = C, the classical ternary Cantor set on [0,1].
3.	[0,l]2andAf=3withIFSmaps,
W2{x,y) =	,
/ X /1	11
The attractor is the "Sierpinski gasket" shown in Fig. 1 belnw
Fig. 1. "Sierpinskigasket"
IFSP: Now let ^{X) denote the set of Borel probability measures on X and du the Monge-Kantorovich metric on this set (referred to as the "Hutchinson metric" in the IFS literature):
dH{ß,v)= sup
feLipi (XjR)
JX
f{x)dß- f f{x)dv J X
where
I f{x,) - f{x2)\ < d{x,,X2),\lx,,X2 G X} .
For 1 < 7 < A'^ let 0 < pj < 1 be a partition of unity associated with the IFS maps Wj, so that E^i Pi = 1. Associated with this N-map IFS with probabilities (IFSP) (w, p) is the so-called Markov operator, M: J^{X) J^{X), the action of which is
n
v{S} = {Mß){S} = J^p,ß{wT\S}), \/Se.^{X). i=i
(4)
Theorem 4 (Hutchinson, 1981) M is a contraction mapping on {^{X), du)
dH{Mß,Mv)<cdH{ß,v),	(5)
where c = maxi<K#Q is the contraction factor
Corollary 2 There exists a unique measure ß G ^{X), the so-called invariant measure of the IFSP (w,p), such that ß = Mß. Moreover, for any ß G ^{X), dH{M"ß,ß)^Oasn^^.
Simple examples:
1.	The 2-map IFS in Example 1 above, with attractor A = [0,1]. When pi = p2, the invariant measure ß is Lebesgue measure on [0,1]. A histogram approximation of the invariant measure for the case Pi = 0.4, p2 = 0.6 is presented in Fig. 2.
2.	The 2-map IFS in Example 2 above, with attractor 4 = C, the classical ternary Cantor set on [0,1]. When Pi = p2 = ß is the classical Cantor-Lebesgue (uniform) measure supported on C.
The reader is referred to Bamsley (1989) for more detailed discussions as well as numerous examples.
In applications, it is most convenient to employ affine IFS maps. In this case, the moments of the invariant measure ß of the Markov operator M satisfy a set of relations that allow them to be computed recursively (Bamsley and Demko, 1985; Bamsley, 1989; Forte and Vrscay, 1995). We illustrate with the one dimensional case, i.e., X= [0,1]. The extension to higher dimensions is quite straightforward.
Fig. 2. Histogram approximation to invariant measure on [0,1] for Example 1 above.
The affine IFS maps will be denoted as follows,
Wi{x) = SiX+ at, 7 = 1,2, • • • , A^.	(6)
We consider the moments of a probability measure ß e ^{X) defined as follows,
gn= f X"dß, /3=0,1,2,--- . J X
(7)
By definition, go = \. Now let V = Mfx. Then, from Eq. 4, the moments of V are given by
hn= j'
Jx
. n JX i^i
Expansion of the binomial followed by an interchange of summation and integration yields the result
J=o Vi,
■ n
J, Pi si a';-i=l

(8)
If we let
g = {go,gi,---V, h = {ho,hi,-
denote the (infinite) moment vectors of ß and V, respectively, then the Markov operator M is seen to induce a linear mapping h = Ag, where A is represented by a lower triangular matrix. This was originally pointed out in Forte and Vrscay (1995).
In fact, the linear operator A is contractive in the following complete metric space of weighted fi moment vectors (Forte and Vrscay, 1995),
^ = = = i|<-}. (10)
The unique fixed point of this operator is the moment vector belonging to the invariant measure ß = Mß of the IFSP. In this special case, the moment vectors g and h in Eq. 9 are equal, i.e., h„ = g„. Eq. 8 can then be rearranged to yield
n
\
i=l
po \JJ
' n
^Pisja"'^ i=i
S J
(11)
This result, originally derived in (Bamsley and Demko, 1985), shows that the moments gn of the invariant measure ß may be computed recursively, starting with
go = 1-
The fact that A is contractive naturally leads to a collage theorem for moments (Forte and Vrscay, 1995). This leads to a formulation of the inverse problem of IFSP-based approximation of measures in terms of moments.
Before concluding this section, we mention that in the historical development of IFS, measures (hence the method of IFSP) were viewed as being potentially more usefial for the representation/approximation of images, because oftheir ability to accomodate shading. As a result, a good deal of work was devoted to the inverse problem of approximation of measures using IFSP - see, for example, (Vrscay and Roehrig, 1989; Vrscay, 1990). In fact, the fractal block image coding method of Jacquin (1992) was originally formulated in terms of measures (Jacquin, 1989), although it is also quite naturally expressed in terms of fianctions, as will be seen below.
IFSM: In this setting we consider a general Sanction space	supported on X. The essential
components of a fractal transform operator are as follows.
1. A set of A'^one-to-one contraction maps Wi\X-with the condition that U^j Wi{X) = X.
X
(9) 2. A set of associated ^re/sca/e ma/35 : R ^ R that are assumed to be Lipschitz on R, i.e., for each there exists aKj>0 such that
IUti) - Hb)I < Kj\ti-t2\, for all ti,t2e R.
In most application, the greyscale maps are assumed to be affine, i.e..
({)i{t) = ait+ßi,
(12)
which automatically condition.
satisfies the Lipschitz
The above two sets of maps are said to comprise an "Iterated Function System with greyscale maps" (IFSM), denoted as (w, O) (Forte and Vrscay, 1998a). For each x G A', this IFSM produces one or more fractal components defined as

0,	otherwise.
If several fractal components exist for an x G X, then they are combined with an operation that is suitable for the space in which we are working (see Forte and Vrscay (1998a) for more details and examples of the various fianction spaces that can be considered). We usually consider the summation operator for .^{X) = LP{X), i.e., for a u G LP{X), the action of the fractal transform is given by
n
v{x) = {Tu){x) = J^gix)
(13)
i=i
Theorem 5 (Forte and Vrscay, 1998a) Let (w, O) be an IFSM as deßned above, with spatial contractions Wi and Lipschitz greyscale maps (pj. Then for p> I andu,vGLP{X),
Tu- Tv\\<
Corollary 3 If c =
w
ley PR, i=i

u— V
(14)
< 1, then T is
contractive in LP{X) with fixed point u G LP{X). The fixed point equation.
n
i=l

indicates that u is "self-similar," i.e., that it can be written as a sum of spatially-contracted and greyscale-modißed copies of itself
Example: [0,1] and 3, with IFS maps
/ X	1	/ X	1	1	/ X	1	2
m (X) = -X, W2(X) = -X+-, W3(X) = -X+-,
and associated (pi^ maps,
(Note that in the I^-sense, the subsets Wi{X) may be considered as nonoverlapping.) The fixed-point Sanction u{x) of this IFSM is the famous "Devil's staircase fianction," sketched in Fig. 3 below. Clearly, ü{x) may be viewed as a union of three contracted copies of itself, with the middle copy being a "flattened" one.
Fig. 3. "Devil's staircase function " on [0,1],
It is also convenient to define IFSM operators with condensation functions. For example, given a set of IFS maps Wi, associated constants «y and condensation fianction b{x),xe X, define the action of the associated operator T as follows: For u^Ü{X),
n
v{x) = {Tu){x) = b{x) + Xaiu{w7\x)) . (15)
i=i
We now have the apparatus to consider the inverse problem of IFS-based approximation of Sanctions. In practice, one normally works with a fixed set of IFS maps Wi, I < i< N, and then finds the optimal associated greyscale maps - optimal in the sense that the collage distance || v— Tv\\ is minimized, where I'is the fianction to be approximated. This is the basis of fractal image coding, which we outline in the next section.
LOCAL SELF-SIMILARITY AND BLOCK FRACTAL IMAGE CODING
In practical applications, it is overly ambitious to expect that a signal or image will display the self-similarity property used above, i.e., that it can be well approximated as a union of spatially-contracted and range-modified copies of itself It is more reasonable to expect that a signal or image be locally self-similar, i.e., it may be well approximated as a union of spatially-contracted and range-modified copies of subsets of itself This is the basis of Jacquin's original fractal block coding method (Jacquin, 1989; 1992) which is also known as the local or partitioned IFS method (Bamsley and Hurd, 1993). We forego a formal mathematical discussion of this method and simply consider the particular case of fractal block coding of images. Here, subblocks of an image are
approximated by contracted and greyscale-modified copies of other subblocks of the image.
A very simple prescription for the fractal coding of an /3 x /3-pixel image u{x) is as follows. Let Rk, k =	,Nr, denote a set of hr x riR-pixel
nonoverlapping range blocks that form a partition of the image. Let Dj, k = 1,2, • • • , Afe be a set of IriR x IriR-^hieX domain blocks that are selected from throughout the image. (In order to keep the size of the domain pool down, but at the expense of some accuracy, we may consider the set of nonoverlapping InR x 2/3;e-pixel blocks that cover the image.)
For each range block, compute the collage errors Ajtj associated with all domain blocks, Dj, i.e.,
A,J = mm II u{R,) - au{Dj) -ß\\, J = 1,2,-■■ ,Nd.
a,ß
(16)
Here, u{Dj) denotes the nR x nR-pixel block image obtained by "decimating" the InR x 2/3;e-pixel domain block image u{D]f). (For digital images, decimation is normally accomplished by replacing the image values over four neighbouring pixels that form a square in D]f by their average value placed on one pixel. Following the decimation, we may consider all eight possible isometries that map one block to another, i.e., four rotations and four reflections.) The block L'j(jt) yielding the lowest collage error Ajtj is chosen to be the domain block associated with Rj.
The above procedure yields a fractal transform T which is defined in terms of the range-domain assignments {kj{k)) (along with isometries i{k) if applicable) and ^-map parameters a^^ßk- These parameters comprise the fractal code of the image u. The action of T may be expressed as follows: For each range block Rk,l<k< Nr,
=	xeRk- (17)
By construction, the fractal transform T minimizes the total squared collage distance
Nr

k=\
%j{ky
(18)
over the nonoverlapping range blocks Rk. (Because the range blocks Rk are nonoverlapping, each approximation can be performed independently.)
The fixed point u of T - the desired approximation to u - is then generated by iteration: Start with any n x /3-pixel "seed" image, uq, for example the blank image uq = 0, and form the iteration sequence u^+i = Tum. (Because the image is discrete, convergence is achieved in a finite number of iterations.) At each step
/33 > 1 of the iteration procedure, each range block image u„,{Rk) of u„, supported on Rk is replaced by the affine scaled image akU„,{Dk) + ßk-
The result of this procedure, as applied to the 512 x 512-pixel, 8 bit-per-pixel test image Boat, is shown in Fig. 4. Here, the range blocks Rk employed in the calculation were the 4096 nonoverlapping 8x8-pixel blocks of the image. The domain blocks D k used were the 1024 nonoverlapping 16 x 16-pixel blocks. The top left of the figure shows the original test image. Moving clockwise in the figure, the iterates u\ and U2 corresponding to the "seed" image uo = 0 (black) are shown. The lower left image is the fixed point u = uio corresponding to the fractal transform T.
In this example, there was no attempt to perform any image compression. As such, the a and ß parameters were stored as real numbers to full machine precision, and not quantized according to any prescribed bit allocation. The so-called "PSNR value" of the fixed-point approximation, a measure of the accuracy of the approximation in terms of L^ error is roughly 25 dB. (The higher the PSNR, the lower the L^ error.) A better approximation to the test image, corresponding to a higher PSNR value, would be achieved if 4 x 4-pixel blocks were used for the range blocks Rk.
For more detailed accounts of fractal coding, the reader is referred to Bamsley and Hurd (1993); Fisher (1995); Lu (2003).
GENERALIZED TRANSFORMS
Iterated Function Systems on Multifunctions:
We now outline a simple IFS-type method on multifianctions, that is, set-valued functions. As a motivation, we suppose that to each pixel of an image is associated an interval which measures the "error" in the greyscale value at that pixel. For simplicity, we examine only the one-dimensional case where the base space is A' = [0,1]. The extension to higher dimensions is straightforward.
Consider the space of multifianctions ß^ = {F : X Jf{Y)} and suppose that F{x) is a compact subinterval of F for all xGX.lt is quite straightforward to prove (La Torre et al., 2009a) that is a complete metric space with respect to the following metrics:
cL(F,G)=supc4(F(x),G(x))
AreX
and
dp{F,G) = ( f dt{F{x),G{x))Pdß{x) \.Jx
\ i/p y
Fig. 4. Starting at upper left and moving clockwise: The original 512x512-pixel, 8 bit/pixel "Boat" test image. The iterates u\ and U2 along with the ßxed point u = u\q of the fractal transform operator T designed to approximate the test image. The "seed " image was uq =0 (black).
where ß is the Lebesgue measure. A fractal transform operator on ^ may now be defined in terms of the following ingredients:
1.	As before, a set of 1-1 contractive IFS maps, Wj:
X^X, \<i<N,
2.	A set of associated constants «y G M, 1 < i < A'^
3.	A set of associated place-dependent probability functions pr.X^ (0,°°), I <i<N,
4.	A "condensation multifunction", ß{x) G For each xeX, ß{x) e Jf is an interval in Y.
These ingredients, which comprise an A'-map "Iterated Function System on Multifiinctions" (IFSMF), are
now used to define the following fractal transform operator, T ■. ^ ,
{TF){x)=ß{x) + J^Pi{x)aiF{wT\x)). i=i
The reader will note that this operator is a multifianction analog of the "normal" IFS with condensation in Eq. 15.
Theorem 6 (La Torre et al, 2009a) The following inequalities hold:
( N
dp{T{F), T{G)) < N^P-'yP J^afsfpP d^iF, G),
d^{T{F)J{G))<
N
sup ^ aiPi{x)

where pj =	and Sj > 0 are such that
dix{wi{x)) < Sjdß{x).
Example: In Fig. 5, approximations to the attractor multifimctions onX= [0,1] are plotted for two IFSMF with two contractive IFS maps Wj. The top image corresponds to the attractor for the following IFSMF:
wi(x) = 0.6x, W2{x) = 0.6x+0A, ß{x) = [0.5,1.0].
«1 = 0.7, «2 = 0.5,
piix)=0.5, pi{x)=0.5,
(19)
The bottom image corresponds to the attractor of the IFSMF with the same Wj maps and Uj and pi constants as above but with the following jS-function:
iS« = [0,1],
ß{x) = [0.5,1.5]
0<x<0.5, 0.5 <^-<1. (20)
Note that the sets w\{X) and W2{X) overlap over the interval [0.4,0.6],
The multifunction attractor	corresponding
to Eq. 19 exhibits tiny jumps at x = 0.4 and x =
0.6,	the endpoints of the region of this region of overlap. Because of the self-similarity of the IFSMF, these jumps will be propagated throughout the multifunction. However, the jumps are quite small because the condensation multifunction j8(x) is the same over X= [0,1]. On the other hand, the condensation multifunction ß{x) demonstrates a significant change at point 0.5. This, along with the jumps associated with the overlapping Wi maps, produces much more irregular upper and lower bounds of the intervals comprising .
An inverse problem for multifunction approximation in the space can be formulated as follows: Given a multifunction F G find a contractive A'-map IFSMF operator T ■. ^ that admits a unique fixed point F G ß such that cL {F, F) is sufficiently small. Once again, we consider the simplification of this problem provided by the Collage Theorem. The inverse problem then becomes one of finding a contractive IFSMF operator T that maps the "target" multifunction F as close to itself as possible,
1.e.,	the collage distance do(F, TF) is made as small as possible. The following inequalities are useful for this approach.
Multifunction attractor F{x) of IFSMF in Eq. 19.
Multifunction attractor F{x) of IFSMF in Eq. 20.
Fig. 5. Pictorial representation of the multifunction attractors F{x) for the two IFSMF in the above Example. In each case, for an x G X = [0,1], F{x) is an interval [a{x), . The lower and upper bounds of these intervals, a(x) and b{x), respectively, are plotted in the ßgures.
Theorem 7 (La Torre et ah, 2009a) The following inequalities hold:
dp{F, TF)P < [[mmF-mmTF[[P
+ ||maxF —maxTFII^
N
d.{F, TF)<J^ A-sup max{4.(x),4-W} , where
Ai{x) = [mmF{x)-mm{ß{x) + aiF{wT^){x))[,
Ai{x) = I maxF{x) - max(jS{x) + wr\x)))[ and
Pi = supPi{wi{x)) .
xeX
Iterated Multifunction Systems: In this section, we describe a multifunction extension of IFS. In what follows, X will once again denote a base space, typically [0,1]". We now consider a set of n multiflinctions Tj : X ^	i G \ ...n (for each
i, TiX G	for all x X). Now construct the
multifunction T: X^ X, where
Tx=\jTjx, \/xeX. i=i
(21)
Assume that the multifunctions 7} are contractions with contractivity factors q G [0,1), that is,
TiX, Tiy) < Cid{x,y), V X (22)
Then there exists an element xG X such that xG Tx (Kunze et al, 2007). The element x is known as a fixed point of T. Note that xis not necessarily unique.
Now given a compact set A e Jf{X) consider the image
T{A)=[jTae^{X).	(23)
aeA
Since T : {X,d) {Jf{X),d},) is a continuous Sanction then T{A) is a compact subset of Jf{X). Therefore, we can construct a multifianction T* :
Jf{X) ^ Jf{X) as follows:
1.	For each	define F (/1) = T{A).
2.	Consider the Hausdorff distance on Jf{X): Given two subsets A,Be Jf{X), define
dhh{A, B) = max{sup inf	sup Mdh{x,y)}.
(24)
It then follows that T* : .^{X) ^ .^{X) and
dtt{r{A),r{B))<cdt{A,B). (25)
Now given a point x ^ X and a compact set 4 C A' we know that the Sanction d{x,a) has at least one minimum point a when a G A. We call a the projection of the point x on the set A and denote it as a = KxA. Obviously a is not unique but we choose one of the minima. We now define the following projection fianction ^associated with a multifunction T defined as P{x) = nx{ Tx). We therefore have the following result, proved in Kunze et al. (2008).
Theorem 8 Let {X, d) be a complete metric space and Ti'.X^ Jif{X) be a ßnite number of contractions with contractivity factors CjG [0,1). Let c = maxjCj. Then
1. For all compact A C X there exists a compact subset A (Z X such that = P{A„) A when n +00,
Z ÄcUiTiiÄ).
As in the previous sections, let ^{X) be the set of probability measures on	and consider
the complete metric space {J^{X),dM)- Given a set of multifianctions T : X ^ X with associated probabilities pi, one can now consider generalized Markov operators on ^{X).
Fractal transforms on measure-valued functions: In
what follows, X = [0,1] will denote the "base space," i.e., the support of the images. R^ C M will denote a compact "greyscale space" of values that our images can assume at any xG Xand B will denote the Borel o algebra on R^ with the Lebesgue measure. Let ^ denote the set of all Borel probability measures on R^ and du the Monge-Kantorovich metric on this set. For a given M> 0, let C ^ be a complete subspace of ^ such that dniß, v) < M for all ju, v G . We now define
Y = {fx: X ^	measurable} (26)
and consider on this space the following metric
C/HM,V)= /	(27)
J X
We observe that dy is well defined, since ß and v are measurable fianctions, c/// is bounded and so the fianction ^ (x) = c///(/x(x), v(x)) is integrable on X.
Theorem 9 (La Torre et al, 2009b) The space ( Y, dy) is complete.
We now construct and analyze a fractal transform operator M on the space {Y^dy) of measure-valued fianctions. We list the ingredients for a fractal transform operator in the space Y. The reader will note that they form a kind of blending of IFS-based methods on measures (IFSP) and fianctions (IFSM). For simplicity, we assume that A'= [0,1]. The extension to [0,1]" is straightforward.
1. A set of A'^ one-to-one contraction affine maps Wi:
X ^ X, wj{x) = SiX+ ai, with the condition that
yjl,wix) = x.
2. A set of A'^ greyscale maps :

assumed
to be Lipschitz, i.e., for each i, there exists a «y > 0 such that
I^Ih) - MWI < aj\ ti-t2l Mti, /2 G R^, (28)
3. For each x G A', a set of probabilities Pi{x), i = 1, • • • , A'^ with the following properties:
-	are measurable as functions of x
-	pi{x) = 0 if X ^ Wi{X) and
-	= 1 forallA-GX
The action of the fractal transform operator M: Y ^Y defined by the above is as follows: For aß eV and any subset 5 c M^,
n
v{x){S) = {Mß{x)){S) = J^p,{x)ß{w7\xm7\S)).
(29)
i=i
Theorem 10 (La Torre et al, 2009b) Let pi = Then forßi,ß2 G Y,
/ n	\
dY{Mßi,Mß2)< J,\si\aiPi dyißußi). (30)
/
Corollary 4 Let pj = sup^^;^ pi{x). Then M is a contraction on (7, dy) if
i=i
(31)
Consequently there exists a measure-valued mapping ß eY, such that ß = Mß.
Examples:
1. The fractal transform M defined by the following two-IFS-map system onX= [0,1]:
Wl{x) = ^x, W2{x) =
hit) =
The sets w\{X) and W2{X) overlap at the single point X = ^ so we let
Pl{x) = \, P2ix)=0 Pl{x)=0, p2{x) = \
Pi (^) =/52(2) = 2

It is easy to confirm that Mis contractive. Its fixed point ß is given by
ß{x) = 8{t-x), xG[0,l], (32)
where 5 denotes the Dirac delta fianction.
2. A "perturbation" of the above fractal transform M that is produced by adding the following IPS and associated greyscale maps:
The sets wi (X) and W3 (X) overlap over the entire subinterval [0, so we let
pi{x) = pi{x) = -,
p2{x) = 0
Pi{x) = pi{x) = Q, p2{x) = \

Once again, it is easy to confirm that M is contractive. Its fixed point ß{x) is sketched in Fig. 6.
At this point, we mention that it is difficult to produce a sketch of the fixed point that will capture all of its detailed structure. First of all, the plot of ß{x) in the figure has the appearance of a (sheared) Sierpinski gasket with vertices at (0,0), (0,1/3) and (1,1). The "gaps" in this gasket reflect regions of low measure. Any attempt to increase the darkness of these regions would remove any idea of the self-similar variations in ß{x) in the x-direction.
The important feature to note in this figure is that the the overlapping of the wi and W3 maps over [0, is responsible for the self-similar "splitting" of the measures ß(x) (hence lighter shading) over this interval, since ^ produces an upward shift in the greyscale direction. Since W2{}^ maps the support [0,1] of the entire measure-valued fianction onto [5,1], the self-similarity of the measure over [0, is carried over to [5,1].
Fig. 6. A sketch of the invariant measure ß(x) for the three-IFSmap fractal transform in Example 2, xgX= [0,1],7GR^=[0,1].
We now show that the moments of measures in the space (V, dy) also satisfy recursion relations when the grey scale maps are affine. We now consider the local or x-dependent moments of a measure G Y, defined as follows,
gn{x)= f s^dß^is), m = 0,\,2,----
(33)
where we use the notation ßj^ = ß{x) in the Lebesgue integral for simplicity. By definition, ^o(^) = 1 for xG X. Obviously the Sanctions g^ are measurable on X (since ß{x) are measurable) and bounded so that gm G L^ {X, J^). We now derive the relations between the moments of a measure ß G Y and the moments of V = Mß where M is the fractal transform operator defined in Eq. 29.
Let h„ denote the moments of v = Mß. Then
hnix) = f ^diMßUs)
jRg


For affine greyscale maps of the form (p{s) = «^5+ßj, we have
I- N
hn{x)= /	+
1
J=0
i=l
where
cnj —
The reader may compare the above result to that of Eq. 8 for the IFSP case. The place-dependent moments h„{x) are related to the moments gn evaluated at the preimages wj^{x). And it is the greyscale (^(s) maps that now "mix" the measures, as opposed to the spatial IPS maps wj{x) in Eq. 8.
In the special case that ß = ß = Mß, the fixed point of M, then hn{x) = gn{x) and we have
gn{x) = X J=0
N
i=i
gji^Y'i^)) ■
In other words, the moments gn{x) satisfy recursion relations that involve moments of all orders up to n evaluated at preimages wj^{x). Note that this does not yield a rearrangement, analogous to Eq. 11, which will permit a simple recursive computation of the moments gn{x). Nevertheless, the moment functions can be computed recursively (see (La Torre etal, 2009b)).
MEASURE-VALUED FUNCTIONS, NONLOCAL IMAGE PROCESSING AND FRACTAL CODING
Nonlocal image processing, the manipulation of the value of an image Sanction u{x) based upon values of »(/jt) elsewhere in the image, has recently received a great deal of attention, due in part to the success of the nonlocal means image denoising method (Buades etal., 2010). (This is in contrast to standard image processing methods which are local in nature, i.e., the points /jt lie in a neighbourhood of x.) Fractal image coding, outlined earlier in this paper, is another example of a nonlocal image processing method. In fact, both of these methods may be viewed under the umbrella of a more general model of affine image self-similarity (Alexander et al, 2008), in which subblocks of an image are approximated by other sublocks of the image.
In La Torre et al. (2009b), we showed how the multifianction/measure-valued representation of images, outlined in the previous section of this paper, may be usefial in nonlocal image processing methods. In these methods, the value of an image function at a point x G X is, replaced by a transformed value Tu{x) which is usually composed by several values of the image function u{yi^) that lie elsewhere in the image. It may be useful to store these values in a measure or distribution ß{x) before performing the final projection of these values in order to produce the transformed value Tu{}dj. For example, the measure ß{x) could be used to characterize the local self-similarity of an image at a point xG X. The measure-valued formalism was used to analyze both the methods of nonlocal means denoising as well as fractal image coding. We now outline briefly its application to the latter.
Historically, most fractal image coding research focussed on its compression capabilities - obtaining acceptable accuracy with the smallest possible domain pool in order to minimize (i) search times and (ii) storage of the fractal code. The fact that range blocks Ri of an image are, in general, well approximated by a good number of domain blocks Dj does not seem to have been emphasized or exploited. Consequently, investigations generally focussed on the results yielded by optimal domain blocks of the pool and not on the possible use of suboptimal ones. The reader will recall that the fractal coding method described earlier in this paper was based on the selection of the best domain block for each range block.
More recently, however, the redundancy of good domain/range pairings has been exploited (Alexander, 2005) in order to perform image denoising. As in the case of nonlocal means denoising, the use of several domain blocks corresponds to an averaging over multiple samples, resulting in a reduction of noise variance. This may be viewed as a multiparent fmctal block coding method.
At this point, it is important to mention that the idea of using several domain blocks for each range block is not new. Some examples of multiparent fractal coding schemes may be found in Gharavi-Alkhansari and Huang (1994); Vines (1995), and Furusawa and Nakagawa (2004).
A simple measure-valued method associated with fractal coding: Here we outline a simple multiparent block fractal coding scheme that results from a modification of the block-based fractal coding method outlined an earlier section. This multiparent scheme lends itself nicely to a measure-based formalism.
For convenience, we consider the same (square) range and domain block arrangement used in the fractal image coding scheme outlined earlier. For each range block Rj, we compute the Ajj approximation errors associated with a//domain blocks Dj, c/ Eq. 16. (Recall that for each range/domain pairing {Rj,Dj) there are eight spatial contraction/decimation maps w^j : Dj ^ Rj, \ < k < Once again, for simplicity of notation, we shall omit the k index.) The optimal greyscale map minimizing the error Aij will be denoted as
^ij{t) = aijt + ßij.	(34)
For this pairing we also assign a weight pij, normalized so that
No J=i
(35)
For each range block it would seem natural to employ higher weights pij for those domain blocks Dj that yield lower collage errors Aij. Here we consider the following weighting scheme.
1
= ^exp
(
~hP
(36)
where P > 0, A > 0 and 4 = Ejexp(-Aj.//?0 is the normalization factor guaranteeing that Xj Pij = 1 for each i. In practice, P is either 1 or 2. Here, we shall employ P = 2, i.e., a Gaussian-type weighting. Regarding the adjustable parameter h >0:
1. In the limit A ^ 0+, the pij with the smallest error Aij will be selected.
2. In the limit h ^ <=°, all pjj become equal, i.e., all
range/domain pairings are employed equally.
The resulting multiparent block transform operator T is then defined as
Nd
v{x) = {Tu){x) = £ Pijaiju{w7l{x))+ßij, J=i
x(^Ri, l<i<NR. (37)
This definition represents a generalization of the fractal transform operator of Eq. 17 since not only one but several, perhaps all, domain blocks Dj G & can contribute to the modification of u{x) for x G Under appropriate restrictions on the a parameters, this transform is attractive, which implies the existence of a fixed point Sanction u{x) which will provide an approximation to the original image being fractally coded.
At this point, we emphasize that the above multiparent fractal transform operator T takes all of the preimages u{wTj^{x)) of an image function value u{x) and from them produces a single value v{x). We now illustrate how the measure-valued image fianction can be used to examine the range of values assumed by these preimages.
First, we associate with the image fianction a corresponding measure-valued image fianction }x{x) G (F, dy) as follows:
=	xeX.
(38)
Here, dt denotes a unit point-mass measure at t e M^. With an eye to Eq. 29, we now define a measure-valued image function v = Mß G 7 as follows: For any measurable set 5 C M^ = [0,1] and any x G Ri, we define
viS) = iMß)ix)iS) j=i
(39)
Given a range block Rj, then at each point/pixel x G Ri, we keep track of all greyscale values of the image fianction u that are mapped to x by a domain/range mapping Wij and modified by the corresponding greyscale map (pjj. These values are weighted and combined to define the probability measure v{x).
This idea is illustrated below for the Boat image shown earlier First of all, we shall concentrate on the row of pixels u{256,J),J= I,-- - ,256 in the image. These are the pixels that run from the midpoint of the left edge to the center of the Boat image. These greyscale values are plotted in Fig. 7 below.
200	250
Fig. 7. Greyscale values of the (normalized) Boat image u(256, j), j = 1, 256.
The "dips" in the above plot corresond to the various masts and the two prominent areas of increased greyscale value/brightness correspond to the lighthouse (left) and the boat's cabin (right). The dark area at the extreme right represents the shaded part of the cabin.
In Fig. 8 we show pictorial representations of the measure-valued functions v(x) for these pixels corresponding to three values of the parameter h, namely h = 0.01, h = 0.1 and h = 1.0.
In these figures, darker regions have higher associated measures. In the leftmost figure, the very small parameter value h = 0.01 concentrates the measures close to the Boa^t image values u(256, j) since only domain blocks with low approximation errors are used to approximate them. As h is increased to 0.1, the measure (middle figure) becomes more diffuse, as blocks with higher errors are admitted. The measure associated with h = 1.0 (right figure) is virtually identical to that of h = 0.1.
In each of the above three cases, however, the measure v(x) at a pixel x represents a "preprocessing" of the fractal coding method, essentially giving a picture of the preimages of the pixel value u(x) that are then used to construct a transformed value v(x).
Even more interesting is the effect of (additive) noise on these measures. As expected, the measures become even more diffuse. This feature, along with a simple associated denoising method, was analyzed in La Torre et al. (2009b).
CONCLUSIONS
Starting with the classical definitions of gener^alized fractal tr^ansforms (GFT), we have reviewed the results of more recent work on the
formulation of GFTs over spaces of multifunctions, including the space of measure-valued functions. These new formalisms may be useful in nonlocal image processing. We plan to explore the further use of these methods in characterizing image self-similarity in future papers.
ACKNOWLEDGEMENTS
This work has been supported in part by a Discovery Grant (ERV) from the Natural Sciences and Engineering Research Council of Canada (NSERC).
REFERENCES
Alexander SK (2005). Multiscale methods in image modelling and image processing. PhD Thesis. Department of Applied Mathematics, University of Waterloo.
Alexander SK, Vrscay ER, Tsurumi S (2008). A simple, general model for the affine self-similarity of images. In: Image Analysis and Recognition, ICIAR 2008. Lect Notes Comput Sci 5112:192-203.
Barnsley MF. Fractals Everywhere. New York: Academic Press.
Barnsley MF, Demko S (1985). Iterated function systems and the global construction of fractals. Proc Roy Soc London Ser A 399:243-75.
Barnsley MF, Ervin V, Hardin D, Lancaster J (1985). Solution of an inverse problem for fractals and other sets. Proc Natl Acad Sci USA 83:1975-7. Barnsley MF, Hurd L (1993). Fractal Image Compression.
Wellesley, Massachussetts: AK Peters. Buades A, Coll B, Morel JM (2010). A review of image denoising algorithms, with a new one. SIAM Rev 52:113-47.
Cabrelli CA, Forte B, Molter UM, Vrscay ER (1992). Iterated fuzzy set systems: a new approach to the inverse problem for fractals and other sets. J Math Anal Appl 171:79-100.
Centore P, Vrscay ER (1994). Continuity of fixed points for attractors and invariant meaures for iterated function systems. Can Math Bull 37:315-29.
Fisher Y (1995). Fractal Image Compression, Theory and
Application. New York: Springer Verlag. Forte B, Mendivil F, Vrscay ER (1999). IFS operators on integral transforms. In: Fractals: Theory and Applications in Engineering. London: Springer Verlag.
Forte B, Vrscay ER (1995). Solving the inverse problem for measures using iterated function systems: A new approach. Adv Appl Prob 27:800-20.
Forte B, Vrscay ER (1998a). Theory of generalized fractal transforms. In: Fisher Y, ed. Fractal Image Encoding
0
h = OM
h=0.\
h= 1.0
Fig. 8. Pictorial representations of the measure v (x) for the u(25 6 ,J),J = 1, • • • 25 6, row of pixels for the Boat image. Three values of the weighting parameter h in Eg. 36 have been used, with P =2. Darker regions have higher measures.
and Analysis, NATO ASl Series F Vol 159. New York: Springer Verlag.
Forte B, Vrscay ER (1998b). Inverse problem methods for generalized fractal transforms. In: Fractal Image Encoding and Analysis, ibid..
Furusawa R, Nakagawa M (2004). Fractal image coding with multiscaling-domain. Electron Comm Jpn 3 87(2):79-87.
Gharavi-Alkhansari M, Huang TS (1994). Generalized image coding using fractal-based methods. Proc Intern Picture Coding Symp PCS'94. Sacramento, California.
Hutchinson J (1981). Fractals and self-similarity. Indiana Univ J Math 30:713-47.
Jacquin A (1989). A Fractal Theoiy of Iterated Markov Operators with Applications to Digital Image Coding. PhD Thesis. School of Mathematics, Georgia Institute of Technology.
Jacquin A (1992). Image coding based on a theoiy of iterated contractive image transformations. IEEE Trans Image Proc 1(1): 18-30.
Kunze H, La Torre D, Vrscay ER (2007). Contractive multifunctions, fixed point inclusions and iterated multifunction systems. J Math Anal Appl 330:157-73.
Kunze H, La Torre D, Vrscay ER (2008). From Iterated Function Systems to Iterated Multifunction Systems. Comm Appl Nonlinear Anal 15(4): 1-15.
Kunze H, Vrscay ER (1999). Solving inverse problems for ordinary differential equations using the Picard contraction mapping. Inverse Probl 15(9):2669-75.
La Torre D, Mendivil F (2008). Iterated function systems on multifunctions and inverse problems. J Math Anal Appl 340(2): 1469-79.
La Torre D, Mendivil F (2009). Union-additive multimeasures and self-similarity. Comm Math
Anal 7(2): 51-61.
La Torre D, Mendivil F, Vrscay ER (2009a). Iterated function systems on multifunctions. In: Aletti G, ed. Math Eveiywhere - Deterministic and Stochastic Modelling in Biomedicine, Economics and Industiy. Heidelberg, Berlin: Springer-Verlag.
La Torre D, Vrscay ER, Ebrahimi M, Bamsley M (2009b). Measure-valued images, associated fractal transforms and the affine self-similarity of images. SIAM J Imag Sei 2(2):470-507.
Lu N (2003). Fractal imaging. New York: Academic Press.
Mandelbrot B (1977). The Fractal Geometry of Nature. New York: WH Freeman.
Mendivil F, Vrscay ER (1997). Correspondence between fractal-wavelet transforms and iterated function systems with grey-level maps. In: Levy-Vehel J, Lutton E, Tricot C, eds. Fractals in Engineering: From Theoiy to Industrial Apphcations. London: Springer Verlag.
Mendivil F, Vrscay ER (2002). Fractal vector measures and vector calculus on planar fractal domains. Chaos Soliton Fract 14:1239-54.
Vines G (1995). Orthogonal basis IFS. In: Fisher Y ed. Fractal Image Compression op cit
Vrscay ER (1990). Moment and collage methods for the inverse problem of fractal constmction with iterated function systems. In: Peitgen HO, Henriques JM, Penedo LF, eds. Fractals in the Fundamental and Applied Sciences. 443-61.
Vrscay ER (1998). A generahzed class of fractal-wavelet transforms for image representation and compression. Can J Elec Comp Eng 23:69-84.
Vrscay ER, Roehrig CJ (1989). Iterated function systems and the inverse problem of fractal construction using moments. In: Kaltofen E, Watt SM, eds. Computers and Mathematics New York: Springer Verlag.