/* The GIMP -- an 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 2 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, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#ifdef __GNUC__
#warning FIXME: GDK_DISABLE_DEPRECATED
#endif

#undef GDK_DISABLE_DEPRECATED

#include <gdk/gdk.h>

#include <libgimp/gimp.h>

#include "ifscompose.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,
	   gint     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 directon 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)
{
  double xt = t->vals[0][0]*x + t->vals[0][1]*y + t->vals[0][2]*z + t->vals[0][3];
  double yt = t->vals[1][0]*x + t->vals[1][1]*y + t->vals[1][2]*z + t->vals[1][3];
  double 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;
      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 perimiter 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 lenghs 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)
{
  int 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,
		  gint         selected,
		  gint         width,
		  gint         height, 
		  GdkDrawable *win,
		  GdkGC       *normal_gc,
		  GdkGC       *selected_gc,
		  GdkFont     *font)
{
  GdkGC *gc;
  gint string_width = gdk_string_width (font,elem->name);
  gint string_height = font->ascent  + font->descent + 2;

  if (selected)
    gc = selected_gc;
  else
    gc = normal_gc;

  gdk_draw_string(win,font,gc,
		  elem->v.x*width-string_width/2,
		  elem->v.y*width+string_height/2,elem->name);

  if (elem->click_boundary != elem->draw_boundary)
    gdk_draw_polygon(win,normal_gc,FALSE,elem->click_boundary->points,
		     elem->click_boundary->npoints);

  gdk_draw_polygon(win,gc,FALSE,elem->draw_boundary->points,
		   elem->draw_boundary->npoints);
}

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 d = sqrt(SQR(i-*brush_offset)+SQR(j-*brush_offset));
	  gdouble pixel = 0.0;

	  if (d-0.5*sqrt(2) > radius)
	    pixel = 0.0;
	  else if (d+0.5*sqrt(2) < 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[(int)(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,
	    gint             preview)
{
  gint i,k;
  gdouble x,y;
  gdouble r,g,b;
  gint ri,gi,bi;
  gint p0,psum;
  gdouble pt;
  guchar *ptr;
  gint *prob;
  gdouble *fprob;
  gint subdivide;

  guchar *brush = NULL;
  gint brush_size;
  gdouble brush_offset;

  if (preview)
    subdivide = 1;
  else
    subdivide = vals->subdivide;
  
  /* compute the probabilities and transforms */
  fprob = g_new(gdouble,num_elements);
  prob = g_new(gint,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 += G_MAXRAND * (fprob[i]/pt);
      prob[i] = psum;
    }
  prob[i-1] = G_MAXRAND;	/* 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;

  /* now run the iteration */
  for (i=0;i<nsteps;i++) 
    {
      if (!preview && !(i % 5000))
	gimp_progress_update ((gdouble) i / (gdouble) nsteps);
      p0 = rand();
      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.999*r);
      gi = (gint)(255.999*g);
      bi = (gint)(255.999*b);
      
      if (preview)
	{
	  if ((x<width) && (y<(band_y+band_height)) && 
	      (x >= 0) && (y >= band_y) &&
	      (ri >= 0) && (ri < 256) &&
	      (gi >= 0) && (gi < 256) &&
	      (bi >= 0) && (bi < 256))
	    {
	      ptr = data + 3 * (((gint)(y-band_y))*width + (gint)x);
	      *ptr++ = ri;
	      *ptr++ = gi;
	      *ptr = bi;
	    }
	}
      else
	if ((ri >= 0) && (ri < 256) &&
	    (gi >= 0) && (gi < 256) &&
	    (bi >= 0) && (bi < 256))
	  {
	    guint m_old;
	    guint m_new;
	    guint m_pix;
	    guint n_hits;
	    guint old_scale;
	    guint pix_scale;
	    
	    gint index;
	    gint ii,jj;
	    gint jj0 = floor(y-brush_offset-band_y*subdivide);
	    gint ii0 = floor(x-brush_offset);
	    gint jjmax,iimax;
	    gint jjmin = 0;
	    gint iimin = 0;

	    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++)
		{
		  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 (brush)
    g_free(brush);
  g_free(prob);
  g_free(fprob);
}