_       _       
  (_) ___ | |_ ____
  | |/ _ \| __|_  /
  | | (_) | |_ / / 
 _/ |\___/ \__/___|
|__/               

Pencil sketch fractals

//
// barnsley.c - written by Ted Burke - last updated 12-July-2023
//
// To compile:
//   gcc barnsley.c -o barnsley -lm
//
// To run:
//   ./barnsley
//
// To convert output file to PNG:
//   convert image.pgm image.png
//

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

#define W 16000
#define H 16000

unsigned char p[H][W] = {0};      // oversized canvas for rendering
unsigned char pp[H/4][W/4][3] = {0}; // reduced canvas for final image

int main()
{
    // Generate image
    long i, j, n;
    double a, b, c, d, e, f;
    double x=-3, y=0, z;
    int xx, yy, q, r;
    int th = 16;

    // The following values will be used when apply a coordinate
    // transformation to each point prior to rendering to scale
    // and rotate the fern in the final image 
    double angle = -(35.0*M_PI/180.0);
    double sin_a, cos_a;
    sin_a = sin(angle);
    cos_a = cos(angle);

    // Iterate many times, plotting one point for each iteration
    for (i=0 ; i<400000000L ; ++i)
    {
        // Randomly select one of the four affine transforms
        n = rand()%100;

        // Leaf parameters
        if (n < 2)
            {th = 16; a=0.00; b=0.00; c=0.00; d=0.16; e=0.00; f=0.00;}
        else if (n < 86)
            {if(th)th--; a=0.85; b=0.02; c=-0.02; d=0.85; e=0.00; f=1.00;}
        else if (n < 93)
            {if(th)th--; a=0.4; b=-0.4; c=0.4; d=0.4; e=0.00; f=0.60;}
        else
            {if(th)th--; a=-0.4; b=0.4; c=0.4; d=0.4; e=0.00; f=0.70;}

        // Calculate next x,y point from the previous one
        z = a*x + b*y + e;
        y = c*x + d*y + f;
        x = z;

        // Map x,y point to pixel coordinates in image
        // (apply coordinate transformation to rotate and translate)
        xx = W/6 + W/10 + (W/7.0)*(cos_a*x - sin_a*y);
        yy = H/10 + H/10 + (H/7.0)*(sin_a*x + cos_a*y);
        yy = H - yy; // flip vertically

        // Mark that pixel (and neighbours if appropriate)
        for (q=yy-th/8 ; q<=yy+th/8 ; ++q)
            for (r=xx-th/8 ; r<=xx+th/8 ; ++r)
                if (0<=q && q<H && 0<=r && r<W) p[q][r] = 255;

        // Show progress updates while generating image
        if (i%10000000 == 0) fprintf(stderr, "%ld ", i/10000000);
    }

    fprintf(stderr, "\n");

    // Resample image to quarter the original resolution
    int px, ii, jj;
    for (ii=0 ; ii<H/4 ; ++ii) for (jj=0 ; jj<W/4 ; ++jj)
    {
        px = 0;
        for (i=0 ; i<4 ; ++i) for (j=0 ; j<4 ; ++j) px += p[4*ii+i][4*jj+j];
        px /= 16;        
        px = 255 - px;

        pp[ii][jj][0] = px;
        pp[ii][jj][1] = px;
        pp[ii][jj][2] = px;
    }

    // Write pixel to PNM file
    FILE *file = fopen("image.pgm", "w");
    fprintf(file, "P6\n%d %d\n255\n", W/4, H/4);
    fwrite(pp, 3*W/4, H/4, file);
    fclose(file);

    return 0;
}