Microcontroller Abuse: using an AVR as a fractal generator

Discussion in 'General Software' started by LINUX, Oct 7, 2010.

  1. LINUX

    LINUX Member

    Joined:
    Sep 18, 2001
    Messages:
    3,062
    Location:
    Dudley, Newcastle
    I've always had a fascination of fractal generators. They've kind of been my version of "hello world!" when learning a new language.

    Recently I've gotten into programming AVRs and they have not been allowed to escape this historical trend. This is a Mandelbrot set rendered on an ATMEGA8L @ 8MHz and printed to Hyper Terminal:

    Code:
    ...............................................................................
    ...............................................................................
    ...............................................................................
    ...............................................................................
    ......................................................##.......................
    .....................................................####......................
    ....................................................#####......................
    ....................................................######.....................
    ....................................................#####......................
    ......................................................#.###....................
    ................................................##############.................
    ...........................................####################..###...........
    ............................................########################...........
    .......................................#..########################.............
    ........................................###########################............
    ........................................############################...........
    ......................................##############################.#.........
    ......................................###############################..........
    .........................#...#.#......###############################..........
    .........................########.....################################.........
    ........................###########..################################..........
    .......................#############.################################..........
    ......................##############.################################..........
    ...................#################################################...........
    ....#..#.......#.#################################################.............
    ...................#################################################...........
    ......................##############.################################..........
    .......................#############.################################..........
    ........................###########..################################..........
    .........................########.....################################.........
    .........................#...#.#......###############################..........
    ......................................###############################..........
    ......................................##############################.#.........
    ........................................############################...........
    ........................................###########################............
    .......................................#..########################.............
    ............................................########################...........
    ...........................................####################..###...........
    ................................................##############.................
    ......................................................#.###....................
    ....................................................#####......................
    ....................................................######.....................
    ....................................................#####......................
    .....................................................####......................
    ......................................................##.......................
    ...............................................................................
    ...............................................................................
    ...............................................................................
    Calculation time was "a few minutes" (wasn't timing).

    Using C made this fairly easy as all the floating point maths is done for you. If anyone's interested here's the code:

    It's written for gcc-avr.

    Code:
    #include <avr/io.h>
    #include <stdint.h>
    
    void usart_string(unsigned char *s);
    void usart_char(unsigned char a);
    void usart_number(unsigned int num);
    
    const unsigned char newline[] = {13, 10, 0};
    int isInside(float x, float y);
    
    int main(void)
    {
    	float x, y;
    	const float top = 1.3, bottom = -1.3, left = -2.0, right = 0.7;
    	float xinc, yinc;
    	const int xres = 78, yres=60;
    	unsigned char xpos, ypos;
    
    	//Setup UART
    	// Set frame format to 8 data bits, no parity, 1 stop bit
    	UCSRC = (0<<USBS)|(1<<UCSZ1)|(1<<UCSZ0);
    	// Enable transmitter
    	UCSRB = (1<<TXEN);
    	UCSRA = (1<<U2X);
    	UBRRH=0x00;
    	UBRRL=8;
    
    	xinc = (right-left)/(float)xres;
    	yinc = (top-bottom)/(float)yres;
    
    	for(ypos = 0, y=top; ypos <= yres; ypos++, y-=yinc)
    	{
    		for(xpos = 0,x=left; xpos<=xres; xpos++, x+=xinc)
    		{
    			if(isInside(x,y))
    				usart_char('#');
    			else
    				usart_char('.');
    		}
    		usart_string(newline);
    	}
    
    	return 1;
    }
    
    int isInside(float x, float y)
    {
    	const int maxIters=100;
    	int iters=0;
    	float xl=0.0, yl=0.0;
    	float xc=0.0, yc=0.0;
    
    	while(((xc*xc+yc*yc) < 4.0) && (iters < maxIters))
    	{
    		xc = xl*xl - yl*yl + x;
    		yc = 2*xl*yl + y;
    		xl=xc;
    		yl=yc;
    		iters++;
    	}
    	if(iters<maxIters)
    		return 0;
    	else
    		return 1;
    }
    
    void usart_string(unsigned char *s)
    {
    	unsigned char x=0;
    
    	while(s[x])
    	{
    		while ( (UCSRA&(1<<UDRE)) == 0 );
    		UDR=s[x];
    		x++;
    	}
    }
    
    void usart_char(unsigned char a)
    {
    	while ( (UCSRA&(1<<UDRE)) == 0 );
    	UDR=a;
    }
     
  2. scorpydude

    scorpydude Member

    Joined:
    Mar 13, 2008
    Messages:
    116
    I've always been curious about these...

    Whats the smallest amount of code to achieve the same in vb.net , so i can get my head around it?
     
  3. OP
    OP
    LINUX

    LINUX Member

    Joined:
    Sep 18, 2001
    Messages:
    3,062
    Location:
    Dudley, Newcastle
    Never used VB.net so I can only give you psudocode:

    Code:
    //Loop over a 2D area in the complex plain
    for y = 2 to -2 (top to bottom)
         for x = -2 to 2 (left to right)
              if the complex point x,y is inside the mandelbrot set draw a black pixel
              else draw a while pixel
         endfor
    endfor
    
    To decide if a point (x,y) is "inside" the mandelbrot set you do the following iteration:

    Code:
    z = z^2 + c
    
    Where z is a complex number x+iy, the formula in real and imaginary parts is:

    Code:
    iterations = 0
    loop
         x_new = x*x - y*y +cx
         y_new = 2*x*y + cy
         x = x_new
         y = y_new
         iterations++
    while (x*x + y*y < 4) AND iterations <= MAX_ITERATIONS//4 is known as the "bailout radius", it's a circle beyond which no point exists in the mandelbrot set.
    
    if(iterations == MAX_ITERATIONS)
      point is within set, return true
    else
      return false
    
    Where cx and cy are the real and imaginary parts of the constant, c. In the Mandelbrot algorithm the constant c is the point in question in the main loop a the start.

    That's a *very* rough overview. You need at least a basic understanding of complex arithmetic and mathematical iteration to programme this type of fractal generator. I'm happy to answer questions though.
     
    Last edited: Oct 8, 2010
  4. TheChemist

    TheChemist R.I.P

    Joined:
    Apr 5, 2007
    Messages:
    1,035
    VB.NET is an arsehole of a language and should never ever be used for anything serious.

    With that said,

    Code:
    Private Sub Form2_Paint(sender As Object, e As PaintEventArgs)
    	Dim s As Double = (4.0 / Me.Height), xl As Double, xc As Double, yl As Double, yc As Double
    	Dim x As Integer, y As Integer, i As Integer
    
    	For x = 0 To Me.Height - 1
    		For y = 0 To Me.Height - 1
    			i = 0
    			xl = 0F
    			yl = 0F
    			xc = 0F
    			yc = 0F
    
    			While ((xc * xc + yc * yc) < 4F) AndAlso (System.Threading.Interlocked.Increment(i) < 255)
    				xc = xl * xl - yl * yl + (x * s - 2F)
    				yc = 2 * xl * yl + (y * s - 2F)
    				xl = xc
    				yl = yc
    			End While
    
    			e.Graphics.FillRectangle(New SolidBrush(System.Drawing.Color.FromArgb(i, i, i)), x, y, 1, 1)
    		Next
    	Next
    End Sub
    
    This bumps it up to 255 iterations and is greyscale.
     
  5. OmNomNomNom

    OmNomNomNom (Banned or Deleted)

    Joined:
    Oct 13, 2010
    Messages:
    11
    Location:
    Australia
    That's pretty cool, but I'm surprised you're still using a tiny ATmega8 and not one of the replacements with more memory such as the ATmega328.
     
  6. OP
    OP
    LINUX

    LINUX Member

    Joined:
    Sep 18, 2001
    Messages:
    3,062
    Location:
    Dudley, Newcastle
    Well, there are 2 reasons: 1) I wasn't aware of the 328 because the Jaycar catalog is my marketing document and 2) Farnell sell the for 8 half the price of the 328. I grabbed a 32A in the same order for playing around with more memory :).

    Edit: The 328 runs up to 20MHz? Hrm...might grab one in the next order, thanks for the heads up about it.
     
    Last edited: Oct 14, 2010
  7. OP
    OP
    LINUX

    LINUX Member

    Joined:
    Sep 18, 2001
    Messages:
    3,062
    Location:
    Dudley, Newcastle
    Tonight I played with this a bit more and modified the Mandelbrot iteration loop to use fixed point integer arithmatic as opposed to floating point. This reduced the calculation time on a 200x60 pixel 100 iteration fractal from 6min 52s to 1min 8s :). System clock was set to 8MHz in both cases (Max for the atmega8L I'm using).

    There are some minor differences in the output (I can post both if people are really interested) but they're basically the same. I haven't optimised the fixed point calculation for max precision before overflow.

    I'll just post the iteration loop, the rest of the code is identical

    Code:
    int isInside(float x, float y)
    {
    	const int maxIters=100;
    	int iters=0;
    	long xl=0, yl=0;
    	long xc=0, yc=0;
    	long xi;
    	long yi;
    	const long div = 10000;
    
    	xi = (long)(div*x);
    	yi = (long)(div*y);
    
    	while( ((((xc*xc)+(yc*yc))) < 4*div*div) && (iters < maxIters))
    	{
    		xc = ((xl*xl) - (yl*yl))/div + xi;
    		yc = (2*xl*yl)/div + yi;
    		xl=xc;
    		yl=yc;
    		iters++;
    	}
    	if(iters<maxIters)
    		return 1;
    	else
    		return 0;
    }
    
    The value of div defines the precision, which in this case is 4 decimal places (cf ~7dp for 32-bit float). I haven't tried increasing this value but I suspect it's optimised for a 32-bit long.

    In order to prevent overflow in the multiplicaion the maximum value of the integer needs to be less than the square root of of the maximum value of the data type. For example if a datatype can hold values as high as 10 the fixed point representation musn't exceed sqrt(10) = 3.2 = 3 in order for a square multiplication to not overflow.

    In the case of 32-bit longs this value is 46 340 (sqrt ( 2^32 / 2 ), div by 2 because it's signed), or 4.6340 in the above fixed point calculations.

    On an x86 machine I played with long long's and found that 9 decimal places was the best that could be achieved in the fractal calculation before overflows destryed it.

    Anway, not a bad little exercise :). Much was learned. The code could be optimised further if the main() loop was all converted to fixed point. This wouldn't result in much of a speed increase as the iteration loop is called *far* more often that the point incrementation but it would cause the code size to shrink dramatically as none of the floating point emulation code would need to be compiled in.
     
  8. Impulse IX

    Impulse IX Member

    Joined:
    Jul 31, 2001
    Messages:
    437
    Location:
    Adelaide/Brisbane
    I ported your code from your first post to my Arduino Duemilanove. 8 Second render.

    Code:
    void setup()
    {
      Serial.begin(9600);
      // 5 Second delay to open serial monitor
      delay(5000);
    }
    
    void loop()
    {
      float x, y;
      const float top = 1.3, bottom = -1.3, left = -2.0, right = 0.7;
      float xinc, yinc;
      const int xres = 78, yres=60;
      unsigned char xpos, ypos;
    
      xinc = (right-left)/(float)xres;
      yinc = (top-bottom)/(float)yres;
    
      for(ypos = 0, y=top; ypos <= yres; ypos++, y-=yinc)
      {
        for(xpos = 0,x=left; xpos<=xres; xpos++, x+=xinc)
        {
          if(isInside(x,y))
            Serial.print('#');
          else
            Serial.print('.');
        }
        Serial.print(10,BYTE);
      }
    
    }
    
    int isInside(float x, float y)
    {
      const int maxIters=100;
      int iters=0;
      float xl=0.0, yl=0.0;
      float xc=0.0, yc=0.0;
    
      while(((xc*xc+yc*yc) < 4.0) && (iters < maxIters))
      {
        xc = xl*xl - yl*yl + x;
        yc = 2*xl*yl + y;
        xl=xc;
        yl=yc;
        iters++;
      }
      if(iters<maxIters)
        return 0;
      else
        return 1;
    }
    
     
  9. OP
    OP
    LINUX

    LINUX Member

    Joined:
    Sep 18, 2001
    Messages:
    3,062
    Location:
    Dudley, Newcastle
    8 Seconds? Cripes! Are you able to get the ASM out of the compiler?

    Edit: And do you know how many bits the default float datatype uses? The Atmega168 on that board can run at 20MHz, but even ~20s (if it's at 8MHz) is nothing compared to the several minutes that the avr-gcc code takes.
     
    Last edited: Oct 14, 2010

Share This Page