Hypnagogic imagery

#include <cassert>
#include <cmath>
#include <cerrno>
#include <cstdio>
#include <cstring>
#include <random>
#include <stdexcept>
#include <tuple>
#include <vector>

extern "C" {
#include <png.h>
}

template<typename T>
T clamp(const T &x, const T &min, const T &max)
{
    if (x < min)
        return min;
    if (x > max)
        return max;

    return x;
}

struct vector {
    double x, y;

    explicit vector(double x, double y):
        x(x), y(y)
    {
    }

    double normsq() const
    {
        return x * x + y * y;
    }

    double norm() const
    {
        return sqrt(normsq());
    }

    double angle() const
    {
        return atan2(y, x);
    }

    void operator+=(const vector &other)
    {
        x += other.x;
        y += other.y;
    }
};

bool operator==(const vector &a, const vector &b)
{
    return a.x == b.x && a.y == b.y;
}

vector operator-(const vector &a, const vector &b)
{
    return vector(a.x - b.x, a.y - b.y);
}

static double bezier_interpolate(double p0, double p1, double t)
{
    return (1 - t) * p0 + t * p1;
}

static double bezier_interpolate(double p0, double p1, double p2, double t)
{
    double tp = (1 - t);
    return tp * tp * p0 + 2 * tp * t * p1 + t * t * p2;
}

static double bezier_interpolate(double p0, double p1, double p2, double p3, double t)
{
    double tp = 1 - t;
    return tp * tp * tp * p0
        + 3 * tp * tp * t * p1
        + 3 * tp * t * t * p2
        + t * t * t * p3;
}

static const unsigned int width = 320;
static const unsigned int height = 240;

static std::tuple<double, double, double> hsl_to_rgb(double H, double S, double L)
{
    double C = (1 - fabs(2 * L - 1)) * S;
    double X = C * (1 - fabs(fmod(H, 2) - 1));

    double r, g, b;
    if (H >= 0 && H < 1) {
        r = C; g = X; b = 0;
    } else if (H >= 1 && H < 2) {
        r = X; g = C; b = 0;
    } else if (H >= 2 && H < 3) {
        r = 0; g = C; b = X;
    } else if (H >= 3 && H < 4) {
        r = 0; g = X; b = C;
    } else if (H >= 4 && H < 5) {
        r = X; g = 0; b = C;
    } else if (H >= 5 && H <= 6) {
        r = C; g = 0; b = X;
    }

    double m = L - C / 2;

    return std::make_tuple(r + m, g + m, b + m);
}

template<unsigned int nr_corners>
struct blob {
    unsigned int x;
    unsigned int y;
    double size;

    double corners[nr_corners];
    double between[nr_corners];
};

template<typename Random, unsigned int nr_corners>
static void output_frame(Random &rand, const char *filename, const blob<nr_corners> &blue)
{
    png_byte data[width * height * 4];
    for (unsigned int y = 0; y < height; ++y) {
        for (unsigned int x = 0; x < width; ++x) {
            png_bytep rgb = &data[4 * (y * width + x)];
            rgb[0] = 0;
            rgb[1] = 0;
            rgb[2] = 0;
            rgb[3] = 255;
        }
    }

#define LESS

    for (unsigned int y = 0; y < height; ++y) {
        for (unsigned int x = 0; x < width; ++x) {
#ifdef LESS
            if (std::uniform_int_distribution<>(0, 100)(rand) < 50)
                continue;
#endif

            double H = 2 * std::uniform_int_distribution<>(0, 2)(rand);
            double S = std::uniform_real_distribution<>(0, 1)(rand);
#ifdef LESS
            double L = clamp(fabs(std::normal_distribution<>(0, 0.05)(rand)), 0., 1.);
#else
            double L = clamp(0.1 * std::exponential_distribution<>()(rand), 0., 1.);
#endif

            double r, g, b;
            std::tie(r, g, b) = hsl_to_rgb(H, S, L);

            png_bytep rgb = &data[4 * (y * width + x)];
            rgb[0] = 255 * r;
            rgb[1] = 255 * g;
            rgb[2] = 255 * b;
            rgb[3] = 255;
        }
    }

    unsigned int cx = blue.x;
    unsigned int cy = blue.y;
    unsigned int size = blue.size;

    for (unsigned int y = 0; y < 2 * size; ++y) {
        if (cy + y < size)
            continue;
        if (cy + y >= height + size)
            continue;

        for (unsigned int x = 0; x < 2 * size; ++x) {
            if (cx + x < size)
                continue;
            if (cx + x >= width + size)
                continue;

            vector p((double) y - size, (double) x - size);

            double corner = fmod(nr_corners * (M_PI + p.angle()) / 2 / M_PI, nr_corners);
            unsigned int corner_i = floor(corner);

            double t = corner - corner_i;
            assert(t >= 0 && t < 1);

            double norm = size * bezier_interpolate(
                blue.corners[corner_i],
                blue.corners[corner_i] + blue.between[corner_i],
                blue.corners[(corner_i + 1) % nr_corners] - blue.between[(corner_i + 1) % nr_corners],
                blue.corners[(corner_i + 1) % nr_corners],
                t);

            if (p.norm() > norm)
                continue;

            if (size > 100) {
                if (std::uniform_real_distribution<>(0, (p.norm() / norm))(rand) < (size - 100) / 500.)
                    continue;
            }

            double H = fmod(std::normal_distribution<>(4.4, .15)(rand), 6);
            double S = clamp(1., 0., 1.);
            double L;
            if (p.norm() > 0.98 * norm)
                L = clamp(std::normal_distribution<>(.5, .1)(rand), 0., 1.);
            else
                L = clamp(std::normal_distribution<>(.3, .1)(rand), 0., 1.);

            double r, g, b;
            std::tie(r, g, b) = hsl_to_rgb(H, S, L);

            png_bytep rgb = &data[4 * ((cy + y - size) * width + (cx + x - size))];
            rgb[0] = 255 * r;
            rgb[1] = 255 * g;
            rgb[2] = 255 * b;
            rgb[3] = 255;
        }
    }

    FILE *fp = fopen(filename, "wb");
    if (!fp)
        throw std::runtime_error(strerror(errno));

    png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, 0, 0, 0);
    if (!png_ptr)
        throw std::runtime_error("png_create_write_struct");

    png_infop info_ptr = png_create_info_struct(png_ptr);
    if (!info_ptr)
        throw std::runtime_error("png_create_info_struct");

    png_init_io(png_ptr, fp);
    png_set_filter(png_ptr, 0, PNG_FILTER_NONE);
    png_set_IHDR(png_ptr, info_ptr, width, height, 8, PNG_COLOR_TYPE_RGBA,
        PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT,
        PNG_FILTER_TYPE_DEFAULT);

    png_bytep rows[height];
    for (unsigned int i = 0; i < height; ++i)
        rows[i] = &data[i * width * 4];

    png_set_rows(png_ptr, info_ptr, rows);
    png_write_png(png_ptr, info_ptr, PNG_TRANSFORM_IDENTITY, NULL);
    fclose(fp);
}

int main(int argc, char *argv[])
{
    std::random_device rd;
    std::mt19937 rand(rd());

    unsigned int frame_i = 0;

    for (unsigned int i = 0; i < 3; ++i) {
        static const unsigned int nr_corners = 6;
        blob<nr_corners> blue;
        blue.x = std::normal_distribution<>(width / 2, width / 8)(rand);
        blue.y = std::normal_distribution<>(height / 2, height / 8)(rand);
        blue.size = 0;

        for (unsigned int i = 0; i < nr_corners; ++i)
            blue.corners[i] = std::uniform_real_distribution<>(.6, 1)(rand);
        for (unsigned int i = 0; i < nr_corners; ++i)
            blue.between[i] = std::uniform_real_distribution<>(.1, .3)(rand);

        double size_inc = 1;
        for (unsigned int j = 0; blue.size < width; ++j) {
            printf("Frame %u (size %.2f)\n", frame_i, blue.size);

            char png_filename[100];
            sprintf(png_filename, "output/frame-%03u.png", frame_i);

            char gif_filename[100];
            sprintf(gif_filename, "output/frame-%03u.gif", frame_i);

            output_frame(rand, png_filename, blue);
            ++frame_i;

            char command[100];
            sprintf(command, "convert %s %s && rm %s", png_filename, gif_filename, png_filename);
            system(command);

            size_inc *= 1.08;
            blue.size += size_inc;
        }
    }

    return 0;
}