/*****************************************************************************
                         YUVPAK  -  a program for
                         FRACTAL IMAGE COMPRESSION
                              Color Version
                          image.tga ==> image.ifs
                               version 2.0
*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <graphics.h>
#include <math.h>
#include <time.h>
#define number_flips 8
#define levels 2
#define max_patch 8
#define max_scale 1.2
long int usual (unsigned char Image[max_patch][max_patch], int size);
int mapping;
int main(int argc, char **argv)
{  /* begin main block */
    FILE *in, *out, *cf;
    char *inf, *outf, rmsstr[13];
    int ysize = 256,xsize = 256, min_patch = 4;
    unsigned char Image[ysize][xsize][3], blur[ysize-1][xsize-1], plotdx, plotdy;
    unsigned char Range[number_flips][max_patch][max_patch], Domain[max_patch][max_patch];
    unsigned char DX[xsize*ysize], DY[xsize*ysize], colormap[256][3], YUV, YYUV, UYUV, VYUV;
    int YUVTEMP,reverse;
    long int Domain_Class[levels][ysize - min_patch + 1][xsize - min_patch + 1][2][2];
    long int Range_Class[levels][64][64][2];
    int i, rx, ry, dx, dy, x, y, besti, bestdx, bestdy, patchsize, xx, yy;
    long int class1, class2, class3, class4, count;
    int inlevel;
    long int offset, best_offset;
    long int sumr, sumd, sumrd, sumr_sq, sumd_sq;
    long int Domain_sums[2][xsize - min_patch + 1][ysize - min_patch + 1];
    float fsumr, fsumd, fsumrd, fsumr_sq, fsumd_sq, fmagica;
    float best_scale, scale, root_mean_sq, best_root_mean_sq, min_variance;
    float mean_sq, root_mean_sq_tolerance, mean_sq_tolerance, best_mean_sq, fpatchsize_sq;
    float temp, variance, UF, VF, bf, gf, rf, uf, vf;
    time_t start_time, finish_time;
    long time_used;

    struct header_t {   /* "should" be a 12 byte header... we'll see */
      long time;        /* 4 bytes for compression time in seconds */
      short rms;        /* 2 bytes for 100.*rms value */
      short add1;       /* 2 bytes to be added later... room for growth */
      long add2;        /* 4 bytes to be added later... room for growth */
      } header[1];

    struct ifs_t {
      unsigned char dx;
      unsigned char dy;
      signed char scale : 7;
      short int offset : 7;
      unsigned short int flip : 3;
      unsigned short int size : 1;
      } ifs[1],
        ifs_table[levels][64][64];

    if (argc < 4)   {
        printf("usage: yuvpak rms infile.ext outfile.ext \n\n");
        printf("YUVPAK Version 2.0, Copyright (C) 1993, WD Young\n");
        printf("YUVPAK comes with ABSOLUTELY NO WARRANTY\n");
        printf("Please see files 'copying.wy' and 'copying' for details\n");
        printf("If these files are missing,\n");
        printf("write: WD Young, P.O. Box 632871, Nacogdoches TX 75963-2871\n");
        return 1;
        }


    root_mean_sq_tolerance = atof(argv[1]);
    header[0].rms = (short)(100.*root_mean_sq_tolerance);

    inf = argv[2]; outf = argv[3];
    if ((in = fopen(inf, "rb")) == NULL)    {
        fprintf(stderr, "Cannot open input file.\n");
        return 1;
        }
    if ((out = fopen(outf, "wb")) == NULL)  {
        fprintf(stderr, "Cannot open output file.\n");
        return 1;
        }
    fclose(out);

    start_time = time(start_time);
    min_variance = mean_sq_tolerance = root_mean_sq_tolerance*root_mean_sq_tolerance;

    GrSetMode(GR_default_graphics);
    
    for (y = 0; y < 64; y++)
        GrSetColor(y,4*y,4*y,4*y);
    for (y = 64; y < 256; y++)
        GrSetColor(y,y,0,0);

    for (y = 0; y < 18; y++) fgetc(in);

    for (y = 0; y < ysize; y++)
    for (x = 0; x < xsize; x++) {
        bf = fgetc(in);
        gf = fgetc(in);
        rf = fgetc(in);
        Image[y][x][0]  =  (unsigned char)(0.30*rf + 0.59*gf + 0.11*bf);
        uf =  0.62*rf - 0.52*gf - 0.10*bf;
        vf = -0.15*rf - 0.29*gf + 0.44*bf;
        Image[y][x][1] =  (unsigned char)(255.*(uf + 158.)/316.);
        Image[y][x][2] =  (unsigned char)(255.*(vf + 112.)/224.);
        }

 fclose(in);
    for (y = 0; y < ysize; y+=4)
    for (x = 0; x < xsize; x+=4) {
        UF = VF = 0;
        for (yy = 0; yy < 4; yy++)
        for (xx = 0; xx < 4; xx++) {
            UF += Image[y + yy][x + xx][1];
            VF += Image[y + yy][x + xx][2];
            }
        Image[y>>2][x>>2][1] = (unsigned char)(UF/16.);
        Image[y>>2][x>>2][2] = (unsigned char)(VF/16.);
        }

 for (YUV = 0; YUV < 3; YUV++) {
    min_patch = 4;
    if (YUV >= 1) {
        xsize = ysize = 64;
        mean_sq_tolerance = 0;
        min_patch = 4;
        }
    for (y = 0; y < ysize; y++)
    for (x = 0; x < xsize; x++) {
        GrPlot(x,y,Image[y][x][YUV]>>2);
        GrPlot(x+300,y,Image[y][x][YUV]>>2);
        }
    for (y = 0; y < ysize - 1; y++)
    for (x = 0; x < xsize - 1; x++)
        blur[y][x] =  (Image[y  ][x  ][YUV]
                     + Image[y  ][x+1][YUV]
                     + Image[y+1][x  ][YUV]
                     + Image[y+1][x+1][YUV])>>2;


    sprintf(rmsstr,"%4.1f",sqrt((double)mean_sq_tolerance));
    GrTextXY(20,290, "rms tolerance", 255, 0);
    GrTextXY(130,290,rmsstr,255,0);


/*************************************************************************/
/*                            Classify Range's                           */
/*************************************************************************/

    inlevel = 0;
    for (ry = 0; ry < ysize; ry+=8) {
        GrLine(0, ry, ysize - 1, ry, 384);
        for (rx = 0; rx < xsize; rx+=8)
            {
            for (y = 0; y < 8; y++)
            for (x = 0; x < 8; x++)
                Range[0][y][x] = Image [ ry  +y ] [ rx  +x ][YUV];
            class1 = usual(Range[0], 8);
            Range_Class[inlevel][ry>>3][rx>>3][0] = class1;
            }
        GrLine(0, ry, xsize - 1, ry, 384);
        }


/*
****************************************************************
*                                                              *
*  Compute Domain_sums array                                   *
*                                                              *
****************************************************************
*/
    inlevel = 1;
    patchsize = min_patch;
    for (dy = 0; dy < ysize - 2*patchsize+1; dy++) {
        GrLine(300, dy, 300 + xsize - 1, dy, 384);
        for (dx = 0; dx < xsize - 2*patchsize+1; dx++)  {
            Domain_sums[0][dy][dx] = Domain_sums[1][dy][dx] = 0;
            for (y = 0; y < 2*patchsize; y+=2)
            for (x = 0; x < 2*patchsize; x+=2)  {
                Domain[y>>1][x>>1] = blur[dy+y][dx+x];
                Domain_sums[0][dy][dx] += Domain[y>>1][x>>1];
                Domain_sums[1][dy][dx] += Domain[y>>1][x>>1]*Domain[y>>1][x>>1];
                }
            variance = (Domain_sums[1][dy][dx]- Domain_sums[0][dy][dx]*Domain_sums[0][dy][dx]/16.)/15.;
            if (variance > min_variance && YUV == 0) {
                class1 = usual(Domain, 4);
                Domain_Class[inlevel][dy][dx][0][0] = class1;
                Domain_Class[inlevel][dy][dx][0][1] = mapping;

                for (y = 0; y < 2*patchsize; y+=2)
                for (x = 0; x < 2*patchsize; x+=2)
                    Domain[y>>1][x>>1] = 255 - blur[dy+y][dx+x];
                class1 = usual(Domain, 4);
                Domain_Class[inlevel][dy][dx][1][0] = class1;
                Domain_Class[inlevel][dy][dx][1][1] = mapping;
                }
            else Domain_Class[inlevel][dy][dx][0][0] =
                 Domain_Class[inlevel][dy][dx][1][0] =-1;
            }
          GrLine(300, dy, 300 + xsize - 1, dy, 384);
        }
    
/*************************************************************************/
/*                            Classify Range's                           */
/*************************************************************************/
    for (ry = 0; ry < ysize; ry+=4) {
        GrLine(0, ry, xsize - 1, ry, 384);
        for (rx = 0; rx < xsize; rx+=4)
            {
            for (y = 0; y < 4; y++)
            for (x = 0; x < 4; x++)
                Range[0][y][x] = Image [ ry  +y ] [ rx  +x ][YUV];

            class1 = usual(Range[0], 4);
            Range_Class[inlevel][ry>>2][rx>>2][0] = class1;
            }
         GrLine(0, ry, xsize - 1, ry, 384);
        }
/*
****************************************************************
*                                                              *
*  Compute Domain_sums array                                   *
*                                                              *
****************************************************************
*/
    patchsize = max_patch;
    inlevel = 0;
    for (dy = 0; dy < ysize - 2*patchsize+1; dy++) {
        GrLine(300, dy, 300 + xsize - 1, dy, 384);
        for (dx = 0; dx < xsize - 2*patchsize+1; dx++)  {
            sumd = sumd_sq = 0;
            for (y = 0; y < 2*patchsize; y+=2*min_patch)
            for (x = 0; x < 2*patchsize; x+=2*min_patch) {
                sumd += Domain_sums[0][dy + y][dx + x];
                sumd_sq += Domain_sums[1][dy +  y][dx + x];
                }
            variance = (sumd_sq - sumd*sumd/64.)/63.;
            if (variance > min_variance && YUV == 0) {
                for (y = 0; y < 2*patchsize; y+=2)
                for (x = 0; x < 2*patchsize; x+=2)
                    Domain[y>>1][x>>1] = blur[dy+y][dx+x];

                class1 = usual(Domain, 8);
                Domain_Class[inlevel][dy][dx][0][0] = class1;
                Domain_Class[inlevel][dy][dx][0][1] = mapping;

                for (y = 0; y < 2*patchsize; y+=2)
                for (x = 0; x < 2*patchsize; x+=2)
                    Domain[y>>1][x>>1] = 255 - blur[dy+y][dx+x];

                class1 = usual(Domain, 8);
                Domain_Class[inlevel][dy][dx][1][0] = class1;
                Domain_Class[inlevel][dy][dx][1][1] = mapping;
                }
            else Domain_Class[inlevel][dy][dx][0][0] =
                 Domain_Class[inlevel][dy][dx][1][0] =-1;
        }
        GrLine(300, dy, 300 + xsize - 1, dy, 384);
        }
/*
****************************************************************
*                                                              *
*  Range Loop                                                  *
*                                                              *
****************************************************************
*/  
    inlevel = 0;
    for (patchsize = max_patch; patchsize >= min_patch; patchsize/=2, inlevel++)
    for (ry = 0; ry < ysize - patchsize + 1; ry+=patchsize)
    for (rx = 0; rx < xsize - patchsize + 1; rx+=patchsize)   {

        GrLine(rx, ry, rx, ry + patchsize, 384);
        GrLine(rx, ry ,rx + patchsize, ry, 384);
        GrLine(rx + patchsize, ry, rx + patchsize, ry + patchsize, 384);
        GrLine(rx, ry + patchsize, rx + patchsize, ry + patchsize, 384);

        sumr = sumr_sq = 0;
        for (y = 0; y < patchsize; y++)
        for (x = 0; x < patchsize; x++)
            {
            Range[0][y][x] =Image [ry              +y] [rx              +x][YUV];
            Range[1][y][x] =Image [ry+patchsize -1 -x] [rx              +y][YUV];
            Range[2][y][x] =Image [ry+patchsize -1 -y] [rx+patchsize -1 -x][YUV];
            Range[3][y][x] =Image [ry              +x] [rx+patchsize -1 -y][YUV];

            Range[4][y][x] =Image [ry+patchsize -1 -y] [rx              +x][YUV];
            Range[5][y][x] =Image [ry+patchsize -1 -x] [rx+patchsize -1 -y][YUV];
            Range[6][y][x] =Image [ry              +y] [rx+patchsize -1 -x][YUV];
            Range[7][y][x] =Image [ry              +x] [rx              +y][YUV];
            sumr+=Range[0][y][x]; sumr_sq+=Range[0][y][x]*Range[0][y][x];
            }
        fsumr=sumr; fsumr_sq=sumr_sq;

/*
****************************************************************
*                                                              *
*  Domain Loop                                                 *
*                                                              *
****************************************************************
*/
        if ((patchsize < max_patch)&&
            (ifs_table[inlevel-1][ry/(2*patchsize)][rx/(2*patchsize)].offset != (-500>>3))) goto cleanup;
        best_mean_sq = 10000000000.;
        count = 0;
        for (dy = 0; dy < ysize - 2*patchsize+1; dy++) {
        for (dx = 0; dx < xsize - 2*patchsize+1; dx++)
            if ((reverse = (Domain_Class[inlevel][dy][dx][0][0] == Range_Class[inlevel][ry/patchsize][rx/patchsize][0]))||
                (Domain_Class[inlevel][dy][dx][1][0] == Range_Class[inlevel][ry/patchsize][rx/patchsize][0])||
                (YUV >= 1))
            {
            GrPlot(dx+300+(patchsize>>1),dy+(patchsize>>1),384);
            DX[count] = dx; DY[count] = dy;
            count++;
            reverse = 1 - reverse;
            for (y = 0; y < (2*patchsize); y+=2)
            for (x = 0; x < (2*patchsize); x+=2)
                Domain[y>>1][x>>1] = blur[dy+y  ][dx+x  ];

            sumd = sumd_sq = 0;
            for (y = 0; y < 2*patchsize; y+=2*min_patch)
            for (x = 0; x < 2*patchsize; x+=2*min_patch) {
                sumd += Domain_sums[0][dy + y][dx + x];
                sumd_sq += Domain_sums[1][dy +  y][dx + x];
                }
            fsumd=sumd; fsumd_sq=sumd_sq;
            fpatchsize_sq = (float)(patchsize*patchsize);
            fmagica = (float)(sumd_sq - sumd*sumd/fpatchsize_sq);
            for (i = 0; i < number_flips; i++) {
                sumrd = 0;
                for (y = 0; y < patchsize; y++)
                for (x = 0; x < patchsize; x++)
                    sumrd += Domain[y][x]*Range[i][y][x];
                fsumrd = sumrd;
                if (fmagica != 0.)
                    scale = (fsumrd - fsumd*fsumr/fpatchsize_sq)/fmagica;
                else scale = 0;
                if (scale*scale < max_scale*max_scale) {
                    scale = (signed char) 63. * scale / max_scale;
                    scale = max_scale * scale / 63.;
                    offset = (long int)(fsumr - scale*fsumd)/fpatchsize_sq;
                    offset = (offset>>3)<<3;
                    mean_sq = (fsumr_sq + scale*(scale*fsumd_sq - 2*fsumrd + 2*offset*fsumd)
                        + offset*(offset*fpatchsize_sq - 2.*fsumr)) / fpatchsize_sq;

                    if (mean_sq < best_mean_sq) {
                        besti = i;  best_mean_sq = mean_sq; bestdx = dx; bestdy = dy;
                        best_scale = scale;  best_offset = offset;
                        }
                    if (mean_sq < mean_sq_tolerance) {
                        goto gotbest;
                        }
                    } /* end of conditional */
               }

            } /* end of Domain loop */
        }
          
gotbest:
        for (i = 0; i < count; i++)
                 GrPlot(DX[i]+300+(patchsize>>1),DY[i]+(patchsize>>1),384);

        sprintf(rmsstr,"%8.4f",sqrt((double)best_mean_sq));
        GrTextXY(560,20,rmsstr,255,0);
        ifs[0].dx = bestdx;
        ifs[0].dy = bestdy;
        ifs[0].flip = besti;
        ifs[0].scale = 63. * best_scale / max_scale;
        ifs[0].offset = (patchsize == min_patch||mean_sq < mean_sq_tolerance) ? best_offset>>3 : -500>>3;
        ifs_table[inlevel][ry/patchsize][rx/patchsize] = ifs[0];

cleanup:
        GrLine(rx, ry, rx, ry + patchsize, 384);
        GrLine(rx, ry ,rx + patchsize, ry, 384);
        GrLine(rx + patchsize, ry, rx + patchsize, ry + patchsize, 384);
        GrLine(rx, ry + patchsize, rx + patchsize, ry + patchsize, 384);

        }

        if ((out = fopen(outf, "ab")) == NULL)
            {
            fprintf(stderr, "Cannot open output file.\n");
            return 1;
            }

    if (YUV == 2) {
        finish_time = time(finish_time);
        time_used = (long)difftime(finish_time, start_time);
        header[0].time = time_used;
        fwrite(header, sizeof(struct header_t), 1, out);
        }

    for (ry = 0; ry < ysize>>3; ry++)
    for (rx = 0; rx < xsize>>3; rx++)  {
        if (ifs_table[0][ry][rx].offset == (-500>>3))
            for (y = 2*ry; y < 2*ry + 2; y++)
            for (x = 2*rx; x < 2*rx + 2; x++) {
                ifs[0] = ifs_table[1][y][x];
                ifs[0].size = 1;
                fwrite(ifs, sizeof(struct ifs_t), 1, out);
            }
        else {
            ifs[0] = ifs_table[0][ry][rx];
            ifs[0].size = 0;
            fwrite(ifs, sizeof(struct ifs_t), 1, out);
            }
        }
            
   fclose(out);
 }  
   GrSetMode(GR_default_text);
/* All done. Whew... */
   return 0;

}  /* end main block */
#include "wdyusual.c"