home *** CD-ROM | disk | FTP | other *** search
/ Club Amiga de Montreal - CAM / CAM_CD_1.iso / files / 309.lha / PBM_PLUS / ppm / ppmquant.c < prev    next >
C/C++ Source or Header  |  1980-12-04  |  15KB  |  534 lines
  1. /* ppmquant.c - quantize the colors in a pixmap down to a specified number
  2. **
  3. ** Copyright (C) 1989 by Jef Poskanzer.
  4. **
  5. ** Permission to use, copy, modify, and distribute this software and its
  6. ** documentation for any purpose and without fee is hereby granted, provided
  7. ** that the above copyright notice appear in all copies and that both that
  8. ** copyright notice and this permission notice appear in supporting
  9. ** documentation.  This software is provided "as is" without express or
  10. ** implied warranty.
  11. */
  12.  
  13. #include <stdio.h>
  14. #ifdef    SYSV
  15. #include <string.h>
  16. #define srandom srand
  17. #define random rand
  18. #else    SYSV
  19. #include <strings.h>
  20. #endif    SYSV
  21. #include "ppm.h"
  22. #include "ppmcmap.h"
  23.  
  24. #define min(a,b) ((a) < (b) ? (a) : (b))
  25. #define max(a,b) ((a) > (b) ? (a) : (b))
  26.  
  27. #define MAXCOLORS 32768
  28. #define CLUSTER_MAXVAL 63
  29.  
  30. main( argc, argv )
  31. int argc;
  32. char *argv[];
  33.     {
  34.     FILE *ifd;
  35.     register pixel **pixels, *pP;
  36.     int argn, rows, cols, row;
  37.     register int col, limitcol;
  38.     pixval maxval;
  39.     int newcolors, colors;
  40.     register int index;
  41.     colorhist_vector chv, colormap, mediancut();
  42.     colorhash_table cht;
  43.     int floyd;
  44.     long *thisrerr, *nextrerr, *thisgerr, *nextgerr, *thisberr, *nextberr,
  45.     *temperr;
  46.     register long sr, sg, sb;
  47. #define FS_SCALE 1024
  48.     int fs_direction;
  49.     pixval err;
  50.     char *usage = "[-floyd|-fs] <newcolors> [ppmfile]";
  51.  
  52.     pm_progname = argv[0];
  53.  
  54.     argn = 1;
  55.     floyd = 0;
  56.  
  57.     if ( argn < argc && argv[argn][0] == '-' )
  58.     {
  59.     if ( strncmp(argv[argn],"-fs",max(strlen(argv[argn]),2)) == 0 ||
  60.          strncmp(argv[argn],"-floyd",max(strlen(argv[argn]),2)) == 0 )
  61.         floyd = 1;
  62.     else
  63.         pm_usage( usage );
  64.     argn++;
  65.     }
  66.  
  67.     if ( argn == argc )
  68.     pm_usage( usage );
  69.     if ( sscanf( argv[argn], "%d", &newcolors ) != 1 )
  70.     pm_usage( usage );
  71.     argn++;
  72.     if ( newcolors <= 1 )
  73.     pm_error( "number of colors must be > 1", 0,0,0,0,0 );
  74.  
  75.     if ( argn != argc )
  76.     {
  77.     ifd = pm_openr( argv[argn] );
  78.     argn++;
  79.     }
  80.     else
  81.     ifd = stdin;
  82.  
  83.     if ( argn != argc )
  84.     pm_usage( usage );
  85.  
  86.     /*
  87.     ** Step 0: read in the image.
  88.     */
  89.     pixels = ppm_readppm( ifd, &cols, &rows, &maxval );
  90.     pm_close( ifd );
  91.  
  92.  
  93.     /*
  94.     ** Step 1: attempt to make a histogram of the colors, unclustered.
  95.     */
  96.     fprintf( stderr, "(Making histogram...  " );
  97.     fflush( stderr );
  98.     chv = ppm_computecolorhist( pixels, cols, rows, MAXCOLORS, &colors );
  99.     if ( chv == (colorhist_vector) 0 )
  100.     {
  101.     fprintf( stderr, "Too many colors.)\n" );
  102.     /*
  103.     ** Step 2: try lowering maxval, to increase color coherence.
  104.     */
  105.     if ( maxval <= CLUSTER_MAXVAL )
  106.         { /* (This is not likely to happen.) */
  107.         fprintf(
  108.         stderr, "(Try recompiling with a smaller CLUSTER_MAXVAL.)\n" );
  109.         exit( 1 );
  110.         }
  111.     fprintf(
  112.         stderr,
  113.         "(Scaling colors from maxval=%d to maxval=%d to improve clustering...  ",
  114.         maxval, CLUSTER_MAXVAL );
  115.     fflush( stderr );
  116.     for ( row = 0; row < rows; row++ )
  117.         for ( col = 0, pP = pixels[row]; col < cols; col++, pP++ )
  118.         PPM_CSCALE( *pP, *pP, maxval, CLUSTER_MAXVAL );
  119.     maxval = CLUSTER_MAXVAL;
  120.     fprintf( stderr, "Done.)\n" );
  121.  
  122.     fprintf( stderr, "(Trying histogram again...  " );
  123.     fflush( stderr );
  124.     chv = ppm_computecolorhist( pixels, cols, rows, MAXCOLORS, &colors );
  125.     if ( chv == (colorhist_vector) 0 )
  126.         {
  127.         fprintf(
  128.         stderr, "Still too many colors - try recompiling with a smaller CLUSTER_MAXVAL.)\n" );
  129.         exit( 1 );
  130.         }
  131.     }
  132.     fprintf( stderr, "Done.  %d colors found.)\n", colors );
  133.  
  134.     /*
  135.     ** Step 3: apply median-cut to histogram, making the new colormap.
  136.     */
  137.     fprintf( stderr, "(Choosing %d colors...  ", newcolors );
  138.     fflush( stderr );
  139.     colormap = mediancut( chv, colors, rows * cols, maxval, newcolors );
  140.     ppm_freecolorhist( chv );
  141.     fprintf( stderr, "Done.)\n" );
  142.  
  143.     /*
  144.     ** Step 4: map the colors in the image to their closest match in the
  145.     ** new colormap, and write 'em out.
  146.     */
  147.     fprintf( stderr, "(Mapping image to new colors...  " );
  148.     fflush( stderr );
  149.     cht = ppm_alloccolorhash( );
  150.     ppm_writeppminit( stdout, cols, rows, maxval );
  151.     if ( floyd )
  152.     {
  153.     /* Initialize Floyd-Steinberg error vectors. */
  154.     thisrerr = (long *) pm_allocrow( cols + 2, sizeof(long) );
  155.     nextrerr = (long *) pm_allocrow( cols + 2, sizeof(long) );
  156.     thisgerr = (long *) pm_allocrow( cols + 2, sizeof(long) );
  157.     nextgerr = (long *) pm_allocrow( cols + 2, sizeof(long) );
  158.     thisberr = (long *) pm_allocrow( cols + 2, sizeof(long) );
  159.     nextberr = (long *) pm_allocrow( cols + 2, sizeof(long) );
  160.     srandom( (int) time( 0 ) );
  161.     for ( col = 0; col < cols + 2; col++ )
  162.         {
  163.         thisrerr[col] = random( ) % ( FS_SCALE * 2 ) - FS_SCALE;
  164.         thisgerr[col] = random( ) % ( FS_SCALE * 2 ) - FS_SCALE;
  165.         thisberr[col] = random( ) % ( FS_SCALE * 2 ) - FS_SCALE;
  166.         /* (random errors in [-1 .. 1]) */
  167.         }
  168.     fs_direction = 1;
  169.     }
  170.     for ( row = 0; row < rows; row++ )
  171.     {
  172.     if ( floyd )
  173.         for ( col = 0; col < cols + 2; col++ )
  174.         nextrerr[col] = nextgerr[col] = nextberr[col] = 0;
  175.     if ( ( ! floyd ) || fs_direction )
  176.         {
  177.         col = 0;
  178.         limitcol = cols;
  179.         pP = pixels[row];
  180.         }
  181.     else
  182.         {
  183.         col = cols - 1;
  184.         limitcol = -1;
  185.         pP = &(pixels[row][col]);
  186.         }
  187.     do
  188.         {
  189.         if ( floyd )
  190.         {
  191.         /* Use Floyd-Steinberg errors to adjust actual color. */
  192.         sr = PPM_GETR(*pP) * FS_SCALE + thisrerr[col + 1];
  193.         sg = PPM_GETG(*pP) * FS_SCALE + thisgerr[col + 1];
  194.         sb = PPM_GETB(*pP) * FS_SCALE + thisberr[col + 1];
  195.         PPM_ASSIGN( *pP, sr / FS_SCALE, sg / FS_SCALE, sb / FS_SCALE );
  196.         }
  197.  
  198.         /* Check hash table to see if we have already matched this color. */
  199.         index = ppm_lookupcolor( cht, *pP );
  200.         if ( index == -1 )
  201.         { /* No; search colormap for closest match. */
  202.         register int i, r1, g1, b1, r2, g2, b2;
  203.         register long dist, newdist;
  204.         r1 = PPM_GETR( *pP );
  205.         g1 = PPM_GETG( *pP );
  206.         b1 = PPM_GETB( *pP );
  207.         dist = 2000000000;
  208.         for ( i = 0; i < newcolors; i++ )
  209.             {
  210.             r2 = PPM_GETR( colormap[i].color );
  211.             g2 = PPM_GETG( colormap[i].color );
  212.             b2 = PPM_GETB( colormap[i].color );
  213.             newdist = ( r1 - r2 ) * ( r1 - r2 ) +
  214.                   ( g1 - g2 ) * ( g1 - g2 ) +
  215.                   ( b1 - b2 ) * ( b1 - b2 );
  216.             if ( newdist < dist )
  217.             {
  218.             index = i;
  219.             dist = newdist;
  220.             }
  221.             }
  222.         ppm_addtocolorhash( cht, *pP, index );
  223.         }
  224.  
  225.         if ( floyd )
  226.         {
  227.         /* Propagate Floyd-Steinberg error terms. */
  228.         if ( fs_direction )
  229.             {
  230.             err = sr - PPM_GETR( colormap[index].color ) * FS_SCALE;
  231.             thisrerr[col + 2] += ( err * 7 ) / 16;
  232.             nextrerr[col    ] += ( err * 3 ) / 16;
  233.             nextrerr[col + 1] += ( err * 5 ) / 16;
  234.             nextrerr[col + 2] += ( err     ) / 16;
  235.             err = sg - PPM_GETG( colormap[index].color ) * FS_SCALE;
  236.             thisgerr[col + 2] += ( err * 7 ) / 16;
  237.             nextgerr[col    ] += ( err * 3 ) / 16;
  238.             nextgerr[col + 1] += ( err * 5 ) / 16;
  239.             nextgerr[col + 2] += ( err     ) / 16;
  240.             err = sb - PPM_GETB( colormap[index].color ) * FS_SCALE;
  241.             thisberr[col + 2] += ( err * 7 ) / 16;
  242.             nextberr[col    ] += ( err * 3 ) / 16;
  243.             nextberr[col + 1] += ( err * 5 ) / 16;
  244.             nextberr[col + 2] += ( err     ) / 16;
  245.             }
  246.         else
  247.             {
  248.             err = sr - PPM_GETR( colormap[index].color ) * FS_SCALE;
  249.             thisrerr[col    ] += ( err * 7 ) / 16;
  250.             nextrerr[col + 2] += ( err * 3 ) / 16;
  251.             nextrerr[col + 1] += ( err * 5 ) / 16;
  252.             nextrerr[col    ] += ( err     ) / 16;
  253.             err = sg - PPM_GETG( colormap[index].color ) * FS_SCALE;
  254.             thisgerr[col    ] += ( err * 7 ) / 16;
  255.             nextgerr[col + 2] += ( err * 3 ) / 16;
  256.             nextgerr[col + 1] += ( err * 5 ) / 16;
  257.             nextgerr[col    ] += ( err     ) / 16;
  258.             err = sb - PPM_GETB( colormap[index].color ) * FS_SCALE;
  259.             thisberr[col    ] += ( err * 7 ) / 16;
  260.             nextberr[col + 2] += ( err * 3 ) / 16;
  261.             nextberr[col + 1] += ( err * 5 ) / 16;
  262.             nextberr[col    ] += ( err     ) / 16;
  263.             }
  264.         }
  265.  
  266.         *pP = colormap[index].color;
  267.  
  268.         if ( ( ! floyd ) || fs_direction )
  269.         {
  270.         col++;
  271.         pP++;
  272.         }
  273.         else
  274.         {
  275.         col--;
  276.         pP--;
  277.         }
  278.         }
  279.     while ( col != limitcol );
  280.  
  281.     if ( floyd )
  282.         {
  283.         temperr = thisrerr;
  284.         thisrerr = nextrerr;
  285.         nextrerr = temperr;
  286.         temperr = thisgerr;
  287.         thisgerr = nextgerr;
  288.         nextgerr = temperr;
  289.         temperr = thisberr;
  290.         thisberr = nextberr;
  291.         nextberr = temperr;
  292.         fs_direction = ! fs_direction;
  293.         }
  294.  
  295.     ppm_writeppmrow( stdout, pixels[row], cols, maxval );
  296.     }
  297.     fprintf( stderr, "Done.)\n" );
  298.  
  299.     exit( 0 );
  300.     }
  301.  
  302. /*
  303. ** Here is the fun part, the median-cut colormap generator.  This is based
  304. ** on Paul Heckbert's paper "Color Image Quantization for Frame Buffer
  305. ** Display", SIGGRAPH '82 Proceedings, page 297.
  306. */
  307.  
  308. typedef struct box *box_vector;
  309. struct box
  310.     {
  311.     int index;
  312.     int colors;
  313.     int sum;
  314.     };
  315.  
  316. colorhist_vector
  317. mediancut( chv, colors, sum, maxval, newcolors )
  318. colorhist_vector chv;
  319. int colors, sum, newcolors;
  320. pixval maxval;
  321.     {
  322.     colorhist_vector colormap;
  323.     box_vector bv;
  324.     register int bi, i;
  325.     int boxes;
  326.     int redcompare(), greencompare(), bluecompare(), sumcompare();
  327.  
  328.     bv = (box_vector) malloc( sizeof(struct box) * newcolors );
  329.     colormap =
  330.     (colorhist_vector) malloc( sizeof(struct colorhist_item) * newcolors );
  331.     if ( bv == (box_vector) 0 || colormap == (colorhist_vector) 0 )
  332.     pm_error( "out of memory", 0,0,0,0,0 );
  333.     for ( i = 0; i < newcolors; i++ )
  334.     PPM_ASSIGN( colormap[i].color, 0, 0, 0 );
  335.  
  336.     /*
  337.     ** Set up the initial box.
  338.     */
  339.     bv[0].index = 0;
  340.     bv[0].colors = colors;
  341.     bv[0].sum = sum;
  342.     boxes = 1;
  343.  
  344.     /*
  345.     ** Main loop: split boxes until we have enough.
  346.     */
  347.     while ( boxes < newcolors )
  348.     {
  349.     register int indx, clrs;
  350.     int sm;
  351.     register int minr, maxr, ming, maxg, minb, maxb, v;
  352.     int halfsum, lowersum;
  353.  
  354.     /*
  355.     ** Find the first splittable box.
  356.     */
  357.     for ( bi = 0; bv[bi].colors < 2 && bi < boxes; bi++ )
  358.         ;
  359.     if ( bi == boxes )
  360.         break;    /* ran out of colors! */
  361.     indx = bv[bi].index;
  362.     clrs = bv[bi].colors;
  363.     sm = bv[bi].sum;
  364.  
  365.     /*
  366.     ** Go through the box finding the minimum and maximum of each
  367.     ** component -- the boundaries of the box.
  368.     */
  369.     minr = maxr = PPM_GETR( chv[indx].color );
  370.     ming = maxg = PPM_GETG( chv[indx].color );
  371.     minb = maxb = PPM_GETB( chv[indx].color );
  372.     for ( i = 1; i < clrs; i++ )
  373.         {
  374.         v = PPM_GETR( chv[indx + i].color );
  375.         if ( v < minr ) minr = v;
  376.         if ( v > maxr ) maxr = v;
  377.         v = PPM_GETG( chv[indx + i].color );
  378.         if ( v < ming ) ming = v;
  379.         if ( v > maxg ) maxg = v;
  380.         v = PPM_GETB( chv[indx + i].color );
  381.         if ( v < minb ) minb = v;
  382.         if ( v > maxb ) maxb = v;
  383.         }
  384.  
  385.     /*
  386.     ** Find the largest dimension, and sort by that component.
  387.     */
  388.     if ( maxr - minr >= maxg - ming && maxr - minr >= maxb - minb )
  389.         qsort(
  390.         (char *) &(chv[indx]), clrs, sizeof(struct colorhist_item),
  391.         redcompare );
  392.     else if ( maxg - ming >= maxb - minb )
  393.         qsort(
  394.         (char *) &(chv[indx]), clrs, sizeof(struct colorhist_item),
  395.         greencompare );
  396.     else
  397.         qsort(
  398.         (char *) &(chv[indx]), clrs, sizeof(struct colorhist_item),
  399.         bluecompare );
  400.     
  401.     /*
  402.     ** Now find the median based on the counts, so that about half the
  403.     ** pixels (not colors, pixels) are in each subdivision.
  404.     */
  405.     lowersum = chv[indx].value;
  406.     halfsum = sm / 2;
  407.     for ( i = 1; i < clrs - 1; i++ )
  408.         {
  409.         if ( lowersum >= halfsum )
  410.         break;
  411.         lowersum += chv[indx + i].value;
  412.         }
  413.  
  414.     /*
  415.     ** Split the box, and sort to bring the biggest box to the top.
  416.     */
  417.     bv[bi].colors = i;
  418.     bv[bi].sum = lowersum;
  419.     bv[boxes].index = indx + i;
  420.     bv[boxes].colors = clrs - i;
  421.     bv[boxes].sum = sm - lowersum;
  422.     boxes++;
  423.     qsort( (char *) bv, boxes, sizeof(struct box), sumcompare );
  424.     }
  425.  
  426.     /*
  427.     ** Ok, we've got enough boxes.  Now choose a representative color for
  428.     ** each box.  There are a number of possible ways to make this choice.
  429.     ** One would be to choose the center of the box; this ignores any structure
  430.     ** within the boxes.  Another method would be to average all the colors in
  431.     ** the box -- this is the method specified in Heckbert's paper.  A third
  432.     ** method is to average all the pixels in the box.
  433.     */
  434. /* #define CENTER_BOX */
  435. /* #define AVERAGE_COLORS */
  436. #define AVERAGE_PIXELS
  437.     for ( bi = 0; bi < boxes; bi++ )
  438.     {
  439. #ifdef CENTER_BOX
  440.     register int indx = bv[bi].index;
  441.     register int clrs = bv[bi].colors;
  442.     register int minr, maxr, ming, maxg, minb, maxb, v;
  443.  
  444.     minr = maxr = PPM_GETR( chv[indx].color );
  445.     ming = maxg = PPM_GETG( chv[indx].color );
  446.     minb = maxb = PPM_GETB( chv[indx].color );
  447.     for ( i = 1; i < clrs; i++ )
  448.         {
  449.         v = PPM_GETR( chv[indx + i].color );
  450.         minr = min( minr, v );
  451.         maxr = max( maxr, v );
  452.         v = PPM_GETG( chv[indx + i].color );
  453.         ming = min( ming, v );
  454.         maxg = max( maxg, v );
  455.         v = PPM_GETB( chv[indx + i].color );
  456.         minb = min( minb, v );
  457.         maxb = max( maxb, v );
  458.         }
  459.     PPM_ASSIGN(
  460.         colormap[bi].color, ( minr + maxr ) / 2, ( ming + maxg ) / 2,
  461.         ( minb + maxb ) / 2 );
  462. #endif CENTER_BOX
  463. #ifdef AVERAGE_COLORS
  464.     register int indx = bv[bi].index;
  465.     register int clrs = bv[bi].colors;
  466.     register long r = 0, g = 0, b = 0;
  467.  
  468.     for ( i = 0; i < clrs; i++ )
  469.         {
  470.         r += PPM_GETR( chv[indx + i].color );
  471.         g += PPM_GETG( chv[indx + i].color );
  472.         b += PPM_GETB( chv[indx + i].color );
  473.         }
  474.     r = r / clrs;
  475.     g = g / clrs;
  476.     b = b / clrs;
  477.     PPM_ASSIGN( colormap[bi].color, r, g, b );
  478. #endif AVERAGE_COLORS
  479. #ifdef AVERAGE_PIXELS
  480.     register int indx = bv[bi].index;
  481.     register int clrs = bv[bi].colors;
  482.     register long r = 0, g = 0, b = 0, sum = 0;
  483.  
  484.     for ( i = 0; i < clrs; i++ )
  485.         {
  486.         r += PPM_GETR( chv[indx + i].color ) * chv[indx + i].value;
  487.         g += PPM_GETG( chv[indx + i].color ) * chv[indx + i].value;
  488.         b += PPM_GETB( chv[indx + i].color ) * chv[indx + i].value;
  489.         sum += chv[indx + i].value;
  490.         }
  491.     r = r / sum;
  492.     if ( r > maxval ) r = maxval;    /* avoid math errors */
  493.     g = g / sum;
  494.     if ( g > maxval ) g = maxval;
  495.     b = b / sum;
  496.     if ( b > maxval ) b = maxval;
  497.     PPM_ASSIGN( colormap[bi].color, r, g, b );
  498. #endif AVERAGE_PIXELS
  499.     }
  500.  
  501.     /*
  502.     ** All done.
  503.     */
  504.     return colormap;
  505.     }
  506.  
  507. int
  508. redcompare( ch1, ch2 )
  509. colorhist_vector ch1, ch2;
  510.     {
  511.     return (int) PPM_GETR( ch1->color ) - (int) PPM_GETR( ch2->color );
  512.     }
  513.  
  514. int
  515. greencompare( ch1, ch2 )
  516. colorhist_vector ch1, ch2;
  517.     {
  518.     return (int) PPM_GETG( ch1->color ) - (int) PPM_GETG( ch2->color );
  519.     }
  520.  
  521. int
  522. bluecompare( ch1, ch2 )
  523. colorhist_vector ch1, ch2;
  524.     {
  525.     return (int) PPM_GETB( ch1->color ) - (int) PPM_GETB( ch2->color );
  526.     }
  527.  
  528. int
  529. sumcompare( b1, b2 )
  530. box_vector b1, b2;
  531.     {
  532.     return b2->sum - b1->sum;
  533.     }
  534.