#include <stdio.h>
#include <errno.h>
#include <unistd.h>
#include <string.h>
#include <stdlib.h>
#include <jack/jack.h>
#include <math.h>

static int mbuf[256][256];

static void draw_metaball(float x, float y, float radius)
{
    int cx,cy;
    float px,py;
    int *p = &mbuf[0][0];
    px = -x;
    py = -y;

    radius *= 400000.0f;

    for(cy=0; cy<256; cy++) {
        for(cx=0; cx<256;cx++) {
            float d = px*px+py*py;
            *p++ += radius/(d+1);
            px++;
        }
        py++;
        px = -x;
    }
}

void render_metaballs(uint8_t out[256][256], float time, float ltime)
{
    int obj = 0;
    int cx, cy;
    float dist1 = 128 + sinf(time * M_PI * 1.0 * 0.8 * 0.4 + 0) * 75;
    float dist2 = 128 + sinf(time * M_PI * 1.2 * 0.8 * 0.4 + 1) * 75;
    float dist3 = 130 + sinf(time * M_PI * 1.3 * 0.8 * 0.4 + 2) * 75;
    float dist4 = 100 + sinf(time * M_PI * 1.4 * 0.8 * 0.4 + 3) * 75;
    float dist5 = 128 + sinf(time * M_PI * 1.5 * 0.5 * 0.4 + 4) * 75;
    float dist6 = 128 + sinf(time * M_PI * 1.6 * 0.5 * 0.4 + 5) * 75;
    float dist7 = 130 + sinf(time * M_PI * 1.7 * 0.5 * 0.4 + 6) * 75;
    float dist8 = 100 + sinf(time * M_PI * 1.8 * 0.5 * 0.4 + 7) * 75;

    if (time < 1.0)
        dist1 -= 256-256*sinf((time) * M_PI / 2);
    if (time < 3.0)
        dist2 += 256-256*sinf((time-2.0) * M_PI / 2);
    if (time < 5.0) {
        dist7 -= 256-256*sinf((time-4.0) * M_PI / 2);
        dist3 += 256-256*sinf((time-4.0) * M_PI / 2);
    }
    if (time < 7.0) {
        dist8 -= 256-256*sinf((time-6.0) * M_PI / 2);
        dist4 -= 256-256*sinf((time-6.0) * M_PI / 2);
    }

    if (ltime < 7.0)
        dist1 -= 256-256*sinf((ltime-6.0) * M_PI / 2);
    if (ltime < 5.0)
        dist2 += 256-256*sinf((ltime-4.0) * M_PI / 2);
    if (ltime < 3.0) {
        dist7 -= 256-256*sinf((ltime-2.0) * M_PI / 2);
        dist3 += 256-256*sinf((ltime-2.0) * M_PI / 2);
    }
    if (ltime < 1.0) {
        dist8 -= 256-256*sinf((ltime) * M_PI / 2);
        dist4 -= 256-256*sinf((ltime) * M_PI / 2);
    }

    memset(mbuf, 0, sizeof mbuf);

    if (ltime > 6.0)
        draw_metaball(dist1, dist5, 45);
    if (time > 2.0 && ltime > 4.0)
        draw_metaball(dist2, dist6, 10);
    if (time > 4.0 && ltime > 2.0)
        draw_metaball(dist7, dist3, 30);
    if (time > 6.0)
        draw_metaball(dist8, dist4, 70);

    for(cy=0; cy<256; cy++) {
        for(cx=0; cx<256;cx++) {
            if(mbuf[cy][cx] > 25500)
                out[cy][cx] = 255;
            else
                out[cy][cx] = mbuf[cy][cx] / 100;
        }
    }
}

static inline int min(int a, int b)
{
    return a > b ? b : a;
}

static inline int max(int a, int b)
{
    return a < b ? b : a;
}

typedef struct {
    int x, y;
} Point;

static int orient2d(const Point a, const Point b, const Point c)
{
    return (b.x-a.x)*(c.y-a.y) - (b.y-a.y)*(c.x-a.x);
}

void render_triangle(float *buffer, int w, int h,
                     int x0, int y0, float z0,
                     int x1, int y1, float z1,
                     int x2, int y2, float z2)
{
    Point v0 = {x0, y0};
    Point v1 = {x1, y1};
    Point v2 = {x2, y2};

    int minX = max(0, min(x0, min(x1, x2)));
    int minY = max(7*h/32, min(y0, min(y1, y2)));
    int maxX = min(w - 1, max(x0, max(x1, x2)));
    int maxY = min(h - 1 - 7*h/32, max(y0, max(y1, y2)));

    int a = orient2d(v0, v1, v2);
    if (a == 0)
        return;
    Point p;
    for (p.y = minY; p.y <= maxY; p.y++) {
        for (p.x = minX; p.x <= maxX; p.x++) {
            // Determine barycentric coordinates
            int w0 = orient2d(v1, v2, p);
            int w1 = orient2d(v2, v0, p);
            int w2 = orient2d(v0, v1, p);

            if (w0 >= 0 && w1 >= 0 && w2 >= 0) {
                float v = z0 * w0 / (float)a +
                    z1 * w1 / (float)a +
                    z2 * w2 / (float)a;

                if (v > buffer[w*p.y + p.x])
                    buffer[w*p.y + p.x] = v;
            }
        }
    }    
}
