Iterated Function Systems

[Ross Hilbert]

David,

Yep, the "log(Pi)/log(Si)=constant" method provides more accurate results than the method from "Chaos and Fractals"! One small problem, if the set of affine transformations includes a transformation with scale=1 (e.g., a rotation), the algorithm fails since the log(Si)=0. In that case, do you simply use a method like that from "Chaos and Fractals" or is there a better solution?

Ross

[cKleinhuis]

you can now use tex for formatting formulas:

[David Makin]

You are correct, the method only works for scales >0 and <1.

If you have scales of "1" for a "whole" transform included as part of the IFS then I don't really know how you're getting a convergent fractal

Having said that you can either use the Barnsley/"Chaos and Fractals" method instead or maybe use say 0.99 as the scale instead of 1 or (and I think I did this) first normalise the scales so the sum of them is 1 then apply the log(Pi)/log(Si)=constant method using the normalised scales as Si.

One thing to mention is that in

log(Pi)/log(Si)=constant

the constant is the fractal dimension if you use the linear scales.
If using the scales directly from the determinants then the constant is dimension/2 for 2D fractals and dimension/3 for 3D fractals.

[Ross Hilbert]

Thanks David,

As to your question:
  If you have scales of "1" for a "whole" transform included as part of the IFS then I don't really know how you're getting a convergent fractal...

See http://classes.yale.edu/Fractals/IntroToFrac/InvProb/SnowIFS/Snow2.html for an example. They state: "The Random Algorithm still will generate the snowflake, so long as on average the transformations applied are contractions."

Just using 0.99 gives results that are far worse than using the Barnsley/"Chaos and Fractals" method. I tried your suggestion of normalizing the scales first and that works fine.

Thanks for your clarification with respect to the relationship of the constant (in log(Pi)/log(Si)=constant) and the fractal dimension, but just to be clear, this is not relevent with respect to the algorithm to calculate the probabilities from the scales, right?

Ross

[Ross Hilbert]

For completeness, here is the updated code using the reference David suggested.
Ross

Code:

'
' Define Affine object.
'
Object Affine { A : B : C : D : E : F }
'
' Normalize the count values in weights[].
' On return, the values will all be between 0 and 1
' and the sum of all values will equal 1.
'
void Math.NormalizeWeights(weights[], count) {
  sum = 0
  
  for (i = 0, i < count, i += 1) {
    sum += weights[i]
  }
  for (i = 0, i < count, i += 1) {
    weights[i] /= sum
  }
}
'
' Return the determinate of an affine transformation matrix.
'
Complex Affine.Determinate(Affine a) = a.A*a.D - a.B*a.C
'
' Affine.TransformPoint applies an affine transformation 
' matrix to a point z and returns the resulting point.
'
' A B E     z.X     A*z.X + B*z.Y + E
' C D F  *  z.Y  =  C*z.X + D*z.Y + F
' 0 0 1     1       1
'
Complex Affine.TransformPoint(Affine a, z) {
  return Complex( \
    a.A*z.x + a.B*z.y + a.E, \
    a.C*z.x + a.D*z.y + a.F  \
  )
}
'
' Generate a random contraction map.
'
' To be a contraction map, the affine transformation
' coefficients must satisfy the following conditions:
'
'   a^2 + c^2 < 1
'   b^2 + d^2 < 1
'   a^2 + b^2 + c^2 + d^2 < 1 + (a*d - c*b)^2
'
Affine Affine.RandomContractionMap(center, radius) {
  while (True) {
    a = Random.NumberInRange(-1, 1)
    w = Sqrt(1 - a^2)
    c = Random.NumberInRange(-w, w)
    
    b = Random.NumberInRange(-1, 1)
    w = Sqrt(1 - b^2)
    d = Random.NumberInRange(-w, w)
  
    if (a^2 + b^2 + c^2 + d^2 < 1 + (a*d - c*b)^2) {
      break
    }
  }
  e = center.x + Random.NumberInRange(-radius, radius)
  f = center.y + Random.NumberInRange(-radius, radius)
  
  return Affine(a, b, c, d, e, f)
}
'
' Initialize weights w[] using the weighted average of
' the absolute value of the determinates of the affine 
' transformations in s[].
' 
' For a discussion of the following algorithm, see:
'   "Chaos and Fractals" 
' by Peitgen, Jurgens, Saupe, Section 6.3 pages 327-329.
'
void Affine.InitializeWeights(Affine s[], w[], count) {
  for (i = 0, i < count, i += 1) {
    w[i] = Max(0.01, Abs(Affine.Determinate(s[i])))
  }
  Math.NormalizeWeights(w[], count)
}
'
' Affine.InitializeWeights2F is used by Affine.InitializeWeights2.
'
Complex Affine.InitializeWeights2F(p, w[], count, imax, wmaxLog) {
  f = p - 1
  
  for (i = 0, i < count, i += 1) {
    if (i <> imax) {
      f += p^(Log(w[i])/wmaxLog)
    }
  }
  return f
}
'
' Affine.InitializeWeights2 initializes the weights
' associated with the count affine transformations
' given in s[]. 
'
' For a discussion of the following algorithm, see:
'   "A Multifractal Analysis of IFSP Invariant Measures 
'    With Application to Fractal Image Generation"
' by J.M. Gutierrez, A. Iglesias, and M.A. Rodriguez
' http://grupos.unican.es/ai/meteo/articulos/fractals96.pdf
' 
void Affine.InitializeWeights2(Affine s[], w[], count) {
  wmax = 0
  imax = 0
  '
  ' Initialize the scale factors in w[].
  '
  for (i = 0, i < count, i += 1) {
    t = Max(0.01, Abs(Affine.Determinate(s[i])))
    
    if (wmax < t) {
      wmax = t
      imax = i
    }
    w[i] = t
  }
  '
  ' Check if any scale factors are >= 1 and if so, renormalize w[].
  '
  if (wmax >= 1) {
    Math.NormalizeWeights(w[], count)
    wmax = w[imax]
  }
  wmaxLog = Log(wmax)
  '
  ' Compute the probability p associated with s[imax].
  '
  #include Math.NewtonMethod2( \
    "p", \
    "Affine.InitializeWeights2F(p, w[], count, imax, wmaxLog)", \
    "0.5" \
  )
  w[imax] = p
  '
  ' Set w[i] to the probability associated with s[i].
  '
  for (i = 0, i < count, i += 1) {
    if (i <> imax) {
      w[i] = p^(Log(w[i])/wmaxLog)
    }
  }
}
'
''''''''''''''''''''''''''''''''''''''''''''''''''''
' N is the number of Affine transformations.
' s[N] are the set of N random Affine transformations.
' w[N] are the weights associated with the transformations.
''''''''''''''''''''''''''''''''''''''''''''''''''''
'
Complex w[N]
Affine s[N]

for (i = 0, i < N, i += 1) {
  s[i] = Affine.RandomContractionMap(0, 1)
}
Affine.InitializeWeights2(s[], w[], N)

[David Makin]

"See http://classes.yale.edu/Fractals/IntroToFrac/InvProb/SnowIFS/Snow2.html for an example. They state: "The Random Algorithm still will generate the snowflake, so long as on average the transformations applied are contractions.""

Yea - I remember now, the chaos game will still work, as will the deterministic method if you include the "artificial" termination of the IFS tree at a given maximum depth in addition to terminating at a given overall scale

With respect to the relevance of the fractal dimension, it's relevant in that it's the constant involved and in the fact that you can use log(Pi)/log(Si)=constant to generate a fractal based on the original so that it has a particular dimension but the same "form" as the original.

e.g. I use:

  if @fixdim
    gs = fd/@dim
    gi = ft
    repeat
      gt = tfs[gi]
      tfs[gi] = tfs[gi]^gs
      gt = tfs[gi]/gt
      tf0[gi] = tf0[gi]*gt
      tf1[gi] = tf1[gi]*gt
      tf2[gi] = tf2[gi]*gt
      tf3[gi] = tf3[gi]*gt
    until (gi=gi-1)<0
    fd = @dim
  endif

Where:
 fd is the original calculated dimension
 @dim is the requested dimension
 ft is "first transform" number (usually number of transforms-1)
 gi is the loop transform number
 tfs[] are the scales
 tf0[] to tf3[] are the 2*2 transformation matrix values.

This method is based on keeping the same probabilities.

[Ross Hilbert]

Cool, I learned a lot. Thanks.

[David Makin]

Just to add one other thing with respect to generating IFS fractals using the chaos game method - you probably know anyway - the statistical randomness of the pseudo-random you use is almost as important as getting the probabilities  "correct" - the default C++/Windows rand() is *terrible*.
I did quite a lot of searching on this and found that the best pseudo-random generator is one called the "Mersenne Twister", I've lost my links but you should be able to find code for it.

[David Makin]

e.g. I use:

  if @fixdim
    gs = fd/@dim
    gi = ft
    repeat
      gt = tfs[gi]
      tfs[gi] = tfs[gi]^gs
      gt = tfs[gi]/gt
      tf0[gi] = tf0[gi]*gt
      tf1[gi] = tf1[gi]*gt
      tf2[gi] = tf2[gi]*gt
      tf3[gi] = tf3[gi]*gt
    until (gi=gi-1)<0
    fd = @dim
  endif

Where:
 fd is the original calculated dimension
 @dim is the requested dimension
 ft is "first transform" number (usually number of transforms-1)
 gi is the loop transform number
 tfs[] are the scales
 tf0[] to tf3[] are the 2*2 transformation matrix values.

This method is based on keeping the same probabilities.

Sorry - important to mention that in the above the scales used are the linear scales i.e. sqrt(abs(determinant)). (Use cuberoot(abs(determinant)) for 3D).

[David Makin]

Note that using the chaos game method to render IFS then the fractal dimension is important for when making your program as "black box" as possible.
I mean that if you render a fractal at one size with a given iteration count for the chaos game method and then wish to render it at say a much larger size then the required iteration count for "the same" image is based on the fractal dimension and ideally your program should be able to automatically adjust the actual iteration count used accordingly.

Normally you'd use something like:

count = (userparam*sqrt(pixelwidth*pixelheight))^fractaldimension

[Ross Hilbert]

Great information David, thanks again.
Ross

[David Makin]

Note that using the chaos game method to render IFS then the fractal dimension is important for when making your program as "black box" as possible.
I mean that if you render a fractal at one size with a given iteration count for the chaos game method and then wish to render it at say a much larger size then the required iteration count for "the same" image is based on the fractal dimension and ideally your program should be able to automatically adjust the actual iteration count used accordingly.

Normally you'd use something like:

count = (userparam*sqrt(pixelwidth*pixelheight))^fractaldimension

Just to add that the above works best when the whole fractal is on-screen, if not then the accuracy will depend on how the density of the fractal varies across it's area and which area of the fractal is in view.