#include "config.h"
#include "common.h"

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

#include <pipi.h>

#define Y(x) (0.299*(x)[0] + 0.587*(x)[1] + 0.114*(x)[2])

#define STEPS 16
#define EPSILON (0.000001)

int main(int argc, char *argv[])
{
    double palette[][3] =
    {
        { 1.0, 0.0, 0.0 }, /* red */
        { 0.0, 1.0, 0.0 }, /* green */
        { 0.0, 0.0, 1.0 }, /* blue */
        { 1.0, 1.0, 1.0 }, /* white */
        { 1.0, 1.0, 0.0 }, /* yellow */
        { 1.0, 0.5, 0.0 }, /* orange */
    };
#define NCOLORS ((int)(sizeof(palette)/sizeof(*palette)))

    double tmp;
    int i, j;

    /*
     * 1. Find the darkest and lightest colours
     */
    double *dark = NULL, *light = NULL;
    double min = 1.0, max = 0.0;
    for(i = 0; i < NCOLORS; i++)
    {
        double p = Y(palette[i]);
        if(p < min)
        {
            dark = palette[i];
            min = p;
        }
        if(p > max)
        {
            light = palette[i];
            max = p;
        }
    }
    printf("lightest (%g,%g,%g)\n", light[0], light[1], light[2]);
    printf("darkest (%g,%g,%g)\n", dark[0], dark[1], dark[2]);

    /*
     * 2. Find two base vectors for the chrominance planes
     * FIXME: this doesn't work in all cases because u can be null
     */
    double y[3], u[3], v[3];
    double ylen;

    y[0] = light[0] - dark[0];
    y[1] = light[1] - dark[1];
    y[2] = light[2] - dark[2];
    ylen = sqrt(y[0] * y[0] + y[1] * y[1] + y[2] * y[2]);

    u[0] = y[1];
    u[1] = -y[0];
    u[2] = 0;
    tmp = sqrt(u[0] * u[0] + u[1] * u[1] + u[2] * u[2]);
    u[0] /= tmp; u[1] /= tmp; u[2] /= tmp;

    v[0] = y[1] * u[2] - y[2] * u[1];
    v[1] = y[2] * u[0] - y[0] * u[2];
    v[2] = y[0] * u[1] - y[1] * u[0];
    tmp = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
    v[0] /= tmp; v[1] /= tmp; v[2] /= tmp;

    printf("gray axis (%g,%g,%g) - length %g\n", y[0], y[1], y[2], ylen);
    printf("u (%g,%g,%g)\n", u[0], u[1], u[2]);
    printf("v (%g,%g,%g)\n", v[0], v[1], v[2]);

    /*
     * 3. Browse the grey axis and do stuff
     */
    double t;
    for(t = 0.; t <= 1.0; t += 1. / STEPS)
    {
        double pts[NCOLORS*NCOLORS/2][3];
        double p0[3];
        int npts = 0;

        p0[0] = dark[0] + t * y[0];
        p0[1] = dark[1] + t * y[1];
        p0[2] = dark[2] + t * y[2];
        printf("%g,%g,%g\n", p0[0], p0[1], p0[2]);

        /*
         * 3.1. Find all edges that intersect the t.y + (u,v) plane
         */
        for(i = 0; i < NCOLORS; i++)
        {
            double k1[3];
            k1[0] = palette[i][0] - dark[0];
            k1[1] = palette[i][1] - dark[1];
            k1[2] = palette[i][2] - dark[2];
            tmp = sqrt(k1[0] * k1[0] + k1[1] * k1[1] + k1[2] * k1[2]);

            /* If k1.y > t.y.y, we don't want this point */
            double yk1 = y[0] * k1[0] + y[1] * k1[1] + y[2] * k1[2];
            if(yk1 > t * ylen * ylen + EPSILON)
                continue;

            for(j = 0; j < NCOLORS; j++)
            {
                if(i == j)
                    continue;

                double k2[3];
                k2[0] = palette[j][0] - dark[0];
                k2[1] = palette[j][1] - dark[1];
                k2[2] = palette[j][2] - dark[2];
                tmp = sqrt(k2[0] * k2[0] + k2[1] * k2[1] + k2[2] * k2[2]);

                /* If k2.y < t.y.y, we don't want this point */
                double yk2 = y[0] * k2[0] + y[1] * k2[1] + y[2] * k2[2];
                if(yk2 < t * ylen * ylen - EPSILON)
                    continue;

                if(yk2 < yk1)
                    continue;

                double s = yk1 == yk2 ?
                           0.5 : (t * ylen * ylen - yk1) / (yk2 - yk1);

                pts[npts][0] = dark[0] + k1[0] + s * (k2[0] - k1[0]);
                pts[npts][1] = dark[1] + k1[1] + s * (k2[1] - k1[1]);
                pts[npts][2] = dark[2] + k1[2] + s * (k2[2] - k1[2]);
                printf("  %i,%i: s = (%g,%g,%g)\n", i, j, pts[npts][0], pts[npts][1], pts[npts][2]);
                npts++;
            }
        }

        /*
         * 3.2. Find the barycentre of these points' convex hull. (TODO)
         */
    }

    return 0;
}