Simple example program for the PIC12F675 microcontroller

The PIC12F675 is a very inexpensive 8-bit microcontroller from Microchip that’s available in an 8-pin DIP package (i.e. it’s suitable for breadboard use). I’ve had a small tube of these on the shelf for a few years, but I’ve only just got around to trying one out. In this post, I present a simple example circuit and program, which flashes an LED on one of the GPIO pins. The example circuit includes connections to the PICkit2 USB programmer that I’m using to program the device. The PICkit2 also supplies power to the circuit.

First, here’s the video of the example program running on the PIC12F675.

Example Code

This is the full C code for the flashing LED example. I built it using Microchip’s XC8 compiler (the compiler command line is shown in the opening comments).

//
// PIC12F675 example: blink an LED on pin GP5
// Written by Ted Burke - 18-2-2017
//
// To compile:
//
//    xc8 --chip=12F675 main.c
//

#include <xc.h>

#pragma config FOSC=INTRCIO,WDTE=OFF,MCLRE=OFF,BOREN=OFF

void main(void)
{
    TRISIO = 0b11011111; // Make pin GP5 a digital output
    
    while(1)
    {
        GP5 = 1;         // Set pin GP5 high
        _delay(500000);  // 0.5 second delay
        GP5 = 0;         // Set pin GP5 low
        _delay(500000);  // 0.5 second delay
    }
}

Example Circuit

pic12f675_circuit

Click here to download editable Inkscape SVG file of circuit diagram.

Here’s my breadboard circuit:

20170217_162941

Posted in Uncategorized | Tagged , , , , , , , , , , , | Leave a comment

A simple way to read and write audio and video files in C using FFmpeg (part 2: video)

In my previous post, I demonstrated how FFmpeg can be used to pipe raw audio samples in and out of a simple C program to or from media files such as WAV files (based on something similar for Python I found on Zulko’s blog). The same idea can be used to perform video processing, as shown in the program below.

Reading and writing a video file

In this example I use two pipes, each connected to its own instance of FFmpeg. Basically, I read frames one at time from the input pipe, invert the colour of every pixel, and then write the modified frames to the output pipe. The input video I’m using is teapot.mp4, which I recorded on my phone. The modified video is saved to a second file, output.mp4. The video resolution is 1280×720, which I checked in advance using the ffprobe utility that comes with FFmpeg.

The full code is below, but first let’s see the original and modified videos:

The original and modified MP4 video files can be downloaded here:

The program I wrote to convert the original video into the modified version is shown below.

//
// Video processing example using FFmpeg
// Written by Ted Burke - last updated 12-2-2017
//

#include <stdio.h>

// Video resolution
#define W 1280
#define H 720

// Allocate a buffer to store one frame
unsigned char frame[H][W][3] = {0};

void main()
{
    int x, y, count;
    
    // Open an input pipe from ffmpeg and an output pipe to a second instance of ffmpeg
    FILE *pipein = popen("ffmpeg -i teapot.mp4 -f image2pipe -vcodec rawvideo -pix_fmt rgb24 -", "r");
    FILE *pipeout = popen("ffmpeg -y -f rawvideo -vcodec rawvideo -pix_fmt rgb24 -s 1280x720 -r 25 -i - -f mp4 -q:v 5 -an -vcodec mpeg4 output.mp4", "w");
    
    // Process video frames
    while(1)
    {
        // Read a frame from the input pipe into the buffer
        count = fread(frame, 1, H*W*3, pipein);
        
        // If we didn't get a frame of video, we're probably at the end
        if (count != H*W*3) break;
        
        // Process this frame
        for (y=0 ; y<H ; ++y) for (x=0 ; x<W ; ++x)
        {
            // Invert each colour component in every pixel
            frame[y][x][0] = 255 - frame[y][x][0]; // red
            frame[y][x][1] = 255 - frame[y][x][1]; // green
            frame[y][x][2] = 255 - frame[y][x][2]; // blue
        }
        
        // Write this frame to the output pipe
        fwrite(frame, 1, H*W*3, pipeout);
    }
    
    // Flush and close input and output pipes
    fflush(pipein);
    pclose(pipein);
    fflush(pipeout);
    pclose(pipeout);
}

The FFmpeg options used for the input pipe are as follows.

FFmpeg option Explanation
-i teapot.mp4 Selects teapot.mp4 as the input file.
-f image2pipe Tells FFmpeg to convert the video into a sequence of frame images (I think!).
-vcodec rawvideo Tells FFmpeg to output raw video data (i.e. plain unencoded pixels).
-pix_fmt rgb24 Sets the pixel format of the raw data produced by FFmpeg to 3-bytes per pixel – one byte for red, one for blue and one for green.
- This final “-” tells FFmpeg to write to stdout, which in this case will send it into our C program via the input pipe.

The FFmpeg options used for the output pipe are as follows.

FFmpeg option Explanation
-y Tells FFmpeg to overwrite the output file if it already exists.
-f rawvideo Sets the input format as raw video data. I’m not too sure about the relationship between this option and the next one!
-vcodec rawvideo Tells FFmpeg to interpret its input as raw video data (i.e. unencoded frames of plain pixels).
-pix_fmt rgb24 Sets the input pixel format to 3-byte RGB pixels – one byte for red, one for blue and one for green.
-s 1280x720 Sets the frame size to 1280×720 pixels. FFmpeg will form the incoming pixel data into frames of this size.
-r 25 Sets the frame rate of the incoming data to 25 frames per second.
-i - Tells FFmpeg to read its input from stdin, which means it will be reading the data out C program writes to its output pipe.
-f mp4 Sets the output file format to MP4.
-q:v 5 This controls the quality of the encoded MP4 file. The numerical range for this option is from 1 (highest quality, biggest file size) to 32 (lowest quality, smallest file size). I arrived at a value of 5 by trial and error. Subjectively, it seemed to me to give roughly the best trade off between file size and quality.
-an Specifies no audio stream in the output file.
-vcodec mpeg4 Tells FFmpeg to use its “mpeg4” encoder. I didn’t try any others.
output.mp4 Specifies output.mp4 as the output file.

Epilogue: concatenating videos and images

As an interesting aside, the video I embedded above showing the original and modified teapot videos was spliced together using a slightly modified version of the example program above. Of course, it’s possible to concatenate (splice together) multiple files of different formats using ffmpeg on its own, but I couldn’t quite figure out the correct command line, so I just wrote my own little program to do it.

The files being joined together are:

The full source code is shown below.

//
// combine.c - Join multiple MP4 videos and PNG images into one video
// Written by Ted Burke - last updated 12-2-2017
//
// To compile:
// 
//    gcc combine.c
// 

#include <stdio.h>

// Video resolution
#define W 1280
#define H 720

// Allocate a buffer to store one frame
unsigned char frame[H][W][3] = {0};

void main()
{
    int count, n;
    FILE *pipein;
    FILE *pipeout;
    
    // Open output pipe
    pipeout = popen("ffmpeg -y -f rawvideo -vcodec rawvideo -pix_fmt rgb24 -s 1280x720 -r 25 -i - -f mp4 -q:v 5 -an -vcodec mpeg4 combined.mp4", "w");
    
    // Write first 50 frames using original video title image from title_original.png
    pipein = popen("ffmpeg -i title_original.png -f image2pipe -vcodec rawvideo -pix_fmt rgb24 -", "r");
    count = fread(frame, 1, H*W*3, pipein);
    for (n=0 ; n<50 ; ++n)
    {
        fwrite(frame, 1, H*W*3, pipeout);
        fflush(pipeout);
    }
    fflush(pipein);
    pclose(pipein);
    
    // Copy all frames from teapot.mp4 to output pipe
    pipein = popen("ffmpeg -i teapot.mp4 -f image2pipe -vcodec rawvideo -pix_fmt rgb24 -", "r");
    while(1)
    {
        count = fread(frame, 1, H*W*3, pipein);
        if (count != H*W*3) break;
        fwrite(frame, 1, H*W*3, pipeout);
        fflush(pipeout);
    }
    fflush(pipein);
    pclose(pipein);

    // Write next 50 frames using modified video title image from title_modified.png
    pipein = popen("ffmpeg -i title_modified.png -f image2pipe -vcodec rawvideo -pix_fmt rgb24 -", "r");
    count = fread(frame, 1, H*W*3, pipein);
    for (n=0 ; n<50 ; ++n)
    {
        fwrite(frame, 1, H*W*3, pipeout);
        fflush(pipeout);
    }
    fflush(pipein);
    pclose(pipein);
    
    // Copy all frames from output.mp4 to output pipe
    pipein = popen("ffmpeg -i output.mp4 -f image2pipe -vcodec rawvideo -pix_fmt rgb24 -", "r");
    while(1)
    {
        count = fread(frame, 1, H*W*3, pipein);
        if (count != H*W*3) break;
        fwrite(frame, 1, H*W*3, pipeout);
        fflush(pipeout);
    }
    fflush(pipein);
    pclose(pipein);
    
    // Flush and close output pipe
    fflush(pipeout);
    pclose(pipeout);
}

Note: I used Inkscape to create the video title images. Click here to download the editable Inkscape SVG file.

Posted in Uncategorized | Tagged , , , , , , , , , , , , , , , , , , , , | Leave a comment

A simple way to read and write audio and video files in C using FFmpeg (part 1: audio)

C is my favourite programming language and the one I use most often. However, I have tended to shy away from using it for quick one-off signal processing tasks where I needed to read or write audio or video files because it always seemed like a lot of hassle figuring out which library to use and what functions to call to actually get access to the raw data. Recently though, I’ve found a new way of dealing with audio and video files that is relatively painless, so I think it’s worth sharing.

I first learned of this approach from an article on Zulko’s fascinating Python oriented blog on GitHub:

The article describes using FFmpeg to read and write video frames in Python and there’s a link to a second article showing the same thing for audio.

The basic idea is to launch FFmpeg via a pipe, which then converts raw samples to the required format (for writing) or decodes the file into raw samples (for reading). It may seem like a bit of a hack, but it’s surprisingly effective and, provided that you can figure out the correct FFmpeg command line, extremely adaptable.

FFmpeg is described on its project web page as “[a] complete, cross-platform solution to record, convert and stream audio and video.” That’s a pretty accurate description. So far, I haven’t needed to read or write any file format that it couldn’t handle, although in some cases it did involve some googling to figure out the required command line arguments.

FFmpeg is free software and cross-platform, so you should be able to install it on Windows and Mac without too much difficulty, but I haven’t actually tried that myself. I’m working on Xubuntu Linux, so to install FFmpeg I just did this:

sudo apt-get install ffmpeg

Ok, time for some practical examples…

Writing a WAV audio file (or MP3, FLAC, or whatever format you like)

The following example creates a 1-second long WAV audio file (PCM, 16-bit signed integer samples, 44.1 kHz sampling frequency). The sound is just a loud 1 kHz sine wave.

//
// writewav.c - Create a wav file by piping raw samples to ffmpeg
// Written by Ted Burke - last updated 10-2-2017
//
// To compile:
//
//    gcc writewav.c -o writewav -lm
//

#include <stdio.h>
#include <stdint.h>
#include <math.h>

#define N 44100

void main()
{
    // Create audio buffer
    int16_t buf[N] = {0}; // buffer
    int n;                // buffer index
    double Fs = 44100.0;  // sampling frequency
    
    // Generate 1 second of audio data - it's just a 1 kHz sine wave
    for (n=0 ; n<N ; ++n) buf[n] = 16383.0 * sin(n*1000.0*2.0*M_PI/Fs);
    
    // Pipe the audio data to ffmpeg, which writes it to a wav file
    FILE *pipeout;
    pipeout = popen("ffmpeg -y -f s16le -ar 44100 -ac 1 -i - beep.wav", "w");
    fwrite(buf, 2, N, pipeout);
    pclose(pipeout);
}

General comments:

  • I used the fixed length int data type int16_t for the audio buffer just to make absolutely sure I got 16-bit signed integer samples. That type is defined in the stdint.h header file. You can use a different data type, provided it matches the format your writing.
  • In the single-line for loop where I generate the sine wave, I used M_PI which is actually just the value of π (the mathematical constant). M_PI is defined in the math.h header file. Because I used M_PI, I had to link to the math library, which is why the suggested compiler command line includes the -lm option.
  • The popen() function launches a separate program which is accessible via a “pipe” using C’s standard file i/o functions. The first argument to popen is the command line for the program being launched (in this case FFmpeg). The second argument to popen specifies whether we’ll be reading data from that program through the pipe or vice versa. In this case, the second argument is "w" which means that we’ll be sending data to FFmpeg through the pipe. The return value from popen is a file descriptor for the pipe.
  • In this example, the fwrite() function is used to send the entire contents of the audio buffer into the pipe.
  • Finally, the pclose() function is used to close the pipe.

The FFmpeg command line:

The first argument to the popen() function is a string containing the full FFmpeg command line. The specified options included in that command line control how FFmpeg will interpret the raw sample data it receives through the pipe and the type of file it will write.

  • The "-y" option tells FFmpeg that it can overwrite the specified output file if it already exists.
  • The "-f s16le" option tells FFmpeg that the format of the audio data it reads (from its standard input, which means via the pipe from our program) is raw PCM, signed integer, 16-bit and little-endian.
  • The "-ar 44100" option tells FFmpeg that the sampling rate (i.e. sampling frequency) of the audio data it reads is 44.1 kHz.
  • The "-ac 1" option tells FFmpeg that the number of channels in the audio data it reads is 1.
  • The "-i -" option tells FFmpeg to read its input from standard input, which in this case means from the pipe, since FFmpeg was launched by popen.
  • Finally, beep.wav is the output filename FFmpeg will use.

If you want to hear the output file, you can download it here: beep.wav

Reading a WAV audio file (or MP3, FLAC or whatever format you like)

The following example reads 20 ms of samples from the beginning of a WAV file called whistle.wav (click it to download the file), then prints the sample values to a CSV file. The WAV file format is 16-bit signed integer samples, mono, with a sampling frequency of 44.1 kHz.

//
// readwav.c - Read samples from a WAV file using FFmpeg via a pipe
// Written by Ted Burke - last updated 10-2-2017
//
// To compile:
//
//    gcc readwav.c -o readwav -lm
//

#include <stdio.h>
#include <stdint.h>

#define N 882

void main()
{
    // Create a 20 ms audio buffer (assuming Fs = 44.1 kHz)
    int16_t buf[N] = {0}; // buffer
    int n;                // buffer index
    
    // Open WAV file with FFmpeg and read raw samples via the pipe.
    FILE *pipein;
    pipein = popen("ffmpeg -i whistle.wav -f s16le -ac 1 -", "r");
    fread(buf, 2, N, pipein);
    pclose(pipein);
    
    // Print the sample values in the buffer to a CSV file
    FILE *csvfile;
    csvfile = fopen("samples.csv", "w");
    for (n=0 ; n<N ; ++n) fprintf(csvfile, "%d\n", buf[n]);
    fclose(csvfile);
}

A screenshot of the resulting CSV file, samples.csv, open for viewing in my text editor is shown below.

samples-csv_geany_screenshot

This example program shares many elements with the previous one, but there are a couple of noteworthy differences:

  • The FFmpeg command line specified in the first argument of the popen() function tells FFmpeg to read its input from the file whistle.wav (the exact format will be detected automatically) and write its output to stdout as raw samples (16-bit signed integers, little-endian).
  • The second argument to popen is "r", which means that our program will read FFmpeg’s standard output via the pipe (it was the other way around in the previous example).

To ensure that the raw samples coming from FFmpeg are the correct size for my audio buffer, I specified "-f s16le" as a command line argument to FFmpeg. This was actually redundant in my case, because the WAV file I happened to be using was already in that format. I decided to explicitly specify the format anyway, just to show that FFmpeg can convert to a specific format if desired.

Modifying a WAV audio file

This example reads, modifies and writes a WAV audio file. It actually launches two instances of FFmpeg. One opens the original WAV file ("12345678.wav") and passes raw 16-bit signed integer samples into this program via a pipe. The other receives modified samples from this program via a second pipe and writes them to another WAV file ("out.wav"). The modification that is applied to the samples is a simple tremolo effect (modulation of the signal amplitude).

//
// modifywav.c - Modify a WAV file using FFmpeg via pipes
// This example adds a crude tremolo effect to the audio
// Written by Ted Burke - last updated 10-2-2017
//
// To compile:
//
//    gcc modifywav.c -o modifywav -lm
//

#include <stdio.h>
#include <stdint.h>
#include <math.h>

void main()
{
    // Launch two instances of FFmpeg, one to read the original WAV
    // file and another to write the modified WAV file. In each case,
    // data passes between this program and FFmpeg through a pipe.
    FILE *pipein;
    FILE *pipeout;
    pipein  = popen("ffmpeg -i 12345678.wav -f s16le -ac 1 -", "r");
    pipeout = popen("ffmpeg -y -f s16le -ar 44100 -ac 1 -i - out.wav", "w");
    
    // Read, modify and write one sample at a time
    int16_t sample;
    int count, n=0;
    while(1)
    {
        count = fread(&sample, 2, 1, pipein); // read one 2-byte sample
        if (count != 1) break;
        ++n;
        sample = sample * sin(n * 5.0 * 2*M_PI / 44100.0);
        fwrite(&sample, 2, 1, pipeout);
    }
    
    // Close input and output pipes
    pclose(pipein);    
    pclose(pipeout);
}

In case you want to hear the original and modified WAV audio files, here they are:

Click here to continue to part 2, where I show how to apply the same approach to video processing.

Posted in Uncategorized | Tagged , , , , , , , , , , , , , , , , , , , , , | Leave a comment

Fraktalismus Outtakes: Sea Monsters

Sea Creature 1

seacreature1

//
// seacreature1.c - Creates a fractal "sea creature"
// Written by Ted Burke, last updated 3-8-2016
//
// gcc seacreature1.c -o seacreature1 -lm
// ./seacreature1
// ffmpeg -framerate 20 -i %03d.pgm seacreature1.gif
// 

#include <stdio.h>
#include <complex.h>
#include <math.h>

#define W 640
#define H 640

void main()
{
    unsigned char pixels[H][W];
    
    // n is used to count iterations of the iterating function
    // x and y and pixel coordinates
    // t is frame number and T is the total number of frames
    // z and c are the complex variables in the iterating function
    // pwx was originally pixel width, but isn't exactly that now
    int n, x, y, t, T=100;
    complex double z, c;
    double pxw = 4.0 / H;
    
    char filename[100];
    FILE *f;
    
    // Generate all frames in the animation
    for (t=0 ; t<T ; ++t)
    {
        // Calculate pixel values for current frame
        for (y=0 ; y<H ; ++y) for (x=0 ; x<W ; ++x)
        {
            c = cexp(I*(t*M_PI*2.0/T)) * cpow(pxw * ((x-W/2) + (y-H/2)*I), -3.0 + 0.0*creal(cexp(I*(2.0*M_PI*t)/(0.5*T))));
            z = -0.25 + 0.25*cexp(I*(2.0*M_PI*t)/(0.5*T));
            for (n=0 ; n<50 && cabs(z)<4.0 ; ++n) z = z*z + c;
            pixels[y][x] = 5*n;
        }
        
        // Write current frame to PGM image file
        sprintf(filename, "%03d.pgm", t);
        fprintf(stderr, "Writing %s...", filename);
        f = fopen(filename, "w");
        fprintf(f, "P5\n%d %d\n255\n", W, H);
        fwrite(pixels, 1, W*H, f);
        fclose(f);
        fprintf(stderr, "DONE\n");
    }
}

Sea Creature 2

seacreature2

//
// seacreature2.c - Creates a fractal "sea creature"
// Written by Ted Burke, last updated 3-8-2016
//
// gcc seacreature2.c -o seacreature2 -lm
// ./seacreature2
// ffmpeg -framerate 30 -i %03d.pgm -vf reverse seacreature2.gif
// 

#include <stdio.h>
#include <complex.h>
#include <math.h>

#define W 640
#define H 640

void main()
{
    unsigned char pixels[H][W];
    
    // n is used to count iterations of the iterating function
    // x and y and pixel coordinates
    // t is frame number and T is the total number of frames
    // z and c are the complex variables in the iterating function
    // pwx was originally pixel width, but isn't exactly that now
    int n, x, y, t, T=100;
    complex double z, c, centre;
    double pxw = 6.0 / H;
    
    char filename[100];
    FILE *f;
    
    centre = -0.625 -0.1*I;

    // Generate all frames in the animation
    for (t=0 ; t<T ; ++t)
    {
        // Calculate pixel values for current frame
        for (y=0 ; y<H ; ++y) for (x=0 ; x<W ; ++x)
        {
            c = 1/(pxw * ((x-W/3) + (y-H/2)*I));
            z = centre + 0.175 + 0.075*cexp(I*(2.0*M_PI*t)/T);
            for (n=0 ; n<50 && cabs(z)<8.0 ; ++n) z = z*z + c;
            pixels[y][x] = n<2 ? 255 : 5*n;
        }
        
        // Write current frame to PGM image file
        sprintf(filename, "%03d.pgm", t);
        fprintf(stderr, "Writing %s...", filename);
        f = fopen(filename, "w");
        fprintf(f, "P5\n%d %d\n255\n", W, H);
        fwrite(pixels, 1, W*H, f);
        fclose(f);
        fprintf(stderr, "DONE\n");
    }
}
Posted in Uncategorized | Tagged , , , , , , , , | 1 Comment

Fraktalismus – my presentation at Dublin Maker 2016

A zip archive of my complete Fraktalismus folder (containing all the content for the presentation) can be download here:

001_title

002_fractal

003_complex_numbers

004_drawing_fractals

005_iterating_functions

006_lets_see_it

fractal_from_scratch_screenshot

fractal

This is the code I used to create the Julia Set fractal zoom animation above:

//
// zoom.c - Julia Set fractal zoom animation
// Written by Ted Burke - last modified 23-7-2016
//
// To compile and run in Linux:
// 
//    gcc zoom.c -o zoom -lm
//    ./zoom
//
// To combine the resulting 100 PNG frames into an mp4 video:
// 
//    ffmpeg -framerate 10 -i %03d.png zoomin.mp4
//    ffmpeg -i zoomin.mp4 -vf reverse zoomout.mp4
//    ffmpeg -f concat -i list.txt -c copy zoom.mp4
// 
// where list.txt contains this...
// 
//    file 'zoomin.mp4'
//    file 'zoomout.mp4'
//    file 'zoomin.mp4'
//    file 'zoomout.mp4'
//

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

#define W 1024
#define H 768

void main()
{
    complex float z, c;
    int i, y, x, n;
    unsigned char pixels[H][W];
    float px=0.003;
    
    FILE *f;
    char filename[100];
    char command[255];
    
    for (i=0 ; i<100 ; ++i)
    {
        px = 0.95 * px;
        
        for (y=0 ; y<H ; ++y)
        {
            for (x=0 ; x<W ; ++x)
            {
                c = -0.625 -0.4*I;
                z = px*((x-W/2) - I*(y-H/2));
                z = z - (-0.03632573427952118+0.027244300709640887*I);
                for (n=0 ; cabs(z)<2 && n<255 ; ++n) z = z*z + c;
                pixels[y][x] = n;
            }
        }
        
        sprintf(filename, "%03d.pgm", i);
        f = fopen(filename, "w");
        fprintf(f, "P5\n%d %d\n255\n", W, H);
        fwrite(pixels, W, H, f);
        fclose(f);
        
        sprintf(command, "mogrify -format png %s", filename);
        system(command);
        sprintf(command, "rm %s", filename);
        system(command);
    }
}

007_pantograph

Pantograph animation is by AlphaZeta, sourced from Wikimedia Commons

fractal_explorer_screenshot

This is the complete Python code for the Julia Set Explorer:

#
# explorer.py - Real time orbit viewer for iterative complex functions
# Written by Ted Burke, last updated 21-1-2016
#
 
from Tkinter import *
import math
import sys
 
def generate_fractal():
    sys.stdout.write('Generating fractal for c = {:.3f} + {:.3f}j...'.format(c.real,c.imag))
    sys.stdout.flush()
    for y in range(h):
        for x in range(w):
            z = centre + complex(((x-(w-1)/2.0)*pxw),(((h-1)/2.0-y)*pxw))
            n = 0
            while abs(z) < limit and n < 51:
                try:
                    z = pow(z,a) + c
                    n = n + 1
                except ZeroDivisionError:
                    z = limit
            pixel = int(255 * (0.5 + 0.5*math.cos(math.pi*n/51.0)))
            img.put('#{:02x}{:02x}{:02x}'.format(pixel,pixel,pixel),(x,y))
    sys.stdout.write('OK\n')
    sys.stdout.flush()
 
def set_z0(event):
    global z0
    z0 = complex((event.x - w/2) * pxw, (h/2 - event.y) * pxw)
    paint()
 
def set_c(event):
    global c
    c = complex((event.x - w/2) * pxw, (h/2 - event.y) * pxw)
    generate_fractal()
    paint()
 
def paint():
    canv.delete("all")
    canv.create_rectangle(0, 0, w, h, fill="white")
    canv.create_image((0,0), anchor="nw", image=img, state="normal")
    z = z0
    x = int(w/2 + z.real/pxw)
    y = int(h/2 - z.imag/pxw)
    canv.create_oval(x-3,y-3,x+4,y+4,fill="green")
    for n in range(100):
        p1 = z
        for m in range(M+2):
            if m<=M:
                p2 = pow(z,1.0+((a-1.0)*m)/M)
            else:
                p2 = p2 + c
            x1 = int(w/2 + p1.real / pxw)
            y1 = int(h/2 - p1.imag / pxw)
            x2 = int(w/2 + p2.real / pxw)
            y2 = int(h/2 - p2.imag / pxw)
            canv.create_line(x1, y1, x2, y2, fill="blue")
            p1 = p2
        z = pow(z,a) + c
        x = int(w/2 + z.real/pxw)
        y = int(h/2 - z.imag/pxw)
        canv.create_oval(x-3,y-3,x+4,y+4,fill="red")
        if (abs(z) > 10): break
 
# Create master Tk widget
master = Tk()
master.title('Julia Set Explorer - by Ted Burke - see http://batchloaf.com')
 
# Dimensions and resolution
w,h = 800,800
centre = 0 + 0j
pxw = 0.005
 
# Create Tk canvas widget
canv = Canvas(master, width=w, height=h)
canv.pack()
canv.bind("<B1-Motion>", set_z0)
canv.bind("<Button-1>",  set_z0)
canv.bind("<Button-2>",  set_c)
canv.bind("<Button-3>",  set_c) # this is actually the right click
 
# Iterating function parameters
a = 2.0
c = 0.325 + 0.285j
z0 = 0 + 0j
limit = 10
 
# Number of steps in plotting each curve within each orbit
M = 40
 
# Create image object to store Julia Set image
img = PhotoImage(width=w, height=h)
 
# Generate initial Julia Set image
generate_fractal()
paint()
 
# Main Tk event loop
mainloop()

The code used to create the Julia Set rotation animation above:

//
// rotating.c - Julia set animation generator
// Written by Ted Burke - last modified 23-7-2016
// 
// To compile and run (in Linux):
//
//    gcc rotating.c -o rotating -lm
//    ./rotating
// 
// To convert the resulting PGM files into PNG:
//
//    mogrify -format png *.pgm
// 
// To combine the PNG frames into an mp4 video:
//
//    ffmpeg -i %03d.png rotating.mp4
//

#include <stdio.h>
#include <complex.h>
#include <math.h>

void main()
{
    complex double z, c, c0;
    int x, y, n, N, w, h;
    unsigned char i;
    double r, px;
    
    px = 0.0025;
    w = 1024;
    h = 768;
    c0 = -0.625 -0.4*I;
    N = 200;
    r = 0.1;
    
    FILE *f;
    char filename[100];
    
    for (n=0 ; n<N ; ++n)
    {
        sprintf(filename, "%03d.pgm", n);
        f = fopen(filename, "w");
        printf("Generating %s\n", filename);
        fprintf(f, "P5\n%d %d\n255\n", w, h);
                
        for (y=0 ; y<h ; ++y)
        {
            for (x=0 ; x<w ; ++x)
            {
                z = px*(x-w/2) + px*(y-h/2)*I;
                c = c0 + r * z * cexp(I*(2*M_PI*n)/N);
                i = 0;
                while(cabs(z)<2 && ++i<127) z = z*z + c;
                fputc(255-2*i, f);
            }
        }
        
        fclose(f);
    }
}

009_slide6

The code used to generate the lifeform fractal animation above:

//
// lifeform.c - Julia Set animation generator
// Written by Ted Burke, last updated 23-7-2016
//
// To compile and run (in Linux):
//
//    gcc lifeform.c -o lifeform -lm
//    ./lifeform
// 
// To combine the resulting PNG frames into an mp4 video:
//
//    ffmpeg -i %03d.png lifeform.mp4
//

#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <stdlib.h>

#define W 1024
#define H 768

void main()
{
    complex float centre = 0 + 2.4I; // complex value at centre of image
    float pxw = 0.0053;              // pixel width
    
    float limit = 4.0;             // once z reaches this value, iteration ceases
    float a, b;
    complex float c, z;
    int t, x, y, n;
    char filename[100];
    FILE *f;

    // Create a buffer to store a row of pixel values
    unsigned char pixels[H][W] = {0};
     
    // Create fractal image files
    for (t=0 ; t<100 ; ++t)
    {
        a = 5.0 + 0.05*cos(M_PI*t/50.0);
        b = 2.5 + 0.25*sin(M_PI*t/50.0);
         
        // Generate image
        for (y=0 ; y<H ; ++y)
        {
            for (x=0 ; x<W ; ++x)
            {
                c = centre + (x-(W-1)/2.0)*pxw + ((H-1)/2.0-y)*pxw*I;
                z = 0.001; // Using exactly zero seems to cause star speckles
                n = 0;
                while (cabs(z) < limit && n < 25)
                {
                    if (cpow(z,b) == c)
                    {
                        z = limit;
                    }
                    else
                    {
                        z = (cpow(z,a) + c)/(cpow(z,b) - c);
                        n = n + 1;
                    }
                }
                pixels[y][x] = (int)(255 * (0.5 - 0.5 * cos(M_PI * n / 25.0)));
            }
        }
         
        // Write image to file
        sprintf(filename, "%03d.pgm", t);
        f = fopen(filename, "w");
        fprintf(stderr, "%s\n", filename);
        fprintf(f, "P5\n%d %d\n255\n", W, H);
        fwrite(pixels, W, H, f);
        fclose(f);
        
        // Convert PGM file to PNG format
        char command[255];
        sprintf(command, "mogrify -format png %s", filename);
        system(command);
        sprintf(command, "rm %s", filename);
        system(command);
    }
}

Posted in Uncategorized | Tagged , , , , , , , , | Leave a comment

Example code from presentation: Ways of Seeing Julia Sets

I gave a presentation at yesterday’s DIT School of Mathematical Sciences lunchtime research seminar titled Ways of Seeing Julia Sets: Visualizing the forces that shape fractal Julia sets. This was the abstract:

Early in the 20th century, work by mathematicians such as Pierre Fatou and Gaston Julia on complex dynamics led to the definition of so-called Julia sets. When a rational complex polynomial function is applied iteratively to a complex number, z, it produces a sequence of complex values called the orbit of z. Depending on the particular function and on the value of z, that orbit may or may not remain bounded. For a given function, the Julia set forms the boundary between those regions of the complex plane where the orbits remain bounded and those where they do not.Intriguingly, even for quite simple iterative complex functions, Julia sets often take on very striking fractal shapes. With the aid of a computer, it is easy to visualize the Julia set of a function but, for most people, understanding why it takes on a fractal shape is difficult. In my own ongoing struggle to gain a more intuitive understanding of fractals, I have written many short computer programs to visualize them in different ways. In this presentation, I will explore some ways of visualizing Julia sets which I found helpful in understanding why they are fractal.

filled_julia_set

My Prezi slides can be viewed here:
http://prezi.com/f5ymufi-bfqc/ways-of-seeing-julia-sets/

During my presentation, I demonstrated two short C programs that generated these fractal images.

This is the complete C code for the Mandelbrot example:

#include <stdio.h>
#include <complex.h>

void main()
{
    complex double z, c;
    int x, y, i;
    
    printf("P2\n4000 4000\n255\n");
    for (y=0 ; y<4000 ; ++y)
    {
        for (x=0 ; x<4000 ; ++x)
        {
            z = 0;
            c = (-2 + x*0.001) + (-2 + y*0.001)*I;
            i = 0;
            while (cabs(z) < 2 && ++i < 255) z = z*z + c;
            printf("%4d", i);
        }
        printf("\n");
    }
}

I compiled and ran the program as follows:

gcc mandelbrot.c -o mandelbrot -lm
./mandelbrot > mandelbrot.pgm

The C code for the Julia set example is almost identical except that it’s now z rather than c which changes from pixel to pixel (i.e. with changing x and y coordinates) and c is constant throughout the image (the value of c you choose determines which Julia set you end up with).

#include <stdio.h>
#include <complex.h>

void main()
{
    complex double z, c;
    int x, y, i;
    
    printf("P2\n4000 4000\n255\n");
    for (y=0 ; y<4000 ; ++y)
    {
        for (x=0 ; x<4000 ; ++x)
        {
            c = -0.625 - 0.4*I;
            z = (-2 + x*0.001) + (-2 + y*0.001)*I;
            i = 0;
            while (cabs(z) < 2 && ++i < 255) z = z*z + c;
            printf("%4d", i);
        }
        printf("\n");
    }
}

I compiled and ran the program as follows:

gcc julia.c -o julia -lm
./julia > julia.pgm
Posted in Uncategorized | Tagged , , , , , , , , | Leave a comment

C or Python? Comparison of execution time for Mandelbrot image generation

Because I’ve been generating a lot of fractal images recently, I’ve found myself spending quite a bit of time sitting around waiting for programs I’ve written to finish long repetitive calculations. In particular, I’ve generated several long sequences of several hundred images that I subsequently combined into animations.

One of the many great things about teaching electrical engineering is that I get to dip in and out of different programming languages all the time, but C remains my primary language and it’s still what I’m most likely to reach for when I’m writing something for my own curiosity. However, over the last couple of months I’ve been using Python for most of my fractal programming because I find its complex number syntax elegant. Furthermore, although complex numbers have been natively supported in C since the C99 standard, I’ve just never really embraced that language feature (until now).

Anyway, because I’ve been spending all this time waiting while numbers are getting crunched, I decided to do a quick comparison of execution time for the same task between Python and C. Naturally, since C is compiled and Python is interpreted, it comes as no surprise that C is significantly faster. Nevertheless, I was still surprised by how much faster it is.

The task I chose for this experiment was the generation of a Mandelbrot Set image of resolution 4000×4000 pixels, spanning a square region of the complex plane from -2 to +2 on the real axis and from -2j to +2j on the imaginary axis. This is the image produced:

mandelbrot

I implemented the same algorithm in C and Python and then performed a rough measurement of the execution time of each program by running the “date” command line tool immediately before and after it, as shown below…

Screenshot_2016-02-02_23-56-08

As you can see, the approximate execution times were as follows:

  • Python implementation: 105 seconds
  • C implementation: 4 seconds

In other words, the C implementation executes more than 26 times faster than the Python implementation! Compiling the program with gcc (also shown in the above screen shot) only takes a fraction of a second too. I’m sure a better programmer could make either of these implementations run a bit faster, but unless the Python version can run dramatically faster than this then the difference is more than enough to persuade me to embrace complex maths in C.

This is the full C code:

//
// mandelbrot.c - written by Ted Burke - 2-Feb-2016
//
// To compile:
//      gcc mandelbrot.c -o mandelbrot -lm
//

#include <complex.h>
#include <stdio.h>

// Create an array of bytes to store pixel values
unsigned char p[4000][4000];

void main()
{
    double complex c, z;
    int w=4000, h=4000, x, y;
    unsigned char px;
    
    // Open output file
    FILE* f = fopen("mandelbrot_c.pgm","w");
    fprintf(f, "P5\n%d %d\n255\n", w, h);
    
    // Calculate pixel values
    for (y=0 ; y<h ; ++y)
    {
        for (x=0 ; x<w ; ++x)
        {
            // Initialise complex values c and z for next
            // iterative pixel value calculation
            c = (-2.0 + x*4.0/w) + (2.0 - y*4.0/h)*I;
            z = 0;
            
            // Iterate complex function while darkening pixel
            px = 255;
            while(px >= 10 && cabs(z) < 4.0)
            {
                z = z*z + c;
                px -= 10;
            }
            
            // Store latest pixel value in byte array
            p[y][x] = px;
        }
    }
    
    // Write the byte array to the output file
    fwrite(p, 1, w*h, f);
    fclose(f);
}

This is the full Python code:

#
# mandelbrot.py - written by Ted Burke - 2-Feb-2016
#

import numpy

# Create an array of bytes to store pixel values
w,h = 4000,4000
p = numpy.zeros((h,w),dtype=numpy.uint8)

# Open output file
f = file('mandelbrot_python.pgm','w')
f.write('P5\n{:d} {:d}\n255\n'.format(w, h))

# Calculate pixel values
for y in range(h):
    for x in range(w):
        # Initialise complex values c and z for next
        # iterative pixel value calculation
        c = (-2.0 + x*4.0/w) + (2.0 - y*4.0/h)*(0+1j)
        z = 0 + 0j
        # Iterate complex function while darkening pixel
        px = 255
        while px >= 10 and abs(z) < 4.0:
            z = z*z + c
            px -= 10
        # Store latest pixel value in byte array
        p[y][x] = px

# Write the byte array to the output file
f.write(bytearray(p))
f.close()
Posted in Uncategorized | 6 Comments