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

/* NTSC artifact simulator, used in the postproduction of http://youtu.be/zQZLK8r1NFs
 * Copyright (C) 2012 Joel Yliluoma, http://iki.fi/bisqwit/
 */

typedef uint_least32_t u32;

static const unsigned Width  = 2048; //640;
static const unsigned Height = 360;
static const unsigned OutWidth = 1280;

static const unsigned Y_collect_width = 12;
static const unsigned U_collect_width = 31*2;
static const unsigned V_collect_width = 86*2;

int main()
{
    char Buf[1024];
    sprintf(Buf,
        "mencoder /mnt/gbatmp/binkvideo3.raw"
        " -demuxer rawvideo -rawvideo w=1920:h=1080:format=bgr24:fps=30"
      //  "mencoder tetre6_comb.avi -ss 20"
        " -sws 9 -aspect %u:%u -vf scale=%u:%u,format=rgb32 -nosound"
        " -ofps 30"
        " -really-quiet -ovc raw -of rawvideo -o -", Width,Height, Width,Height);
    FILE* fp = popen(Buf, "r");
    double colorburst = 0;

    // Create crude lowpass filters
    auto makeCoeff = [] (int size, double* buffer)
    {
        auto sinc = [] (double x) -> double
            { if(x == 0.) return 1.; x *= M_PI; return sin(x)/x; };
        double density = 0.;
        for(int q=-size/2, p=0; p<size; ++p,++q)
        {
            double x = q * 1. / size;
            buffer[p] = (size - 2*(q<0?-q:q)) + size*0.5 + 0.25*size*sinc(4*x);
            if(size == 12) buffer[p] = 1;
            density += buffer[p];
        }
        for(int p=0; p<size; ++p) buffer[p] /= density;
    };
    double coeffY[Y_collect_width]; makeCoeff(Y_collect_width, coeffY);
    double coeffU[U_collect_width]; makeCoeff(U_collect_width, coeffU);
    double coeffV[V_collect_width]; makeCoeff(V_collect_width, coeffV);
    auto lowpass = [] (double buffer[], int center, int width, const double coeff[]) -> double
    {
        int firsts = center-width/2;
        int firstp = 0;
        int limitp = width;
        if(firsts < 0)
            firstp = -firsts;
        if(firsts + limitp > (int)Width)
            limitp = Width - firsts;
        double result = 0.;
        for(int p = firstp; p < limitp; ++p)
            result += buffer[p+firsts] * coeff[p];
        return result;
    };

    for(;;)
    {
        static u32 Frame[Width*Height];
        unsigned bytes = sizeof(Frame), got=0;
        while(got < bytes)
        {
            int r = fread( ((char*)Frame) + got, 1, bytes-got, fp);
            if(r <= 0) goto done;
            got += r;
        }
        #define Gamma 1.0
        for(unsigned y=0; y<Height; ++y)
        {
            static double samples[Width];
            static double samplesU[Width];
            static double samplesV[Width];
            // Convert input scanline into NTSC signal.
            double SigFreq = (M_PI*2)/12.0;
            double offset = 2*M_PI / -6;
            for(unsigned x=0; x<Width; ++x)
            {
                u32 pix = Frame[x + y*Width];
                // Read RGB pixel.
                double r = pow(((pix >>  0) & 0xFF) / 255.0, 1/Gamma);
                double g = pow(((pix >>  8) & 0xFF) / 255.0, 1/Gamma);
                double b = pow(((pix >> 16) & 0xFF) / 255.0, 1/Gamma);

                // Convert to YIQ (the variables still call it YUV because I used YUV before).
                double Y = r * 0.299    + g * 0.587 + b * 0.114;
                double U = r * 0.596000 + g * -0.27400 + b * -0.3220;
                double V = r * 0.211000 + g * -0.52300 + b * 0.31200;

                //if(b > 0.5 && r < 0.5 && g < 0.5) { Y *= 2; }

                // Determine the UV angle.
                double angle = atan2(U,V), power = sqrt(U*U + V*V);
                // Phase-modulate the NTSC chroma signal.
                double phase = colorburst + angle + x*SigFreq;
                // Generate signal sample.
                // This is probably wrong, but I don't know how exactly it should be done.
                samples[x] = Y*((1-power) + power*cos(phase+offset));// < 0. ? -Y : Y;

                // Demodulate the chroma right away, so it can be used in lowpass().
                samplesU[x] = samples[x] * sin(colorburst + x*SigFreq);
                samplesV[x] = samples[x] * cos(colorburst + x*SigFreq);
            }

            // Convert NTSC signal into output scanline.
            static u32 out[OutWidth];
            for(unsigned x=0; x<OutWidth; ++x)
            {
                // Determine range for each component
                int center = x*Width/OutWidth;
            #if 1
                double Y = lowpass(samples, center, Y_collect_width, coeffY);
                double U = lowpass(samplesU,center, U_collect_width, coeffU);
                double V = lowpass(samplesV,center, V_collect_width, coeffV);
            #else
                int Ymin = center - Y_collect_width/2;
                int Ymax = center + (Y_collect_width-Y_collect_width/2);
                int Umin = center - U_collect_width/2;
                int Umax = center + (U_collect_width-U_collect_width/2);
                int Vmin = center - V_collect_width/2;
                int Vmax = center + (V_collect_width-V_collect_width/2);
                // Clamp ranges
                if(Ymin < 0) Ymin = 0;
                if(Umin < 0) Umin = 0;
                if(Vmin < 0) Vmin = 0;
                if(Ymax > int(Width)) Ymax = Width;
                if(Umax > int(Width)) Umax = Width;
                if(Vmax > int(Width)) Vmax = Width;
                // Collect samples
                double Y=0, U=0, V=0;
                double Ydiv = (Ymax-Ymin);
                double Udiv = (Umax-Umin);
                double Vdiv = (Vmax-Vmin);
                for(int p=Ymin; p<Ymax; ++p) Y += samples[p];
                for(int p=Umin; p<Umax; ++p) U += samples[p] * sin(colorburst + p*OutputFreq + offset);
                for(int p=Vmin; p<Vmax; ++p) V += samples[p] * cos(colorburst + p*OutputFreq + offset);
                Y /= Ydiv;
                U /= Udiv;
                V /= Vdiv;
            #endif
                U *= 3;
                V *= 3;
                // Now we have YIQ (even though the variables call it YUV). Convert to RGB.
                auto gammafix = [](double d) { return d<=0 ? 0 : pow(d, Gamma); };
                auto clamp = [](int v) { return v<0 ? 0 : v>255 ? 255 : v; };
                unsigned r = clamp(255.95 * gammafix(Y +  0.95620*U + -0.62140*V));
                unsigned g = clamp(255.95 * gammafix(Y + -0.27270*U + -0.64680*V));
                unsigned b = clamp(255.95 * gammafix(Y + -1.10370*U +  1.70060*V));
                out[x] = r*0x10000 + g*0x100 + b;
            }
            // Update scanline colorburst for signal.
            colorburst = fmod(colorburst + M_PI*2/3.0, M_PI*2);
            // Write scanline to output.
            bytes = sizeof(out); got = 0;
            while(got < bytes)
            {
                int r = fwrite( ((char*)out) + got, 1, bytes-got, stdout);
                if(r <= 0) goto done;
                got += r;
            }
        }
        colorburst += M_PI;
    }
done:;
}