Skip to content

Instantly share code, notes, and snippets.

@rmitton
Created August 30, 2020 02:27
Show Gist options
  • Save rmitton/114b732e48b9fec4024cb253561cb94d to your computer and use it in GitHub Desktop.
Save rmitton/114b732e48b9fec4024cb253561cb94d to your computer and use it in GitHub Desktop.
Mandelbrot set rendered from a first-person viewpoint
// I don't know what I'm looking at here.
#define _USE_MATH_DEFINES
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#define OUTPUT "hipparchus.pgm"
#define WIDTH 1024
#define HEIGHT 256
#define NPARTICLES 10000000
#define NITER 100
#define EXPOSURE 0.0001
uint64_t seed = 1;
double random( double a, double b ) // this is xorshift64, but any old random number source will work
{
seed ^= seed << 13;
seed ^= seed >> 7;
seed ^= seed << 17;
return a + seed * (b-a) / (double)(UINT64_MAX);
}
int main( int argc, char *argv[] )
{
int* histogram = calloc( WIDTH*HEIGHT, sizeof(int) );
for ( int i = 0; i<NPARTICLES; i++ )
{
again:;
// pick a random point on the 2D complex plane
double cx = random( -2, +2 );
double cy = random( -2, +2 );
// iterate it until we get bored
double zx = cx, zy = cy;
for ( int iter = 0; iter<NITER; iter++ )
{
// z' = z*z + c
double zx2 = zx*zx - zy*zy + cx;
double zy2 = zx*zy + zy*zx + cy;
zx = zx2;
zy = zy2;
// cartesian -> polar unwrap
double x = atan2( zy, zx );
double y = sqrt( zx*zx + zy*zy );
if ( y > 2 )
goto again; // reject non-stable orbits
// find pixel coordinates
int ix = (int)floor( x / ( 2.0*M_PI ) * WIDTH );
int iy = (int)floor( ( 1.0 - y / 2.0 ) * HEIGHT );
if ( ix < 0 )
ix += WIDTH;
// plot every point we pass through on our stable orbit
histogram[ iy*WIDTH+ix ]++;
}
}
// save image to PGM format
FILE* fp = fopen( OUTPUT, "wb" );
if ( !fp ) {
perror( OUTPUT );
return 1;
}
fprintf( fp, "P5 %i %i 255\n", WIDTH, HEIGHT );
for ( int i = 0; i<WIDTH*HEIGHT; i++ )
{
// apply simple camera exposure curve
int lum = (int)( 255 * ( 1 - exp( histogram[ i ] * -EXPOSURE ) ) );
fputc( lum, fp );
}
fclose( fp );
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment