Iterated function system fractal

**OnionKnight** · 01-12-2007

I read a book about fractals and it had something that seemed to be a simple way of drawing fractals called IFS so I was intrigued and wanted to try it out but I seem to have gotten it wrong and can't find any good information on how to implement it.
As I understood it you start at a point of origin, which i chose to be the middle of the screen in my program, and then apply a transformation on that point using a randomly selected row from a matrix that works like:
X' = aX + bY + e
Y' = cX + dY + f
and then you plot out this coordinate and repeat a bunch of times.
In my program I'm using a matrix that should generate a slanted sierpinski triangle but all I get is dots in a diagonal way. I'd post a screenshot but I can't right now because I'm VNCing and can't take screenshots.

Code:

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <mt19937ar.h>
#include <SDL/SDL.h>

#define WIDTH  640
#define HEIGHT 480
#define ITER 9000
#define COLOR 0x00ff00
#define PITCH  (screen->pitch/4)
#define PIXELS ((unsigned int*) screen->pixels)

typedef struct {
	double a, b, c, d, e, f;
} ifs_elem;

void draw_frac (ifs_elem* ifs, size_t n);

SDL_Surface* screen;

int main (int argc, char** argv)
{
	SDL_Event event;
	ifs_elem sierp[3] = {{0.5, 0.0, 0.0, 0.5, 0.0, 0.0},
	                     {0.5, 0.0, 0.0, 0.5, 0.0, 0.5},
	                     {0.5, 0.0, 0.0, 0.5, 0.5, 0.5}};

	if (SDL_Init(SDL_INIT_EVERYTHING) < 0) {
		fprintf(stderr, "Error: Unable to init SDL: %s\n", SDL_GetError());
		return 1;
	}
	atexit(SDL_Quit);
	if ((screen = SDL_SetVideoMode(WIDTH, HEIGHT, 32, SDL_HWSURFACE)) == NULL) {
		fprintf(stderr, "Error: Unable to init SDL: %s\n", SDL_GetError());
		return 1;
	}
	init_genrand(time(NULL));

	draw_frac(sierp, sizeof(sierp)/sizeof(*sierp));
	while (SDL_WaitEvent(&event)) {
		switch (event.key.keysym.sym) {
			case SDLK_q:
			case SDLK_ESCAPE:
				return 0;
		}
		if (event.type == SDL_QUIT)
			return 0;
	}

	return 0;
}

void draw_frac (ifs_elem* ifs, size_t n)
{
	int i, j;
	int x, y, oldx;

	x = WIDTH/2;
	y = HEIGHT/2;
	if (SDL_MUSTLOCK(screen))
		SDL_LockSurface(screen);
	for (i = 0; i < ITER; i++) {
		PIXELS[x + y*PITCH] = COLOR;
		j = genrand_int32()%n;
		oldx = x;
		x = ifs[j].a*x + ifs[j].b*y + ifs[j].e;
		y = ifs[j].c*oldx + ifs[j].d*y + ifs[j].f;
	}
	if (SDL_MUSTLOCK(screen))
		SDL_UnlockSurface(screen);
	while (SDL_Flip(screen) < 0);
}

Uses SDL and MT.

**coder8137** · 01-12-2007

Have you tried:

Code:

ifs_elem sierp[3] = {{0.5, 0.0, 0.0, 0.5, 0.0, 0.0},
	                     {0.5, 0.0, 0.0, 0.5, 1.0, 0.0},
	                     {0.5, 0.0, 0.0, 0.5, 0.5, 0.5}};

This should be the transformations to draw the Sierpinski triangle, but with 0 slant.

EDIT: Also, the first point you start with should be (0,0) as far as I know. Not sure if this will affect the algorithm much. (If you use (0,0) you probably need to scale it up before plotting it.)

**OnionKnight** · 01-12-2007

Have you tried:

That's the one the Wikipedia article uses. It didn't make a difference.

EDIT: Also, the first point you start with should be (0,0) as far as I know. Not sure if this will affect the algorithm much. (If you use (0,0) you probably need to scale it up before plotting it.)

(0, 0) in math is usually the middle of a graph (negative values going down and left from it and positive to the right and up) whereas (0, 0) on the screen is the top left corner so setting origo to the middle of the screen should be correct.
Trying to set X and Y to 0 gave me a black screen (no output).

**coder8137** · 01-12-2007

(0,0) will mean that all points will fall between 0 <= x < 2 and 0 <= y < 1. Therefore your program would have probably output 2 pixels, since you don't use floats. Hence the need to scale it up.

Anyways, I tried it in OpenGL. It looks like your initial transformations produce the triangle at 45 degrees, like you said. The ones I suggested work too, with 0 slant. Starting point doesn't seem to affect point distribution notably either. I plotted 1500000 points for the screen shot. I am not familiar with SDL, but all the math seems to be right. Here's the code (keep in mind that my OpenGL knowledge is still beginner level):

Code:

#include <stdlib.h>
#include <time.h>
#include <GL/glut.h>

#define NOS_POINTS 1500000
#define WIN_WIDTH 1024
#define WIN_HEIGHT 512


void display()
{
  int i;
  double matrixes[3][6] = {{0.5, 0.0, 0.0, 0.5, 0.0, 0.0},
	                   {0.5, 0.0, 0.0, 0.5, 0.0, 0.5},
	                   {0.5, 0.0, 0.0, 0.5, 0.5, 0.5}};
  double p[2] = {1, 0.5};
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  glBegin(GL_POINTS);
  for(i=0; i<NOS_POINTS; i++) {
    int j;
    double old_p0;
    glVertex2dv(p); /* draw the point p */
    j = rand()/((double)RAND_MAX + 1)*3;
    old_p0 = p[0];
    p[0] = p[0]*matrixes[j][0] + p[1]*matrixes[j][1] + matrixes[j][4];
    p[1] = old_p0*matrixes[j][2] + p[1]*matrixes[j][3] + matrixes[j][5];
  }
  glEnd();
  glFlush();
}

void myinit() {
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  glOrtho(0, 2, 0, 1, -1, 1);      /* Set up viewing area */
  glMatrixMode(GL_MODELVIEW);
  glPointSize(WIN_WIDTH*WIN_HEIGHT/NOS_POINTS);
  glColor3f(0.0, 0.0, 1.0);
}

int main(int argc, char **argv)
{
    srand(time(NULL));
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB | GLUT_DEPTH);
    glutInitWindowSize(WIN_WIDTH, WIN_HEIGHT);
    glutCreateWindow("Sierpinski Triangle");
    glutDisplayFunc(display);
    glEnable(GL_DEPTH_TEST);
    glClearColor (1.0, 1.0, 1.0, 1.0);
    myinit();
    glutMainLoop();
    return 0;
}

And the output (Note, you can see the command I used to make it on a linux machine. The code should work on Win too, if GLUT / freeglut is installed, but the command is different of course.):

**coder8137** · 01-12-2007

Actually, I'll take that back about your math being correct.
For example, you start with:

Code:

x=WIDTH/2 = 640/2 = 320

Now, taking the transformation 3, the next few x points you plot will be:

Code:

x = 320
x = 320 * 0.5 + 0.5 = 160.5
x = 160.5 * 0.5 + 0.5 = 80.75
x = 40.876
x = 20.9375
x = 10.96875
/* same thing happens with y */
/* this does not allow for all the truncation you would get because you use ints */

The result is that there is no way for the points (regardless of which transformation is chose randomly), to do anything else but go towards 0 way too fast.

You need to either use (0, 0) as a start point and plot the area 0<=x < 2 and 0<=y<=1, OR you need to increase the 0.5 you use as a constant. I.e. the value of e and f in the structure must be alot bigger. Probably 640 / 4 for e and for f: 480 / 2.

**quzah** · 01-12-2007

You're compiling as root? Hang on, I want to whip up a piece of code for you...

Quzah.

**OnionKnight** · 01-13-2007

Changing the constants so that they multiply with the width and height did the trick. The initial values doesn't seem to make a difference so I set them to 0.

Code:

	ifs_elem sierp[] = {{0.5, 0.0, 0.0, 0.5, 0.0*WIDTH, 0.0*HEIGHT},
	                    {0.5, 0.0, 0.0, 0.5, 0.0*WIDTH, 0.5*HEIGHT},
	                    {0.5, 0.0, 0.0, 0.5, 0.5*WIDTH, 0.5*HEIGHT}};

But, when I try to draw the non-slanted version I get an unexpected result as shown in attatched screenshot.

Code:

	ifs_elem sierp[] = {{0.5, 0.0, 0.0, 0.5, 0.0*WIDTH, 0.0*HEIGHT},
	                    {0.5, 0.0, 0.0, 0.5, 1.0*WIDTH, 0.0*HEIGHT},
	                    {0.5, 0.0, 0.0, 0.5, 0.5*WIDTH, 0.5*HEIGHT}};

Drawing something else, like the fern gives some crazy result but I'm not gonna tackle that one before I know what's wrong with the triangle first. The problems are probably linked.

**coder8137** · 01-13-2007

You're compiling as root?

hehe, didn't feel like switching users then.

Anyways, the transformations should probably be:

Code:

ifs_elem sierp[] = {{0.5, 0.0, 0.0, 0.5, 0.0*WIDTH / 2, 0.0*HEIGHT},
	                    {0.5, 0.0, 0.0, 0.5, 1.0*WIDTH / 2, 0.0*HEIGHT},
	                    {0.5, 0.0, 0.0, 0.5, 0.5*WIDTH / 2, 0.5*HEIGHT}};

EDIT: Forgot to mention that your triangle is still upside down. Try plotting it like this:

Code:

PIXELS[x + (HEIGHT - y)*PITCH] = COLOR;

**OnionKnight** · 01-13-2007

Deviding with 2 solved it. But I don't understand why. The reason I left it out before was that it made the first triangle (the slanted one) to take up only half the screen. I'm not sure why it becomes upside down either.

As for the fern I tried using the IFS schemes found on this page: http://local.wasp.uwa.edu.au/~pbourk...ls/ifs_fern_a/
However they both get weird. The first one seems to be able to draw the skeleton or whatever it's called correctly. I noticed that they both go off screen as some of the constants are 1.6, which would translate to 1.6*HEIGHT in my program. I tried dividing by 2 but that didn't change much. I guess I should ask what these constans actually represent and what their ranges are so I can scale them.
Earlier you said

(0,0) will mean that all points will fall between 0 <= x < 2 and 0 <= y < 1.

but I don't really understand that, especially why X has double the range of Y, unless that was specific to your program and not IFS.

**coder8137** · 01-13-2007

The easiest way to understand everything, is if you consider the program as 2 seperate steps. At the moment, you have been mixing them together, making things a lot harder.

Step 1: Calculating the Points

This is just about applying the formulas and getting a bunch of x and y values, which we can then use for Step 2. For this step, you should really be using the original transformation matrixes, which we had at the start (i.e. use 0.5, rather than 0.5 * HEIGHT). Everything I'm saying from now on refers to the NON-slanted triangle, using constants ranging from 0 - 1.

Now, you can theoretically start with any starting point, and the result should be the same (at least for our equations). However, for Step 2, you need to know (roughly) the area to plot, in terms of Xmin, Xmax, Ymin, Ymax. In order to find out what these values are, it's probably a good idea (though not necessary) to put (0,0) into the transformations, iterate a few times, and see what happens. For example, to slove for Xmax, consider using transformation 2 repeatedly, as this is the one that will increase x the most:

Code:

/* Transformation 2 for non-slanted: */
X(n+1) = 0.5 * X(n) + 1

/* Our Starting Point                */
X(0) = 0

/* Iterate                           */
X(1) = 0.5 * 0 + 1 == 1
X(2) = 0.5 * 1 + 1 == 1.5
X(3)               == 1.75
........
X(.)               == 1.992
X(. + 1)           == 1.996
/* you can keep going, but X(n) will remain less than 2,
   until you run out of precision                          */

Now, you could do the same for Ymax (< 1) and then Xmin (== 0) and Ymin (== 0).

This was the trial and error way of finding a good plotting area. Of course, you probably realized that, in this particular case, solving:

Code:

x           = 0.5 * x + 1
/* i.e. INPUT == OUTPUT. Recursion values will never change once this x is reached */

Solve for x:
x - 0.5 * x = 1
x           = 2

is the same, and a lot quicker. If you have more complex equations, checking just a single one may not suffice.

For this whole step, considering our equations, you should be using doubles, not ints, because you'll never get anything bigger than 2.

Step 2: Plotting the Points
Now, that you've created the points, you should treat plotting seperatly. We've checked that all our points fall in the range 0 <= x < 2 and 0<= y < 1 (for the 3 equations for the NON-slanted triangle (slanted one's range is slightly different)). Therefore, this is generally also the area you want to plot. Note that Xmin/max Ymin/max in this section refer to the area we intend to plot; this might sometimes be different to what we determined above (often we zoom in, cropping out points). The plotting step simply involves translating a point in the range we just mentioned, to a pixel in the range 0 <= p_X < WINDOW_WIDTH and 0 <= p_Y < WINDOW_HEIGHT. You can just use a linear translation. Thus:

Code:

p_X =  (WINDOW_WIDTH / Xrange) * (-Xmin + x)
    == (WINDOW_WIDTH / (Xmax - Xmin)) * (-Xmin + x)
    == WINDOW_WIDHT / 2 * x

p_Y =  (WINDOW_HEIGHT / Yrange) * (-Ymin + y)
    == (WINDOW_HEIGHT / (Ymax - Ymin)) * (-Ymin + y)
    == WINDOW_HEIGHT / 1 * y

/* Anything other than p_X and p_Y is in the coordinate system of Step 1. I.e. Xrange is the
   actual range of values you want to plot. It is NOT related to WINDOW_WIDTH                  */

Now, all we need to consider, is what p_X and p_Y actually refer to (i.e what are they relative to). Originally, in Step 1, we had x and y points relative to an origin, which is (0, 0), which is the center of the (general) coordinate system in math. However, the screen doesn't know anything about negative values, so what we did above, is translate points from Step 1, to measurements from the _BOTTOM_ left corner of the screen in pixels. So, p_X and p_Y refer to number of pixels from the bottom left corner. I translated the points like this, because in a normal math graph, without negatives, origin would be bottom left.

So, lastly, we need to convert p_X and p_Y to coordinates that the Windowing system will recognize. In other words, since a lot of windowing systems expect coordinates with respect to the _TOP_ left corner, we need to allow for that. If you wish, you can change the equation above, or, you can leave it as a seperate step.

Code:

p_Y_top = WINDOW_HEIGHT - p_Y

Master Caution
Like I said, I don't know much about SDL, BUT, it seems to be simply setting pixels in an array. Often times, you will want to have points outside a certain range "cropped" out (i.e. when you zoom in). In order to prevent incorrect display and also nasty buffer overruns, you probably need to check each point manually, if there is any chance that it doesn't fall into the range 0 <= p_X < WINDOW_WIDTH and 0 <= p_Y_top < WINDOW_HEIGHT.

Code:

/* This is the final part of plotting, also ensuring clipping is done */
if(p_Y_top >= 0 && p_Y_top < WINDOW_HEIGHT && p_X >= 0 && p_X < WINDOW_WIDTH)
  PIXELS[p_X + p_Y_top*PITCH] = COLOR;

**Furthermore, you may have to use WINDOW_WIDTH/HEIGHT-1, wherever I used WINDOW_WIDTH/HEIGHT, dpending on array limits.**

Furthermore, don't forget to cast macros to doubles, if integer divison would otherwise result.

Other Considerations
This should be enough info to understand why the transformation equation had to use WIDTH / 2. There's nothing really wrong with just changing the original equations, but it can soon get confusing, expecially if you also want to zoom in. I'd recommend sticking with a 2 Step approach.

The Fern
Quite simple now. Just leave the equation as per web page; design your app to handle changing of plotting area (i.e. plotting Xmax/min Ymax/min) easily and use some trial and error. This means you never really have to figure out the exact values of Xmin/max/Ymin/max in Step 1. Changing plotting Xmin/max/Ymin/max using trial & error, will eventually get you the rough values for the real Xmin/max/Ymin/max. Start with say Ymin = -10, Ymax = 10, Xmin = -10, Xmax = 10 for example. See what it looks like and chane it from there. If you feel like manual calculating the range values, that's possible too, but the equations are a little more complex than for the triangle.

Make sure to factor the probability some way or another. I extended the array and used a loop; probably should have used a bunch of else ifs, but that wuld have made the post even longer, hehe:

Code:

void display()
{
  int i;
  double matrixes[][7] = {{0.0,   0.0,   0.0,  0.16, 0.0, 0.0,  0.1},
	                  {0.2,  -0.26,  0.23, 0.22, 0.0, 1.6,  0.08},
	                  {-0.15, 0.28,  0.26, 0.24, 0.0, 0.44, 0.08},
                          {0.75,  0.04, -0.04, 0.85, 0.0, 1.6,  0.74}};
  float colors[][3]    = {{0.9, 0.1,  0.3},
                          {0.1, 0.4,  0.8},
                          {0.3, 0.75, 0.1},
                          {0.2, 0.3,  0.2}};
     /* colors array MUST have >= nos of elements as matrixes array       */
  const int nos_of_equations = sizeof(matrixes)/sizeof(matrixes[0]);
  double p[2] = {0.0, 0.0};                    /* start point             */
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  glBegin(GL_POINTS);
  for(i=0; i<NOS_POINTS; i++) {
    double old_p0 = p[0];
    double j = rand()/((double)RAND_MAX + 1); /* bet 0 and 1 (exc. 1)     */
    double temp_prob_sum;
    int i;
    for(i=0, temp_prob_sum = matrixes[0][6]; ;) {
      if(j < temp_prob_sum) {
        glVertex2dv(p);                            /* draw the point p        */
        p[0] = p[0]*matrixes[i][0] + p[1]*matrixes[i][1] + matrixes[i][4];
        p[1] = old_p0*matrixes[i][2] + p[1]*matrixes[i][3] + matrixes[i][5];
        glColor3fv(colors[i]);                     /* set color for next pt   */
        break;
      }
      if(i < nos_of_equations - 1)
        temp_prob_sum += matrixes[++i][6];
      else
        temp_prob_sum = 1.1;
        /* force the 1st if above true. We shouldn't get here. Check equation probs */
    }
  }
  glEnd();
  glFlush();
}

**OnionKnight** · 01-14-2007

I see, so the IFS codes doesn't have any set dimensions and each of them have their own min and max.

If I want the min and max values I go through each row of the matrix and calculate the min and max values for X and Y and at the end I pick out the highest maximum values I can find and the lowest minimum values?

I tried setting up an equation for a general solution to the maximums,

Code:

Xmax = a*Xmax + b*Ymax + e
Ymax = c*Xmax + d*Ymax + f

(1 - a)*Xmax = b*Ymax + e
(1 - d)*Ymax = c*Xmax + f

Xmax = (b*Ymax + e)/(1 - a)
Ymax = (c*Xmax + f)/(1 - d)

Xmax = (b*((c*Xmax + f)/(1 - d)) + e)/(1 - a)
Ymax = (c*((b*Ymax + e)/(1 - a)) + f)/(1 - d)

//I cheated using Mathematica here to solve the mess
Xmax = (e - de + bf)/(1 - a - bc - d + ad)
Ymax = (ce + f - af)/(1 - a - bc - d + ad)

If I apply that to the non-slanted sierpinski triangle I get

Code:

Set 1: [0.5, 0.0, 0.0, 0.5, 0.0, 0.0]
Set 2: [0.5, 0.0, 0.0, 0.5, 1.0, 0.0]
Set 3: [0.5, 0.0, 0.0, 0.5, 0.5, 0.5]

Set 1: [Xmax = 0, Ymax = 0]
Set 2: [Xmax = 2, Ymax = 0]
Set 3: [Xmax = 1, Ymax = 1]

The highest maximums are Xmax = 2 and Ymax = 1 and for the slanted one I get the expected Xmax = 1.
So is this is a foolproof solution? I'm not sure how to get the minimum values though, and I'm not sure how to handle a possible division by zero.

I have already made changes so that the program uses probability. I also fixed the part that handles the boundschecking. It previously set the X and Y variables so they would be inside but that did of course cause those strange images I posted before. The ones I get now look more like the ones I should be getting, but at a much larger scale. I'll fix the code some more later if I can figure out the minimums.

Code:

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <mt19937ar.h>
#include <SDL/SDL.h>
#include <windows.h>

#define WIDTH  640
#define HEIGHT 480
#define ITER 90000
#define COLOR     0x00ff00
#define FILLCOLOR 0x000000
#define PITCH  (screen->pitch/4)
#define PIXELS ((unsigned int*) screen->pixels)

typedef struct {
	double a, b, c, d, e, f;
	double p;
} ifs_elem;

void draw_frac (ifs_elem* ifs);

SDL_Surface* screen;

int main (int argc, char** argv)
{
	SDL_Event event;
	ifs_elem sierp_slant[] = {{0.5, 0.0, 0.0, 0.5, 0.0*WIDTH, 0.0*HEIGHT, 1.0/3.0},
	                          {0.5, 0.0, 0.0, 0.5, 0.0*WIDTH, 0.5*HEIGHT, 1.0/3.0},
	                          {0.5, 0.0, 0.0, 0.5, 0.5*WIDTH, 0.5*HEIGHT, 1.0/3.0}};
	ifs_elem sierp[] = {{0.5, 0.0, 0.0, 0.5, 0.0*WIDTH/2, 0.0*HEIGHT, 1.0/3.0},
	                    {0.5, 0.0, 0.0, 0.5, 1.0*WIDTH/2, 0.0*HEIGHT, 1.0/3.0},
	                    {0.5, 0.0, 0.0, 0.5, 0.5*WIDTH/2, 0.5*HEIGHT, 1.0/3.0}};
	ifs_elem fern[] = {{ 0.0,   0.0,   0.0,  0.16, 0.0*WIDTH, 1.6*HEIGHT,  0.1},
	                   { 0.2,  -0.26,  0.23, 0.22, 0.0*WIDTH, 0.44*HEIGHT, 0.08},
	                   {-0.15,  0.28,  0.26, 0.24, 0.0*WIDTH, 1.6*HEIGHT,  0.08},
	                   { 0.75,  0.04, -0.04, 0.85, 0.0*WIDTH, 0.0*HEIGHT,  0.74}};
	ifs_elem fern2[] = {{ 0.0,   0.0,   0.0,  0.16, 0.0*WIDTH, 0.0*HEIGHT,  0.1},
	                    { 0.2,  -0.26,  0.23, 0.22, 0.0*WIDTH, 1.6*HEIGHT,  0.07},
	                    {-0.15,  0.28,  0.26, 0.24, 0.0*WIDTH, 0.44*HEIGHT, 0.07},
	                    { 0.85,  0.04, -0.04, 0.85, 0.0*WIDTH, 1.6*HEIGHT,  0.85}};
	ifs_elem* fracs[] = {sierp_slant, sierp, fern, fern2};
	int i = 0;

	if (SDL_Init(SDL_INIT_EVERYTHING) < 0) {
		fprintf(stderr, "Error: Unable to init SDL: %s\n", SDL_GetError());
		return 1;
	}
	atexit(SDL_Quit);
	if ((screen = SDL_SetVideoMode(WIDTH, HEIGHT, 32, SDL_HWSURFACE)) == NULL) {
		fprintf(stderr, "Error: Unable to init SDL: %s\n", SDL_GetError());
		return 1;
	}
	init_genrand(time(NULL));

	draw_frac(fracs[i]);
	while (SDL_WaitEvent(&event)) {
		if (event.type == SDL_KEYDOWN) {
			switch (event.key.keysym.sym) {
				case SDLK_q:
				case SDLK_ESCAPE:
					return 0;
	
				case SDLK_f:
					i++;
					if (i == sizeof(fracs)/sizeof(*fracs))
						i = 0;
					draw_frac(fracs[i]);
					break;
			}
		}
		else if (event.type == SDL_QUIT)
			return 0;
	}

	return 0;
}

void draw_frac (ifs_elem* ifs)
{
	int i, j;
	double x, y, oldx;
	double rnd;

	x = y = 0;
	SDL_FillRect(screen, NULL, FILLCOLOR);
	if (SDL_MUSTLOCK(screen))
		SDL_LockSurface(screen);
	for (i = 0; i < ITER; i++) {
		rnd = genrand_real1();
		for (j = 0; 1; j++) {
			if (ifs[j].p >= rnd)
				break;
			rnd -= ifs[j].p;
		}
		oldx = x;
		x = ifs[j].a*x    + ifs[j].b*y + ifs[j].e;
		y = ifs[j].c*oldx + ifs[j].d*y + ifs[j].f;
		if (x < 0 || x >= WIDTH || y > HEIGHT || y < 1)
			continue;
		PIXELS[(int) x + (HEIGHT - (int) y)*PITCH] = COLOR;
	}
	if (SDL_MUSTLOCK(screen))
		SDL_UnlockSurface(screen);
	while (SDL_Flip(screen) < 0);
}

**coder8137** · 01-14-2007

So is this is a foolproof solution?

Most probably not. It might work with simple equations, but not with the ones for the Fern (as far as I can see). For example, try it on the 2nd equation for the fern, given on the website. I'm getting Xmax = -0.61, yet, starting with (0,0) and iterating once immediately results in an X value of 0 (> Xmax). This may be because you assume:

Code:

Xmax = a*Xmax + b*Ymax + e
/* What if b is negative? Then it should be: */
Xmax = a*Xmax +b*Ymin +e
/* How about if a is negative... */
/* Also, would the equations, if applied properly, ever even result in a point
   (Xmax, Ymax).                                                                                                            */

I don't know if there is a simple equation to solve for xmin/xmax/ymin/ymax, but I have my doubts. Personally, I'd just iterate the equations a few million times and calculate the (approximate) values like that.

Code:

double xmin=0, xmax=0, ymin=0, ymax=0;
for(i=0; i<NOS_POINTS; i++) {
    /* calculate p[0], p[1] based on equations and probability */

    if(p[0] < xmin)
      xmin = p[0];
    else if(p[0] > xmax)
      xmax = p[0];

    if(p[1] < ymin)
      ymin = p[1];
    else if(p[1] > ymax)
      ymax = p[1];
}

Even, after you have obtained Xmin/max/Ymin/Ymax, you still have to modify the equations correctly, and this won't be as easy as it was for the triangle either. For a start, you might try:

Code:

/* matrixes[i][4] represents e.
              [i][5] represents f.       */

  for(i=0; i<nos_of_equations; i++) {
    matrixes[i][4] = matrixes[i][4] / (xmax - xmin) * ((double)WIDTH-1);
        /* modify X constant (e) - We haven't accounted for xmin. Everything is too
                                   far left                                         */
    matrixes[i][5] = matrixes[i][5] / (ymax - ymin) * ((double)HEIGHT-1);
        /* modify Y constant (f)  - We haven't accounted for ymin... */
  }

I can't get the right shift working atm, so, using the above, you will only see 3/4 of the fern.

If possible, in my opinion, you should just change your code to 2 Steps, to make life easier.

**OnionKnight** · 01-16-2007

If possible, in my opinion, you should just change your code to 2 Steps, to make life easier.

That's what I'm trying to do. As I understood it Step 1 is about finding the maximums and minimums so that the display can be adjusted and Step 2 is plotting everything.
Algebraically seemed to be a bit tougher than what I imagined so I'll go with the "brute force" method.

[EDIT] Strange typos.

**coder8137** · 01-16-2007

That's not quite what I meant. The problem with the way you're implementing it is that you can't change xmin/xmax/ymax/ymin freely. Most changes will result in the graph being completly wrong, rather than just the viewing area changing.

Here's an example of what I meant (xmin, xmax, ymin, ymax represent the plotting area and can be set to generic values):

Code:

/* WINDOW SPECS */
#define WIDTH  640
#define HEIGHT 480

void display(void)
{
  /* PLOTTING AREA - change at will */
  const double xmin  = -10.0;
  const double xmax = 10.0;
  const double ymin  = -10.0;
  const double ymax = 10.0;

  /* Any equations... left unchanged */
  ifs_elem ifs[] = {{0.0,   0.0,   0.0,  0.16, 0.0, 0.0,  0.1},
                       {0.2,  -0.26,  0.23, 0.22, 0.0, 1.6,  0.08},
                       {-0.15, 0.28,  0.26, 0.24, 0.0, 0.44, 0.08},
                       {0.75,  0.04, -0.04, 0.85, 0.0, 1.6,  0.74}};

  /* Vars for Step 1 - These are in our coordinate system */
  double x, y, oldx;

  /* Vars for Step 2 - These are in pixels */
  int p_x, p_y, p_y_top;

  /* Other Vars */
  int i;

  /* Start Coordinates - Probably could be any value */
  x = 0.0;
  y = 0.0;

  for (i = 0; i < ITER; i++) {
    /* STEP 1 - Calculate Points */
    /* --- snip ---> Calculating probability here. j indicates equation to use */
    oldx = x;
    x = ifs[j].a*x    + ifs[j].b*y + ifs[j].e;
    y = ifs[j].c*oldx + ifs[j].d*y + ifs[j].f;

    /* STEP 2 - Converting Points and Plotting */
    p_x = ((double)WIDTH - 1.0) / (xmax - xmin) * (-xmin + x) + 0.5;
    p_y = ((double)HEIGHT - 1.0) / (ymax - ymin) * (-ymin + y) + 0.5;
                       /* Note: 0.5 is just for rounding, not really necessary */
    p_y_top = HEIGHT - p_y -1; /* EDIT: Substracted 1 (from height) for consistancy) */
    if(p_y_top >= 0 && p_y_top < HEIGHT && p_x >= 0 && p_x < WIDTH)
      PIXELS[p_x + p_y_top*PITCH] = COLOR;
  }
}

Let me know if it works. Step 1 is just about setting x and y, which will be plotted depending on xmin/xmax/ymin/ymax. Step 2 should be generic for all plotting tasks (it could be optimized slightly if needed).

**OnionKnight** · 01-22-2007

So step 1 is just translation. I was thinking of something that would make sure that the each fractal would take up the whole screen but that has proven to be pretty complicated. I guess you could do one loop where you find the limits and then for the second loop where you plot out evertything you set the seed of the RNG to be the same as the one used in the first loop, that should make the boundaries will stay the same.

But I'm not all that interested in having it fit to the screen. I added some simple keys for moving around and zooming in on the fractal so I'm happy. Might add some more functionality if I get the time to toy around some more with this.

Thread: Iterated function system fractal

Thread Tools

Search Thread

Display

Iterated function system fractal

Similar Threads

Compiling sample DarkGDK Program

Screwy Linker Error - VC2005

measuring system resources used by a function

[question]Simulating the Reliability of a Component-Based System [with my code]

I need help with passing pointers in function calls