Adding upstream version 2.10.36.upstream/2.10.36

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

1 files changed, 1092 insertions, 0 deletions

diff --git a/plug-ins/ifs-compose/ifs-compose-utils.c b/plug-ins/ifs-compose/ifs-compose-utils.c
new file mode 100644
index 0000000..32c05cc
--- /dev/null
+++ b/plug-ins/ifs-compose/ifs-compose-utils.c
@@ -0,0 +1,1092 @@
+/* GIMP - The GNU Image Manipulation Program
+ * Copyright (C) 1995 Spencer Kimball and Peter Mattis
+ *
+ * IfsCompose is a interface for creating IFS fractals by
+ * direct manipulation.
+ * Copyright (C) 1997 Owen Taylor
+ *
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#include "config.h"
+
+#include <stdlib.h>
+#include <string.h>
+
+#include <gdk/gdk.h>
+
+#include <libgimp/gimp.h>
+
+#include "ifs-compose.h"
+
+
+typedef struct
+{
+  GdkPoint point;
+  gdouble angle;
+} SortPoint;
+
+
+/* local functions */
+static void     aff_element_compute_click_boundary (AffElement     *elem,
+                                                    gint            num_elements,
+                                                    gdouble        *points_x,
+                                                    gdouble        *points_y);
+static guchar * create_brush                       (IfsComposeVals *ifsvals,
+                                                    gint           *brush_size,
+                                                    gdouble        *brush_offset);
+
+
+void
+aff2_translate (Aff2    *naff,
+                gdouble  x,
+                gdouble  y)
+{
+  naff->a11 = 1.0;
+  naff->a12 = 0;
+  naff->a21 = 0;
+  naff->a22 = 1.0;
+  naff->b1 = x;
+  naff->b2 = y;
+}
+
+void
+aff2_rotate (Aff2    *naff,
+             gdouble  theta)
+{
+  naff->a11 = cos(theta);
+  naff->a12 = sin(theta);
+  naff->a21 = -naff->a12;
+  naff->a22 = naff->a11;
+  naff->b1 = 0;
+  naff->b2 = 0;
+}
+
+void
+aff2_scale (Aff2    *naff,
+            gdouble  s,
+            gboolean flip)
+{
+  if (flip)
+    naff->a11 = -s;
+  else
+    naff->a11 = s;
+
+  naff->a12 = 0;
+  naff->a21 = 0;
+  naff->a22 = s;
+  naff->b1 = 0;
+  naff->b2 = 0;
+}
+
+/* Create a unitary transform with given x-y asymmetry and shear */
+void
+aff2_distort (Aff2    *naff,
+              gdouble  asym,
+              gdouble  shear)
+{
+  naff->a11 = asym;
+  naff->a22 = 1/asym;
+  naff->a12 = shear;
+  naff->a21 = 0;
+  naff->b1 = 0;
+  naff->b2 = 0;
+}
+
+/* Find a pure stretch in some direction that brings xo,yo to xn,yn */
+void
+aff2_compute_stretch (Aff2    *naff,
+                      gdouble  xo,
+                      gdouble  yo,
+                      gdouble  xn,
+                      gdouble  yn)
+{
+  gdouble denom = xo*xn + yo*yn;
+
+  if (denom == 0.0)     /* singular */
+    {
+      naff->a11 = 1.0;
+      naff->a12 = 0.0;
+      naff->a21 = 0.0;
+      naff->a22 = 1.0;
+    }
+  else
+    {
+      naff->a11 = (SQR(xn) + SQR(yo)) / denom;
+      naff->a22 = (SQR(xo) + SQR(yn)) / denom;
+      naff->a12 = naff->a21 = (xn * yn - xo * yo) / denom;
+    }
+
+  naff->b1 = 0.0;
+  naff->b2 = 0.0;
+}
+
+void
+aff2_compose (Aff2 *naff,
+              Aff2 *aff1,
+              Aff2 *aff2)
+{
+  naff->a11 = aff1->a11 * aff2->a11 + aff1->a12 * aff2->a21;
+  naff->a12 = aff1->a11 * aff2->a12 + aff1->a12 * aff2->a22;
+  naff->b1  = aff1->a11 * aff2->b1  + aff1->a12 * aff2->b2 + aff1->b1;
+  naff->a21 = aff1->a21 * aff2->a11 + aff1->a22 * aff2->a21;
+  naff->a22 = aff1->a21 * aff2->a12 + aff1->a22 * aff2->a22;
+  naff->b2  = aff1->a21 * aff2->b1  + aff1->a22 * aff2->b2 + aff1->b2;
+}
+
+/* Returns the identity matrix if the original matrix was singular */
+void
+aff2_invert (Aff2 *naff,
+             Aff2 *aff)
+{
+  gdouble det = aff->a11 * aff->a22 - aff->a12 * aff->a21;
+
+  if (det==0)
+    {
+      aff2_scale (naff, 1.0, 0);
+    }
+  else
+    {
+      naff->a11 = aff->a22 / det;
+      naff->a22 = aff->a11 / det;
+      naff->a21 = - aff->a21 / det;
+      naff->a12 = - aff->a12 / det;
+      naff->b1 = - naff->a11 * aff->b1 - naff->a12 * aff->b2;
+      naff->b2 = - naff->a21 * aff->b1 - naff->a22 * aff->b2;
+    }
+}
+
+void
+aff2_apply (Aff2    *aff,
+            gdouble  x,
+            gdouble  y,
+            gdouble *xf,
+            gdouble *yf)
+{
+  gdouble xt = aff->a11 * x + aff->a12 * y + aff->b1;
+  gdouble yt = aff->a21 * x + aff->a22 * y + aff->b2;
+
+  *xf = xt;
+  *yf = yt;
+}
+
+/* Find the fixed point of an affine transformation
+   (Will return garbage for pure translations) */
+
+void
+aff2_fixed_point (Aff2    *aff,
+                  gdouble *xf,
+                  gdouble *yf)
+{
+  Aff2 t1,t2;
+
+  t1.a11 = 1 - aff->a11;
+  t1.a22 = 1 - aff->a22;
+  t1.a12 = -aff->a12;
+  t1.a21 = -aff->a21;
+  t1.b1  = 0;
+  t1.b2  = 0;
+
+  aff2_invert (&t2, &t1);
+  aff2_apply (&t2, aff->b1, aff->b2, xf, yf);
+}
+
+void
+aff3_apply (Aff3    *t,
+            gdouble  x,
+            gdouble  y,
+            gdouble  z,
+            gdouble *xf,
+            gdouble *yf,
+            gdouble *zf)
+{
+  gdouble xt = (t->vals[0][0] * x +
+                t->vals[0][1] * y +
+                t->vals[0][2] * z + t->vals[0][3]);
+  gdouble yt = (t->vals[1][0] * x +
+                t->vals[1][1] * y +
+                t->vals[1][2] * z + t->vals[1][3]);
+  gdouble zt = (t->vals[2][0] * x +
+                t->vals[2][1] * y +
+                t->vals[2][2] * z + t->vals[2][3]);
+
+  *xf = xt;
+  *yf = yt;
+  *zf = zt;
+}
+
+static int
+ipolygon_sort_func (const void *a,
+                    const void *b)
+{
+  if (((SortPoint *)a)->angle < ((SortPoint *)b)->angle)
+    return -1;
+  else if (((SortPoint *)a)->angle > ((SortPoint *)b)->angle)
+    return 1;
+  else
+    return 0;
+}
+
+/* Return a newly-allocated polygon which is the convex hull
+   of the given polygon.
+
+   Uses the Graham scan. see
+   http://www.cs.curtin.edu.au/units/cg201/notes/node77.html
+
+   for a description
+*/
+
+IPolygon *
+ipolygon_convex_hull (IPolygon *poly)
+{
+  gint       num_new     = poly->npoints;
+  GdkPoint  *new_points  = g_new (GdkPoint, num_new);
+  SortPoint *sort_points = g_new (SortPoint, num_new);
+  IPolygon  *new_poly    = g_new (IPolygon, 1);
+
+  gint       i, j;
+  gint       x1, x2, y1, y2;
+  gint       lowest;
+  GdkPoint   lowest_pt;
+
+  new_poly->points = new_points;
+  if (num_new <= 3)
+    {
+      memcpy (new_points, poly->points, num_new * sizeof (GdkPoint));
+      new_poly->npoints = num_new;
+      g_free (sort_points);
+      return new_poly;
+    }
+
+  /* scan for the lowest point */
+  lowest_pt = poly->points[0];
+  lowest = 0;
+
+  for (i = 1; i < num_new; i++)
+    if (poly->points[i].y < lowest_pt.y)
+      {
+        lowest_pt = poly->points[i];
+        lowest = i;
+      }
+
+  /* sort by angle from lowest point */
+
+  for (i = 0, j = 0; i < num_new; i++, j++)
+    {
+      if (i==lowest)
+        j--;
+      else
+        {
+          gdouble dy = poly->points[i].y - lowest_pt.y;
+          gdouble dx = poly->points[i].x - lowest_pt.x;
+
+          if (dy == 0 && dx == 0)
+            {
+              j--;
+              num_new--;
+              continue;
+            }
+          sort_points[j].point = poly->points[i];
+          sort_points[j].angle = atan2 (dy, dx);
+        }
+    }
+
+  qsort (sort_points, num_new - 1, sizeof (SortPoint), ipolygon_sort_func);
+
+  /* now ensure that all turns as we trace the perimeter are
+     counter-clockwise */
+
+  new_points[0] = lowest_pt;
+  new_points[1] = sort_points[0].point;
+  x1 = new_points[1].x - new_points[0].x;
+  y1 = new_points[1].y - new_points[0].y;
+
+  for (i = 1, j = 2; j < num_new; i++, j++)
+    {
+      x2 = sort_points[i].point.x - new_points[j - 1].x;
+      y2 = sort_points[i].point.y - new_points[j - 1].y;
+
+      if (x2 == 0 && y2 == 0)
+        {
+          num_new--;
+          j--;
+          continue;
+        }
+
+      while (x1 * y2 - x2 * y1 < 0) /* clockwise rotation */
+        {
+          num_new--;
+          j--;
+          x1 = new_points[j - 1].x - new_points[j - 2].x;
+          y1 = new_points[j - 1].y - new_points[j - 2].y;
+          x2 = sort_points[i].point.x - new_points[j - 1].x;
+          y2 = sort_points[i].point.y - new_points[j - 1].y;
+        }
+      new_points[j] = sort_points[i].point;
+      x1 = x2;
+      y1 = y2;
+    }
+
+  g_free (sort_points);
+
+  new_poly->npoints = num_new;
+
+  return new_poly;
+}
+
+/* Determines whether a specified point is in the given polygon.
+   Based on
+
+   inpoly.c by Bob Stein and Craig Yap.
+
+   (Linux Journal, Issue 35 (March 1997), p 68)
+   */
+
+gint
+ipolygon_contains (IPolygon *poly,
+                   gint      xt,
+                   gint      yt)
+{
+  gint xnew, ynew;
+  gint xold, yold;
+  gint x1,y1;
+  gint x2,y2;
+
+  gint i;
+  gint inside = 0;
+
+  if (poly->npoints < 3)
+    return 0;
+
+  xold=poly->points[poly->npoints - 1].x;
+  yold=poly->points[poly->npoints - 1].y;
+  for (i = 0; i < poly->npoints; i++)
+    {
+      xnew = poly->points[i].x;
+      ynew = poly->points[i].y;
+      if (xnew > xold)
+        {
+          x1 = xold;
+          x2 = xnew;
+          y1 = yold;
+          y2 = ynew;
+        }
+      else
+        {
+          x1 = xnew;
+          x2 = xold;
+          y1 = ynew;
+          y2 = yold;
+        }
+      if ((xnew < xt) == (xt <= xold) &&
+          (yt - y1)*(x2 - x1) < (y2 - y1)*(xt - x1))
+        inside = !inside;
+      xold = xnew;
+      yold = ynew;
+    }
+  return inside;
+}
+
+void
+aff_element_compute_color_trans (AffElement *elem)
+{
+  int i, j;
+
+  if (elem->v.simple_color)
+    {
+      gdouble mag2;
+
+      mag2  = SQR (elem->v.target_color.r);
+      mag2 += SQR (elem->v.target_color.g);
+      mag2 += SQR (elem->v.target_color.b);
+
+      /* For mag2 == 0, the transformation blows up in general
+         but is well defined for hue_scale == value_scale, so
+         we assume that special case. */
+      if (mag2 == 0)
+        for (i = 0; i < 3; i++)
+          {
+            for (j = 0; j < 4; j++)
+              elem->color_trans.vals[i][j] = 0.0;
+
+            elem->color_trans.vals[i][i] = elem->v.hue_scale;
+          }
+      else
+        {
+          /*  red  */
+          for (j = 0; j < 3; j++)
+            {
+              elem->color_trans.vals[0][j] = elem->v.target_color.r
+                / mag2 * (elem->v.value_scale - elem->v.hue_scale);
+            }
+
+          /*  green  */
+          for (j = 0; j < 3; j++)
+            {
+              elem->color_trans.vals[1][j] = elem->v.target_color.g
+                / mag2 * (elem->v.value_scale - elem->v.hue_scale);
+            }
+
+          /*  blue  */
+          for (j = 0; j < 3; j++)
+            {
+              elem->color_trans.vals[2][j] = elem->v.target_color.g
+                / mag2 * (elem->v.value_scale - elem->v.hue_scale);
+            }
+
+          elem->color_trans.vals[0][0] += elem->v.hue_scale;
+          elem->color_trans.vals[1][1] += elem->v.hue_scale;
+          elem->color_trans.vals[2][2] += elem->v.hue_scale;
+
+          elem->color_trans.vals[0][3] =
+            (1 - elem->v.value_scale) * elem->v.target_color.r;
+          elem->color_trans.vals[1][3] =
+            (1 - elem->v.value_scale) * elem->v.target_color.g;
+          elem->color_trans.vals[2][3] =
+            (1 - elem->v.value_scale) * elem->v.target_color.b;
+
+          }
+
+
+      aff3_apply (&elem->color_trans, 1.0, 0.0, 0.0,
+                  &elem->v.red_color.r,
+                  &elem->v.red_color.g,
+                  &elem->v.red_color.b);
+      aff3_apply (&elem->color_trans, 0.0, 1.0, 0.0,
+                  &elem->v.green_color.r,
+                  &elem->v.green_color.g,
+                  &elem->v.green_color.b);
+      aff3_apply (&elem->color_trans, 0.0, 0.0, 1.0,
+                  &elem->v.blue_color.r,
+                  &elem->v.blue_color.g,
+                  &elem->v.blue_color.b);
+      aff3_apply (&elem->color_trans, 0.0, 0.0, 0.0,
+                  &elem->v.black_color.r,
+                  &elem->v.black_color.g,
+                  &elem->v.black_color.b);
+    }
+  else
+    {
+      elem->color_trans.vals[0][0] =
+        elem->v.red_color.r - elem->v.black_color.r;
+      elem->color_trans.vals[1][0] =
+        elem->v.red_color.g - elem->v.black_color.g;
+      elem->color_trans.vals[2][0] =
+        elem->v.red_color.b - elem->v.black_color.b;
+
+      elem->color_trans.vals[0][1] =
+        elem->v.green_color.r - elem->v.black_color.r;
+      elem->color_trans.vals[1][1] =
+        elem->v.green_color.g - elem->v.black_color.g;
+      elem->color_trans.vals[2][1] =
+        elem->v.green_color.b - elem->v.black_color.b;
+
+      elem->color_trans.vals[0][2] =
+        elem->v.blue_color.r - elem->v.black_color.r;
+      elem->color_trans.vals[1][2] =
+        elem->v.blue_color.g - elem->v.black_color.g;
+      elem->color_trans.vals[2][2] =
+        elem->v.blue_color.b - elem->v.black_color.b;
+
+      elem->color_trans.vals[0][3] = elem->v.black_color.r;
+      elem->color_trans.vals[1][3] = elem->v.black_color.g;
+      elem->color_trans.vals[2][3] = elem->v.black_color.b;
+    }
+}
+
+void
+aff_element_compute_trans (AffElement *elem,
+                           gdouble     width,
+                           gdouble     height,
+                           gdouble     center_x,
+                           gdouble     center_y)
+{
+  Aff2 t1, t2, t3;
+
+  /* create the rotation, scaling and shearing part of the transform */
+  aff2_distort (&t1, elem->v.asym, elem->v.shear);
+  aff2_scale (&t2, elem->v.scale, elem->v.flip);
+  aff2_compose (&t3, &t2, &t1);
+  aff2_rotate (&t2, elem->v.theta);
+  aff2_compose (&t1, &t2, &t3);
+
+  /* now create the translational part */
+  aff2_translate (&t2, -center_x*width, -center_y*width);
+  aff2_compose (&t3, &t1, &t2);
+  aff2_translate (&t2, elem->v.x*width, elem->v.y*width);
+  aff2_compose (&elem->trans, &t2, &t3);
+}
+
+void
+aff_element_decompose_trans (AffElement *elem,
+                             Aff2       *aff,
+                             gdouble     width,
+                             gdouble     height,
+                             gdouble     center_x,
+                             gdouble     center_y)
+{
+  Aff2    t1, t2;
+  gdouble det, scale, sign;
+
+  /* pull of the translational parts */
+  aff2_translate (&t1,center_x * width, center_y * width);
+  aff2_compose (&t2, aff, &t1);
+
+  elem->v.x = t2.b1 / width;
+  elem->v.y = t2.b2 / width;
+
+  det = t2.a11 * t2.a22 - t2.a12 * t2.a21;
+
+  if (det == 0.0)
+    {
+      elem->v.scale = 0.0;
+      elem->v.theta = 0.0;
+      elem->v.asym  = 1.0;
+      elem->v.shear = 0.0;
+      elem->v.flip  = 0;
+    }
+  else
+    {
+      if (det >= 0)
+        {
+          scale = elem->v.scale = sqrt (det);
+          sign = 1;
+          elem->v.flip = 0;
+        }
+      else
+        {
+          scale = elem->v.scale = sqrt (-det);
+          sign = -1;
+          elem->v.flip = 1;
+        }
+
+      elem->v.theta = atan2 (-t2.a21, t2.a11);
+
+      if (cos (elem->v.theta) == 0.0)
+        {
+          elem->v.asym = - t2.a21 / scale / sin (elem->v.theta);
+          elem->v.shear = - sign * t2.a22 / scale / sin (elem->v.theta);
+        }
+      else
+        {
+          elem->v.asym = sign * t2.a11 / scale / cos (elem->v.theta);
+          elem->v.shear = sign *
+            (t2.a12/scale - sin (elem->v.theta)/elem->v.asym)
+            / cos (elem->v.theta);
+        }
+    }
+}
+
+static void
+aff_element_compute_click_boundary (AffElement *elem,
+                                    int         num_elements,
+                                    gdouble    *points_x,
+                                    gdouble    *points_y)
+{
+  gint    i;
+  gdouble xtot = 0;
+  gdouble ytot = 0;
+  gdouble xc, yc;
+  gdouble theta;
+  gdouble sth, cth;              /* sin(theta), cos(theta) */
+  gdouble axis1, axis2;
+  gdouble axis1max, axis2max, axis1min, axis2min;
+
+  /* compute the center of mass of the points */
+  for (i = 0; i < num_elements; i++)
+    {
+      xtot += points_x[i];
+      ytot += points_y[i];
+    }
+  xc = xtot / num_elements;
+  yc = ytot / num_elements;
+
+  /* compute the sum of the (x+iy)^2, and take half the the resulting
+     angle (xtot+iytot = A*exp(2i*theta)), to get an average direction */
+
+  xtot = 0;
+  ytot = 0;
+  for (i = 0; i < num_elements; i++)
+    {
+      xtot += SQR (points_x[i] - xc) - SQR (points_y[i] - yc);
+      ytot += 2 * (points_x[i] - xc) * (points_y[i] - yc);
+    }
+  theta = 0.5 * atan2 (ytot, xtot);
+  sth = sin (theta);
+  cth = cos (theta);
+
+  /* compute the minimum rectangle at angle theta that bounds the points,
+     1/2 side lengths left in axis1, axis2, center in xc, yc */
+
+  axis1max = axis1min = 0.0;
+  axis2max = axis2min = 0.0;
+  for (i = 0; i < num_elements; i++)
+    {
+      gdouble proj1 =  (points_x[i] - xc) * cth + (points_y[i] - yc) * sth;
+      gdouble proj2 = -(points_x[i] - xc) * sth + (points_y[i] - yc) * cth;
+      if (proj1 < axis1min)
+        axis1min = proj1;
+      if (proj1 > axis1max)
+        axis1max = proj1;
+      if (proj2 < axis2min)
+        axis2min = proj2;
+      if (proj2 > axis2max)
+        axis2max = proj2;
+    }
+  axis1 = 0.5 * (axis1max - axis1min);
+  axis2 = 0.5 * (axis2max - axis2min);
+  xc += 0.5 * ((axis1max + axis1min) * cth - (axis2max + axis2min) * sth);
+  yc += 0.5 * ((axis1max + axis1min) * sth + (axis2max + axis2min) * cth);
+
+  /* if the the rectangle is less than 10 pixels in any dimension,
+     make it click_boundary, otherwise set click_boundary = draw_boundary */
+
+  if (axis1 < 8.0 || axis2 < 8.0)
+    {
+      GdkPoint *points = g_new (GdkPoint, 4);
+
+      elem->click_boundary = g_new (IPolygon, 1);
+      elem->click_boundary->points = points;
+      elem->click_boundary->npoints = 4;
+
+      if (axis1 < 8.0) axis1 = 8.0;
+      if (axis2 < 8.0) axis2 = 8.0;
+
+      points[0].x = xc + axis1 * cth - axis2 * sth;
+      points[0].y = yc + axis1 * sth + axis2 * cth;
+      points[1].x = xc - axis1 * cth - axis2 * sth;
+      points[1].y = yc - axis1 * sth + axis2 * cth;
+      points[2].x = xc - axis1 * cth + axis2 * sth;
+      points[2].y = yc - axis1 * sth - axis2 * cth;
+      points[3].x = xc + axis1 * cth + axis2 * sth;
+      points[3].y = yc + axis1 * sth - axis2 * cth;
+    }
+  else
+    elem->click_boundary = elem->draw_boundary;
+}
+
+void
+aff_element_compute_boundary (AffElement  *elem,
+                              gint         width,
+                              gint         height,
+                              AffElement **elements,
+                              gint         num_elements)
+{
+  gint      i;
+  IPolygon  tmp_poly;
+  gdouble  *points_x;
+  gdouble  *points_y;
+
+  if (elem->click_boundary && elem->click_boundary != elem->draw_boundary)
+    g_free (elem->click_boundary);
+  if (elem->draw_boundary)
+    g_free (elem->draw_boundary);
+
+  tmp_poly.npoints = num_elements;
+  tmp_poly.points  = g_new (GdkPoint, num_elements);
+  points_x = g_new (gdouble, num_elements);
+  points_y = g_new (gdouble, num_elements);
+
+  for (i = 0; i < num_elements; i++)
+    {
+      aff2_apply (&elem->trans,
+                  elements[i]->v.x * width, elements[i]->v.y * width,
+                  &points_x[i],&points_y[i]);
+      tmp_poly.points[i].x = (gint)points_x[i];
+      tmp_poly.points[i].y = (gint)points_y[i];
+    }
+
+  elem->draw_boundary = ipolygon_convex_hull (&tmp_poly);
+  aff_element_compute_click_boundary (elem, num_elements, points_x, points_y);
+
+  g_free (tmp_poly.points);
+}
+
+void
+aff_element_draw (AffElement  *elem,
+                  gboolean     selected,
+                  gint         width,
+                  gint         height,
+                  cairo_t     *cr,
+                  GdkColor    *color,
+                  PangoLayout *layout)
+{
+  PangoRectangle rect;
+  gint           i;
+
+  pango_layout_set_text (layout, elem->name, -1);
+  pango_layout_get_pixel_extents (layout, NULL, &rect);
+
+  gdk_cairo_set_source_color (cr, color);
+
+  cairo_move_to (cr,
+                 elem->v.x * width - rect.width  / 2,
+                 elem->v.y * width + rect.height / 2);
+  pango_cairo_show_layout (cr, layout);
+  cairo_fill (cr);
+
+  cairo_set_line_width (cr, 1.0);
+
+  if (elem->click_boundary != elem->draw_boundary)
+    {
+      cairo_move_to (cr,
+                     elem->click_boundary->points[0].x,
+                     elem->click_boundary->points[0].y);
+
+      for (i = 1; i < elem->click_boundary->npoints; i++)
+        cairo_line_to (cr,
+                       elem->click_boundary->points[i].x,
+                       elem->click_boundary->points[i].y);
+
+      cairo_close_path (cr);
+
+      cairo_stroke (cr);
+    }
+
+  if (selected)
+    cairo_set_line_width (cr, 3.0);
+
+  cairo_move_to (cr,
+                 elem->draw_boundary->points[0].x,
+                 elem->draw_boundary->points[0].y);
+
+  for (i = 1; i < elem->draw_boundary->npoints; i++)
+    cairo_line_to (cr,
+                   elem->draw_boundary->points[i].x,
+                   elem->draw_boundary->points[i].y);
+
+  cairo_close_path (cr);
+
+  cairo_stroke (cr);
+}
+
+AffElement *
+aff_element_new (gdouble  x,
+                 gdouble  y,
+                 GimpRGB *color,
+                 gint     count)
+{
+  AffElement *elem = g_new (AffElement, 1);
+  gchar buffer[16];
+
+  elem->v.x = x;
+  elem->v.y = y;
+  elem->v.theta = 0.0;
+  elem->v.scale = 0.5;
+  elem->v.asym = 1.0;
+  elem->v.shear = 0.0;
+  elem->v.flip = 0;
+
+  elem->v.red_color   = *color;
+  elem->v.blue_color  = *color;
+  elem->v.green_color = *color;
+  elem->v.black_color = *color;
+
+  elem->v.target_color = *color;
+  elem->v.hue_scale = 0.5;
+  elem->v.value_scale = 0.5;
+
+  elem->v.simple_color = TRUE;
+
+  elem->draw_boundary = NULL;
+  elem->click_boundary = NULL;
+
+  aff_element_compute_color_trans (elem);
+
+  elem->v.prob = 1.0;
+
+  sprintf (buffer,"%d", count);
+  elem->name = g_strdup (buffer);
+
+  return elem;
+}
+
+void
+aff_element_free (AffElement *elem)
+{
+  if (elem->click_boundary != elem->draw_boundary)
+    g_free (elem->click_boundary);
+
+  g_free (elem->draw_boundary);
+  g_free (elem);
+}
+
+#ifdef DEBUG_BRUSH
+static brush_chars[] = {' ',':','*','@'};
+#endif
+
+static guchar *
+create_brush (IfsComposeVals *ifsvals,
+              gint           *brush_size,
+              gdouble        *brush_offset)
+{
+  gint     i, j;
+  gint     ii, jj;
+  guchar  *brush;
+#ifdef DEBUG_BRUSH
+  gdouble  totpix = 0.0;
+#endif
+
+  gdouble radius = ifsvals->radius * ifsvals->subdivide;
+
+  *brush_size = ceil (2 * radius);
+  *brush_offset = 0.5 * (*brush_size - 1);
+
+  brush = g_new (guchar, SQR (*brush_size));
+
+  for (i = 0; i < *brush_size; i++)
+    {
+      for (j = 0; j < *brush_size; j++)
+        {
+          gdouble pixel = 0.0;
+          gdouble d     = sqrt (SQR (i - *brush_offset) +
+                                SQR (j - *brush_offset));
+
+          if (d - 0.5 * G_SQRT2 > radius)
+            pixel = 0.0;
+          else if (d + 0.5 * G_SQRT2 < radius)
+            pixel = 1.0;
+          else
+            for (ii = 0; ii < 10; ii++)
+              for (jj = 0; jj < 10; jj++)
+                {
+                  d = sqrt (SQR (i - *brush_offset + ii * 0.1 - 0.45) +
+                            SQR (j - *brush_offset + jj * 0.1 - 0.45));
+                  pixel += (d < radius) / 100.0;
+                }
+
+          brush[i * *brush_size + j] = 255.999 * pixel;
+
+#ifdef DEBUG_BRUSH
+          putchar(brush_chars[(gint)(pixel * 3.999)]);
+          totpix += pixel;
+#endif /* DEBUG_BRUSH */
+        }
+#ifdef DEBUG_BRUSH
+      putchar('\n');
+#endif /* DEBUG_BRUSH */
+    }
+#ifdef DEBUG_BRUSH
+  printf ("Brush total / area = %f\n", totpix / SQR (ifsvals->subdivide));
+#endif /* DEBUG_BRUSH */
+  return brush;
+}
+
+void
+ifs_render (AffElement     **elements,
+            gint             num_elements,
+            gint             width,
+            gint             height,
+            gint             nsteps,
+            IfsComposeVals  *vals,
+            gint             band_y,
+            gint             band_height,
+            guchar          *data,
+            guchar          *mask,
+            guchar          *nhits,
+            gboolean         preview)
+{
+  gint     i, k, n;
+  gdouble  x, y;
+  gdouble  r, g, b;
+  gint     ri, gi, bi;
+  guint32  p0, psum;
+  gdouble  pt;
+  guchar  *ptr;
+  guint32 *prob;
+  gdouble *fprob;
+  gint     subdivide;
+  guchar  *brush = NULL;
+  gint     brush_size   = 1;
+  gdouble  brush_offset = 0.0;
+
+  if (preview)
+    subdivide = 1;
+  else
+    subdivide = vals->subdivide;
+
+  /* compute the probabilities and transforms */
+  fprob = g_new (gdouble, num_elements);
+  prob = g_new (guint32, num_elements);
+  pt = 0.0;
+
+  for (i = 0; i < num_elements; i++)
+    {
+      aff_element_compute_trans(elements[i],
+                                width * subdivide,
+                                height * subdivide,
+                                vals->center_x,
+                                vals->center_y);
+      fprob[i] = fabs(
+          elements[i]->trans.a11 * elements[i]->trans.a22
+        - elements[i]->trans.a12 * elements[i]->trans.a21);
+
+      /* As a heuristic, if the determinant is really small, it's
+         probably a line element, so increase the probability so
+         it gets rendered */
+
+      /* FIXME: figure out what 0.01 really should be */
+      if (fprob[i] < 0.01)
+        fprob[i] = 0.01;
+
+      fprob[i] *= elements[i]->v.prob;
+
+      pt += fprob[i];
+    }
+
+  psum = 0;
+  for (i = 0; i < num_elements; i++)
+    {
+      psum += (guint32) -1 * (fprob[i] / pt);
+      prob[i] = psum;
+    }
+
+  prob[i - 1] = (guint32) -1;  /* make sure we don't get bitten by roundoff */
+
+  /* create the brush */
+  if (!preview)
+    brush = create_brush (vals, &brush_size, &brush_offset);
+
+  x = y = 0;
+  r = g = b = 0;
+
+  /* n is used to limit the number of progress updates */
+  n = nsteps / 32;
+
+  /* now run the iteration */
+  for (i = 0; i < nsteps; i++)
+    {
+      if (!preview && ((i % n) == 0))
+        gimp_progress_update ((gdouble) i / (gdouble) nsteps);
+
+      p0 = g_random_int ();
+      k = 0;
+
+      while (p0 > prob[k])
+        k++;
+
+      aff2_apply (&elements[k]->trans, x, y, &x, &y);
+      aff3_apply (&elements[k]->color_trans, r, g, b, &r, &g, &b);
+
+      if (i < 50)
+        continue;
+
+      ri = (gint) (255.0 * r + 0.5);
+      gi = (gint) (255.0 * g + 0.5);
+      bi = (gint) (255.0 * b + 0.5);
+
+      if ((ri < 0) || (ri > 255) ||
+          (gi < 0) || (gi > 255) ||
+          (bi < 0) || (bi > 255))
+        continue;
+
+      if (preview)
+        {
+          if ((x < width) && (y < (band_y + band_height)) &&
+              (x >= 0) && (y >= band_y))
+            {
+              ptr = data + 3 * (((gint) (y - band_y)) * width + (gint) x);
+
+              *ptr++ = ri;
+              *ptr++ = gi;
+              *ptr   = bi;
+            }
+        }
+      else
+        {
+          if ((x < width * subdivide) && (y < height * subdivide) &&
+              (x >= 0) && (y >= 0))
+            {
+              gint ii;
+              gint jj;
+              gint jj0   = floor (y - brush_offset - band_y * subdivide);
+              gint ii0   = floor (x - brush_offset);
+              gint jjmin = 0;
+              gint iimin = 0;
+              gint jjmax;
+              gint iimax;
+
+              if (ii0 < 0)
+                iimin = - ii0;
+              else
+                iimin = 0;
+
+              if (jj0 < 0)
+                jjmin = - jj0;
+              else
+                jjmin = 0;
+
+              if (jj0 + brush_size >= subdivide * band_height)
+                jjmax = subdivide * band_height - jj0;
+              else
+                jjmax = brush_size;
+
+              if (ii0 + brush_size >= subdivide * width)
+                iimax = subdivide * width - ii0;
+              else
+                iimax = brush_size;
+
+              for (jj = jjmin; jj < jjmax; jj++)
+                for (ii = iimin; ii < iimax; ii++)
+                  {
+                    guint m_old;
+                    guint m_new;
+                    guint m_pix;
+                    guint n_hits;
+                    guint old_scale;
+                    guint pix_scale;
+                    gint  index = (jj0 + jj) * width * subdivide + ii0 + ii;
+
+                    n_hits = nhits[index];
+                    if (n_hits == 255)
+                      continue;
+
+                    m_pix = brush[jj * brush_size + ii];
+                    if (!m_pix)
+                      continue;
+
+                    nhits[index] = ++n_hits;
+                    m_old = mask[index];
+                    m_new = m_old + m_pix - m_old * m_pix / 255;
+                    mask[index] = m_new;
+
+                    /* relative probability that old colored pixel is on top */
+                    old_scale = m_old * (255 * n_hits - m_pix);
+
+                    /* relative probability that new colored pixel is on top */
+                    pix_scale = m_pix * ((255 - m_old) * n_hits + m_old);
+
+                    ptr = data + 3 * index;
+
+                    *ptr = ((old_scale * (*ptr) + pix_scale * ri) /
+                            (old_scale + pix_scale));
+                    ptr++;
+
+                    *ptr = ((old_scale * (*ptr) + pix_scale * gi) /
+                            (old_scale + pix_scale));
+                    ptr++;
+
+                    *ptr = ((old_scale * (*ptr) + pix_scale * bi) /
+                            (old_scale + pix_scale));
+                  }
+            }
+        }
+    } /* main iteration */
+
+  if (!preview )
+    gimp_progress_update (1.0);
+
+  g_free (brush);
+  g_free (prob);
+  g_free (fprob);
+}