'
' 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)