home *** CD-ROM | disk | FTP | other *** search
/ DP Tool Club 9 / CD_ASCQ_09_1193.iso / news / 4441 / yuvpak / yuvpak.c < prev    next >
C/C++ Source or Header  |  1993-01-30  |  17KB  |  417 lines
  1. /*****************************************************************************
  2.                          YUVPAK  -  a program for
  3.                          FRACTAL IMAGE COMPRESSION
  4.                               Color Version
  5.                           image.tga ==> image.ifs
  6.                                version 2.0
  7. *****************************************************************************/
  8. #include <stdio.h>
  9. #include <stdlib.h>
  10. #include <graphics.h>
  11. #include <math.h>
  12. #include <time.h>
  13. #define number_flips 8
  14. #define levels 2
  15. #define max_patch 8
  16. #define max_scale 1.2
  17. long int usual (unsigned char Image[max_patch][max_patch], int size);
  18. int mapping;
  19. int main(int argc, char **argv)
  20. {  /* begin main block */
  21.     FILE *in, *out, *cf;
  22.     char *inf, *outf, rmsstr[13];
  23.     int ysize = 256,xsize = 256, min_patch = 4;
  24.     unsigned char Image[ysize][xsize][3], blur[ysize-1][xsize-1], plotdx, plotdy;
  25.     unsigned char Range[number_flips][max_patch][max_patch], Domain[max_patch][max_patch];
  26.     unsigned char DX[xsize*ysize], DY[xsize*ysize], colormap[256][3], YUV, YYUV, UYUV, VYUV;
  27.     int YUVTEMP,reverse;
  28.     long int Domain_Class[levels][ysize - min_patch + 1][xsize - min_patch + 1][2][2];
  29.     long int Range_Class[levels][64][64][2];
  30.     int i, rx, ry, dx, dy, x, y, besti, bestdx, bestdy, patchsize, xx, yy;
  31.     long int class1, class2, class3, class4, count;
  32.     int inlevel;
  33.     long int offset, best_offset;
  34.     long int sumr, sumd, sumrd, sumr_sq, sumd_sq;
  35.     long int Domain_sums[2][xsize - min_patch + 1][ysize - min_patch + 1];
  36.     float fsumr, fsumd, fsumrd, fsumr_sq, fsumd_sq, fmagica;
  37.     float best_scale, scale, root_mean_sq, best_root_mean_sq, min_variance;
  38.     float mean_sq, root_mean_sq_tolerance, mean_sq_tolerance, best_mean_sq, fpatchsize_sq;
  39.     float temp, variance, UF, VF, bf, gf, rf, uf, vf;
  40.     time_t start_time, finish_time;
  41.     long time_used;
  42.  
  43.     struct header_t {   /* "should" be a 12 byte header... we'll see */
  44.       long time;        /* 4 bytes for compression time in seconds */
  45.       short rms;        /* 2 bytes for 100.*rms value */
  46.       short add1;       /* 2 bytes to be added later... room for growth */
  47.       long add2;        /* 4 bytes to be added later... room for growth */
  48.       } header[1];
  49.  
  50.     struct ifs_t {
  51.       unsigned char dx;
  52.       unsigned char dy;
  53.       signed char scale : 7;
  54.       short int offset : 7;
  55.       unsigned short int flip : 3;
  56.       unsigned short int size : 1;
  57.       } ifs[1],
  58.         ifs_table[levels][64][64];
  59.  
  60.     if (argc < 4)   {
  61.         printf("usage: yuvpak rms infile.ext outfile.ext \n\n");
  62.         printf("YUVPAK Version 2.0, Copyright (C) 1993, WD Young\n");
  63.         printf("YUVPAK comes with ABSOLUTELY NO WARRANTY\n");
  64.         printf("Please see files 'copying.wy' and 'copying' for details\n");
  65.         printf("If these files are missing,\n");
  66.         printf("write: WD Young, P.O. Box 632871, Nacogdoches TX 75963-2871\n");
  67.         return 1;
  68.         }
  69.  
  70.  
  71.     root_mean_sq_tolerance = atof(argv[1]);
  72.     header[0].rms = (short)(100.*root_mean_sq_tolerance);
  73.  
  74.     inf = argv[2]; outf = argv[3];
  75.     if ((in = fopen(inf, "rb")) == NULL)    {
  76.         fprintf(stderr, "Cannot open input file.\n");
  77.         return 1;
  78.         }
  79.     if ((out = fopen(outf, "wb")) == NULL)  {
  80.         fprintf(stderr, "Cannot open output file.\n");
  81.         return 1;
  82.         }
  83.     fclose(out);
  84.  
  85.     start_time = time(start_time);
  86.     min_variance = mean_sq_tolerance = root_mean_sq_tolerance*root_mean_sq_tolerance;
  87.  
  88.     GrSetMode(GR_default_graphics);
  89.     
  90.     for (y = 0; y < 64; y++)
  91.         GrSetColor(y,4*y,4*y,4*y);
  92.     for (y = 64; y < 256; y++)
  93.         GrSetColor(y,y,0,0);
  94.  
  95.     for (y = 0; y < 18; y++) fgetc(in);
  96.  
  97.     for (y = 0; y < ysize; y++)
  98.     for (x = 0; x < xsize; x++) {
  99.         bf = fgetc(in);
  100.         gf = fgetc(in);
  101.         rf = fgetc(in);
  102.         Image[y][x][0]  =  (unsigned char)(0.30*rf + 0.59*gf + 0.11*bf);
  103.         uf =  0.62*rf - 0.52*gf - 0.10*bf;
  104.         vf = -0.15*rf - 0.29*gf + 0.44*bf;
  105.         Image[y][x][1] =  (unsigned char)(255.*(uf + 158.)/316.);
  106.         Image[y][x][2] =  (unsigned char)(255.*(vf + 112.)/224.);
  107.         }
  108.  
  109.  fclose(in);
  110.     for (y = 0; y < ysize; y+=4)
  111.     for (x = 0; x < xsize; x+=4) {
  112.         UF = VF = 0;
  113.         for (yy = 0; yy < 4; yy++)
  114.         for (xx = 0; xx < 4; xx++) {
  115.             UF += Image[y + yy][x + xx][1];
  116.             VF += Image[y + yy][x + xx][2];
  117.             }
  118.         Image[y>>2][x>>2][1] = (unsigned char)(UF/16.);
  119.         Image[y>>2][x>>2][2] = (unsigned char)(VF/16.);
  120.         }
  121.  
  122.  for (YUV = 0; YUV < 3; YUV++) {
  123.     min_patch = 4;
  124.     if (YUV >= 1) {
  125.         xsize = ysize = 64;
  126.         mean_sq_tolerance = 0;
  127.         min_patch = 4;
  128.         }
  129.     for (y = 0; y < ysize; y++)
  130.     for (x = 0; x < xsize; x++) {
  131.         GrPlot(x,y,Image[y][x][YUV]>>2);
  132.         GrPlot(x+300,y,Image[y][x][YUV]>>2);
  133.         }
  134.     for (y = 0; y < ysize - 1; y++)
  135.     for (x = 0; x < xsize - 1; x++)
  136.         blur[y][x] =  (Image[y  ][x  ][YUV]
  137.                      + Image[y  ][x+1][YUV]
  138.                      + Image[y+1][x  ][YUV]
  139.                      + Image[y+1][x+1][YUV])>>2;
  140.  
  141.  
  142.     sprintf(rmsstr,"%4.1f",sqrt((double)mean_sq_tolerance));
  143.     GrTextXY(20,290, "rms tolerance", 255, 0);
  144.     GrTextXY(130,290,rmsstr,255,0);
  145.  
  146.  
  147. /*************************************************************************/
  148. /*                            Classify Range's                           */
  149. /*************************************************************************/
  150.  
  151.     inlevel = 0;
  152.     for (ry = 0; ry < ysize; ry+=8) {
  153.         GrLine(0, ry, ysize - 1, ry, 384);
  154.         for (rx = 0; rx < xsize; rx+=8)
  155.             {
  156.             for (y = 0; y < 8; y++)
  157.             for (x = 0; x < 8; x++)
  158.                 Range[0][y][x] = Image [ ry  +y ] [ rx  +x ][YUV];
  159.             class1 = usual(Range[0], 8);
  160.             Range_Class[inlevel][ry>>3][rx>>3][0] = class1;
  161.             }
  162.         GrLine(0, ry, xsize - 1, ry, 384);
  163.         }
  164.  
  165.  
  166. /*
  167. ****************************************************************
  168. *                                                              *
  169. *  Compute Domain_sums array                                   *
  170. *                                                              *
  171. ****************************************************************
  172. */
  173.     inlevel = 1;
  174.     patchsize = min_patch;
  175.     for (dy = 0; dy < ysize - 2*patchsize+1; dy++) {
  176.         GrLine(300, dy, 300 + xsize - 1, dy, 384);
  177.         for (dx = 0; dx < xsize - 2*patchsize+1; dx++)  {
  178.             Domain_sums[0][dy][dx] = Domain_sums[1][dy][dx] = 0;
  179.             for (y = 0; y < 2*patchsize; y+=2)
  180.             for (x = 0; x < 2*patchsize; x+=2)  {
  181.                 Domain[y>>1][x>>1] = blur[dy+y][dx+x];
  182.                 Domain_sums[0][dy][dx] += Domain[y>>1][x>>1];
  183.                 Domain_sums[1][dy][dx] += Domain[y>>1][x>>1]*Domain[y>>1][x>>1];
  184.                 }
  185.             variance = (Domain_sums[1][dy][dx]- Domain_sums[0][dy][dx]*Domain_sums[0][dy][dx]/16.)/15.;
  186.             if (variance > min_variance && YUV == 0) {
  187.                 class1 = usual(Domain, 4);
  188.                 Domain_Class[inlevel][dy][dx][0][0] = class1;
  189.                 Domain_Class[inlevel][dy][dx][0][1] = mapping;
  190.  
  191.                 for (y = 0; y < 2*patchsize; y+=2)
  192.                 for (x = 0; x < 2*patchsize; x+=2)
  193.                     Domain[y>>1][x>>1] = 255 - blur[dy+y][dx+x];
  194.                 class1 = usual(Domain, 4);
  195.                 Domain_Class[inlevel][dy][dx][1][0] = class1;
  196.                 Domain_Class[inlevel][dy][dx][1][1] = mapping;
  197.                 }
  198.             else Domain_Class[inlevel][dy][dx][0][0] =
  199.                  Domain_Class[inlevel][dy][dx][1][0] =-1;
  200.             }
  201.           GrLine(300, dy, 300 + xsize - 1, dy, 384);
  202.         }
  203.     
  204. /*************************************************************************/
  205. /*                            Classify Range's                           */
  206. /*************************************************************************/
  207.     for (ry = 0; ry < ysize; ry+=4) {
  208.         GrLine(0, ry, xsize - 1, ry, 384);
  209.         for (rx = 0; rx < xsize; rx+=4)
  210.             {
  211.             for (y = 0; y < 4; y++)
  212.             for (x = 0; x < 4; x++)
  213.                 Range[0][y][x] = Image [ ry  +y ] [ rx  +x ][YUV];
  214.  
  215.             class1 = usual(Range[0], 4);
  216.             Range_Class[inlevel][ry>>2][rx>>2][0] = class1;
  217.             }
  218.          GrLine(0, ry, xsize - 1, ry, 384);
  219.         }
  220. /*
  221. ****************************************************************
  222. *                                                              *
  223. *  Compute Domain_sums array                                   *
  224. *                                                              *
  225. ****************************************************************
  226. */
  227.     patchsize = max_patch;
  228.     inlevel = 0;
  229.     for (dy = 0; dy < ysize - 2*patchsize+1; dy++) {
  230.         GrLine(300, dy, 300 + xsize - 1, dy, 384);
  231.         for (dx = 0; dx < xsize - 2*patchsize+1; dx++)  {
  232.             sumd = sumd_sq = 0;
  233.             for (y = 0; y < 2*patchsize; y+=2*min_patch)
  234.             for (x = 0; x < 2*patchsize; x+=2*min_patch) {
  235.                 sumd += Domain_sums[0][dy + y][dx + x];
  236.                 sumd_sq += Domain_sums[1][dy +  y][dx + x];
  237.                 }
  238.             variance = (sumd_sq - sumd*sumd/64.)/63.;
  239.             if (variance > min_variance && YUV == 0) {
  240.                 for (y = 0; y < 2*patchsize; y+=2)
  241.                 for (x = 0; x < 2*patchsize; x+=2)
  242.                     Domain[y>>1][x>>1] = blur[dy+y][dx+x];
  243.  
  244.                 class1 = usual(Domain, 8);
  245.                 Domain_Class[inlevel][dy][dx][0][0] = class1;
  246.                 Domain_Class[inlevel][dy][dx][0][1] = mapping;
  247.  
  248.                 for (y = 0; y < 2*patchsize; y+=2)
  249.                 for (x = 0; x < 2*patchsize; x+=2)
  250.                     Domain[y>>1][x>>1] = 255 - blur[dy+y][dx+x];
  251.  
  252.                 class1 = usual(Domain, 8);
  253.                 Domain_Class[inlevel][dy][dx][1][0] = class1;
  254.                 Domain_Class[inlevel][dy][dx][1][1] = mapping;
  255.                 }
  256.             else Domain_Class[inlevel][dy][dx][0][0] =
  257.                  Domain_Class[inlevel][dy][dx][1][0] =-1;
  258.         }
  259.         GrLine(300, dy, 300 + xsize - 1, dy, 384);
  260.         }
  261. /*
  262. ****************************************************************
  263. *                                                              *
  264. *  Range Loop                                                  *
  265. *                                                              *
  266. ****************************************************************
  267. */  
  268.     inlevel = 0;
  269.     for (patchsize = max_patch; patchsize >= min_patch; patchsize/=2, inlevel++)
  270.     for (ry = 0; ry < ysize - patchsize + 1; ry+=patchsize)
  271.     for (rx = 0; rx < xsize - patchsize + 1; rx+=patchsize)   {
  272.  
  273.         GrLine(rx, ry, rx, ry + patchsize, 384);
  274.         GrLine(rx, ry ,rx + patchsize, ry, 384);
  275.         GrLine(rx + patchsize, ry, rx + patchsize, ry + patchsize, 384);
  276.         GrLine(rx, ry + patchsize, rx + patchsize, ry + patchsize, 384);
  277.  
  278.         sumr = sumr_sq = 0;
  279.         for (y = 0; y < patchsize; y++)
  280.         for (x = 0; x < patchsize; x++)
  281.             {
  282.             Range[0][y][x] =Image [ry              +y] [rx              +x][YUV];
  283.             Range[1][y][x] =Image [ry+patchsize -1 -x] [rx              +y][YUV];
  284.             Range[2][y][x] =Image [ry+patchsize -1 -y] [rx+patchsize -1 -x][YUV];
  285.             Range[3][y][x] =Image [ry              +x] [rx+patchsize -1 -y][YUV];
  286.  
  287.             Range[4][y][x] =Image [ry+patchsize -1 -y] [rx              +x][YUV];
  288.             Range[5][y][x] =Image [ry+patchsize -1 -x] [rx+patchsize -1 -y][YUV];
  289.             Range[6][y][x] =Image [ry              +y] [rx+patchsize -1 -x][YUV];
  290.             Range[7][y][x] =Image [ry              +x] [rx              +y][YUV];
  291.             sumr+=Range[0][y][x]; sumr_sq+=Range[0][y][x]*Range[0][y][x];
  292.             }
  293.         fsumr=sumr; fsumr_sq=sumr_sq;
  294.  
  295. /*
  296. ****************************************************************
  297. *                                                              *
  298. *  Domain Loop                                                 *
  299. *                                                              *
  300. ****************************************************************
  301. */
  302.         if ((patchsize < max_patch)&&
  303.             (ifs_table[inlevel-1][ry/(2*patchsize)][rx/(2*patchsize)].offset != (-500>>3))) goto cleanup;
  304.         best_mean_sq = 10000000000.;
  305.         count = 0;
  306.         for (dy = 0; dy < ysize - 2*patchsize+1; dy++) {
  307.         for (dx = 0; dx < xsize - 2*patchsize+1; dx++)
  308.             if ((reverse = (Domain_Class[inlevel][dy][dx][0][0] == Range_Class[inlevel][ry/patchsize][rx/patchsize][0]))||
  309.                 (Domain_Class[inlevel][dy][dx][1][0] == Range_Class[inlevel][ry/patchsize][rx/patchsize][0])||
  310.                 (YUV >= 1))
  311.             {
  312.             GrPlot(dx+300+(patchsize>>1),dy+(patchsize>>1),384);
  313.             DX[count] = dx; DY[count] = dy;
  314.             count++;
  315.             reverse = 1 - reverse;
  316.             for (y = 0; y < (2*patchsize); y+=2)
  317.             for (x = 0; x < (2*patchsize); x+=2)
  318.                 Domain[y>>1][x>>1] = blur[dy+y  ][dx+x  ];
  319.  
  320.             sumd = sumd_sq = 0;
  321.             for (y = 0; y < 2*patchsize; y+=2*min_patch)
  322.             for (x = 0; x < 2*patchsize; x+=2*min_patch) {
  323.                 sumd += Domain_sums[0][dy + y][dx + x];
  324.                 sumd_sq += Domain_sums[1][dy +  y][dx + x];
  325.                 }
  326.             fsumd=sumd; fsumd_sq=sumd_sq;
  327.             fpatchsize_sq = (float)(patchsize*patchsize);
  328.             fmagica = (float)(sumd_sq - sumd*sumd/fpatchsize_sq);
  329.             for (i = 0; i < number_flips; i++) {
  330.                 sumrd = 0;
  331.                 for (y = 0; y < patchsize; y++)
  332.                 for (x = 0; x < patchsize; x++)
  333.                     sumrd += Domain[y][x]*Range[i][y][x];
  334.                 fsumrd = sumrd;
  335.                 if (fmagica != 0.)
  336.                     scale = (fsumrd - fsumd*fsumr/fpatchsize_sq)/fmagica;
  337.                 else scale = 0;
  338.                 if (scale*scale < max_scale*max_scale) {
  339.                     scale = (signed char) 63. * scale / max_scale;
  340.                     scale = max_scale * scale / 63.;
  341.                     offset = (long int)(fsumr - scale*fsumd)/fpatchsize_sq;
  342.                     offset = (offset>>3)<<3;
  343.                     mean_sq = (fsumr_sq + scale*(scale*fsumd_sq - 2*fsumrd + 2*offset*fsumd)
  344.                         + offset*(offset*fpatchsize_sq - 2.*fsumr)) / fpatchsize_sq;
  345.  
  346.                     if (mean_sq < best_mean_sq) {
  347.                         besti = i;  best_mean_sq = mean_sq; bestdx = dx; bestdy = dy;
  348.                         best_scale = scale;  best_offset = offset;
  349.                         }
  350.                     if (mean_sq < mean_sq_tolerance) {
  351.                         goto gotbest;
  352.                         }
  353.                     } /* end of conditional */
  354.                }
  355.  
  356.             } /* end of Domain loop */
  357.         }
  358.           
  359. gotbest:
  360.         for (i = 0; i < count; i++)
  361.                  GrPlot(DX[i]+300+(patchsize>>1),DY[i]+(patchsize>>1),384);
  362.  
  363.         sprintf(rmsstr,"%8.4f",sqrt((double)best_mean_sq));
  364.         GrTextXY(560,20,rmsstr,255,0);
  365.         ifs[0].dx = bestdx;
  366.         ifs[0].dy = bestdy;
  367.         ifs[0].flip = besti;
  368.         ifs[0].scale = 63. * best_scale / max_scale;
  369.         ifs[0].offset = (patchsize == min_patch||mean_sq < mean_sq_tolerance) ? best_offset>>3 : -500>>3;
  370.         ifs_table[inlevel][ry/patchsize][rx/patchsize] = ifs[0];
  371.  
  372. cleanup:
  373.         GrLine(rx, ry, rx, ry + patchsize, 384);
  374.         GrLine(rx, ry ,rx + patchsize, ry, 384);
  375.         GrLine(rx + patchsize, ry, rx + patchsize, ry + patchsize, 384);
  376.         GrLine(rx, ry + patchsize, rx + patchsize, ry + patchsize, 384);
  377.  
  378.         }
  379.  
  380.         if ((out = fopen(outf, "ab")) == NULL)
  381.             {
  382.             fprintf(stderr, "Cannot open output file.\n");
  383.             return 1;
  384.             }
  385.  
  386.     if (YUV == 2) {
  387.         finish_time = time(finish_time);
  388.         time_used = (long)difftime(finish_time, start_time);
  389.         header[0].time = time_used;
  390.         fwrite(header, sizeof(struct header_t), 1, out);
  391.         }
  392.  
  393.     for (ry = 0; ry < ysize>>3; ry++)
  394.     for (rx = 0; rx < xsize>>3; rx++)  {
  395.         if (ifs_table[0][ry][rx].offset == (-500>>3))
  396.             for (y = 2*ry; y < 2*ry + 2; y++)
  397.             for (x = 2*rx; x < 2*rx + 2; x++) {
  398.                 ifs[0] = ifs_table[1][y][x];
  399.                 ifs[0].size = 1;
  400.                 fwrite(ifs, sizeof(struct ifs_t), 1, out);
  401.             }
  402.         else {
  403.             ifs[0] = ifs_table[0][ry][rx];
  404.             ifs[0].size = 0;
  405.             fwrite(ifs, sizeof(struct ifs_t), 1, out);
  406.             }
  407.         }
  408.             
  409.    fclose(out);
  410.  }  
  411.    GrSetMode(GR_default_text);
  412. /* All done. Whew... */
  413.    return 0;
  414.  
  415. }  /* end main block */
  416. #include "wdyusual.c"
  417.