#include <stdio.h>
#include <stdlib.h>
#include "pnginput.h"
#include "pngoutput.h"
#include "colourspace.h"
#include "triangulation.h"
#include "mcrandom.h"

#define USE_LUV
// #define SUPERSAMPLE
// #define SUPERDEPTH 4

void read_image( double *data, double *pixeldata,
                 pnginput_t *pnginput, int width, int height )
{
    int i, j;

    for( i = 0; i < height; i++ ) {
        unsigned char *curscanline = pnginput_get_scanline( pnginput, i );

        for( j = 0; j < width; j++ ) {
            double bp = ( (double) curscanline[ 0 ] ) / 255.0;
            double gp = ( (double) curscanline[ 1 ] ) / 255.0;
            double rp = ( (double) curscanline[ 2 ] ) / 255.0;
            double r, g, b, x, y, z, lstar, ustar, vstar;

            nonlinear_to_linear_srgb( rp, gp, bp, &r, &g, &b );
            srgb_to_xyz( r, g, b, &x, &y, &z );
            xyz_to_lupvp( x, y, z, &lstar, &ustar, &vstar );

            *data++ = lstar;
#ifdef USE_LUV
            *pixeldata++ = lstar;
            *pixeldata++ = ustar;
            *pixeldata++ = vstar;
#else
            *pixeldata++ = x;
            *pixeldata++ = y;
            *pixeldata++ = z;
#endif
            curscanline += 3;
        }
    }
}

void supersample_point( triangulation_t *image, int depth,
                        double curx, double cury, double w, double h,
                        double owidth, double oheight,
                        double *a, double *b, double *c, int maxdepth )
{
    double cura, curb, curc;
    double left, right, top;

    if( depth > maxdepth ) {
        double xmod = ( ( random01() * 2.0 ) - 1.0 ) / ( 2.0 * ((double) depth) );
        double ymod = ( ( random01() * 2.0 ) - 1.0 ) / ( 2.0 * ((double) depth) );
        interpolate_value( image, ((curx + xmod) / w) * owidth,
                                  ((cury + ymod) / h) * oheight,
                                  a, b, c );
        return;
    }

    *a = *b = *c = 0.0;

    left  = curx - ( 1.0 / ( 2.0 * ((double) depth) ) );
    right = curx + ( 1.0 / ( 2.0 * ((double) depth) ) );
    top   = cury + ( 1.0 / ( 2.0 * ((double) depth) ) );
    
    supersample_point( image, depth + 1,
                       left + ( ( curx - left ) / 2.0 ),
                       cury + ( ( top - cury ) / 2.0 ),
                       w, h, owidth, oheight, &cura, &curb, &curc, maxdepth );
    *a += cura; *b += curb; *c += curc;

    supersample_point( image, depth + 1,
                       curx + ( ( right - curx ) / 2.0 ),
                       cury + ( ( top - cury ) / 2.0 ),
                       w, h, owidth, oheight, &cura, &curb, &curc, maxdepth );
    *a += cura; *b += curb; *c += curc;

    supersample_point( image, depth + 1,
                       left + ( ( curx - left ) / 2.0 ),
                       cury - ( ( top - cury ) / 2.0 ),
                       w, h, owidth, oheight, &cura, &curb, &curc, maxdepth );
    *a += cura; *b += curb; *c += curc;

    supersample_point( image, depth + 1,
                       curx + ( ( right - curx ) / 2.0 ),
                       cury - ( ( top - cury ) / 2.0 ),
                       w, h, owidth, oheight, &cura, &curb, &curc, maxdepth );
    *a += cura; *b += curb; *c += curc;

    *a = *a * 0.25; *b = *b * 0.25; *c = *c * 0.25;
}

void output_image_rgb( pngoutput_t *pngoutput, triangulation_t *image,
                       int width, int height )
{
    double owidth = (double) (get_triangulation_width( image ) - 1);
    double oheight = (double) (get_triangulation_height( image ) - 1);
    unsigned char *scanline = (unsigned char *) malloc( width * 3 );
    int i, j;

    for( i = 0; i < height; i++ ) {
        unsigned char *pixeldata = scanline;
        fprintf( stderr, "\rscale: Output scanline: %4d", i );
        for( j = 0; j < width; j++ ) {
            double curx = (((double) j) / ((double) (width-1))) * owidth;
            double cury = (((double) i) / ((double) (height-1))) * oheight;
            double lstar, ustar, vstar, x, y, z, r, g, b, rp, gp, bp;

#ifdef USE_LUV
#ifdef SUPERSAMPLE
            supersample_point( image, 1, (double) j, (double) i,
                               (double) (width-1), (double) (height-1),
                               owidth, oheight, &lstar, &ustar, &vstar,
                               SUPERDEPTH );
#else
            interpolate_value( image, curx, cury, &lstar, &ustar, &vstar );
#endif /* SUPERSAMPLE */
            lupvp_to_xyz( lstar, ustar, vstar, &x, &y, &z );
#else  /* !USE_LUV */
            interpolate_value( image, curx, cury, &x, &y, &z );
#endif /* USE_LUV */
            xyz_to_srgb( x, y, z, &r, &g, &b );
            linear_to_nonlinear_srgb( r, g, b, &rp, &gp, &bp );

            *pixeldata++ = (unsigned char) ( ( bp * 255.0 ) + 0.5 );
            *pixeldata++ = (unsigned char) ( ( gp * 255.0 ) + 0.5 );
            *pixeldata++ = (unsigned char) ( ( rp * 255.0 ) + 0.5 );
        }
        pngoutput_scanline( pngoutput, scanline );
    }
    fprintf( stderr, "\n" );
    free( scanline );
}

int main( int argc, char **argv )
{
    pngoutput_t *pngoutput;
    pnginput_t *pnginput;
    int width, height;
    double *data;
    double *pixeldata;
    triangulation_t *image;
    int factor;

    pnginput = pnginput_new( argv[ 1 ] );
    if( !pnginput ) {
        fprintf( stderr, "scale: can't open input %s\n", argv[ 1 ] );
        return 1;
    }
    width = pnginput_get_width( pnginput );
    height = pnginput_get_height( pnginput );
    data = (double *) malloc( width * height * sizeof( double ) );
    if( !data ) {
        fprintf( stderr, "scale: can't allocate input data\n" );
        return 1;
    }
    pixeldata = (double *) malloc( 3 * width * height * sizeof( double ) );
    if( !pixeldata ) {
        fprintf( stderr, "scale: can't allocate pixel data\n" );
        return 1;
    }

    read_image( data, pixeldata, pnginput, width, height );
    pnginput_delete( pnginput );

    factor = atoi( argv[ 4 ] );
    fprintf( stderr, "scale: input is %dx%d\n", width, height );
    fprintf( stderr, "scale: scaling by a factor of %d: output %dx%d\n",
             factor, width * factor, height * factor );

    pngoutput = pngoutput_new( argv[ 2 ], width * factor, height * factor, 0.45 );
    if( !pngoutput ) {
        fprintf( stderr, "scale: can't open output %s\n", argv[ 2 ] );
        return 1;
    }

    image = initialize_image( data, pixeldata, width, height );
    load_triangulation( image, argv[ 3 ] );
    output_image_rgb( pngoutput, image, width * factor, height * factor );
    pngoutput_delete( pngoutput );

    return 0;
}