home *** CD-ROM | disk | FTP | other *** search
/ Photo CD Demo 1 / Demo.bin / formats / gif / gen_code / median.c < prev    next >
Encoding:
Internet Message Format  |  1988-05-26  |  28.4 KB
  1. From: falk@sun.uucp (Ed Falk)
  2. Subject: Re: 24bitcolor --> 8bit colormap
  3. Date: 26 Mar 88 02:12:04 GMT
  4.  
  5. In article <221@abvax.UUCP>, gfs@abvax.UUCP (Greg F. Shay) writes:
  6. > 
  7. > Relating to the question of how to choose a 256 entry colormap to best
  8. > represent a 24bit full color image, the algorithm to use was published
  9. > by Paul Heckbert in the following:
  10. > 
  11. >     "Color image quantization for frame buffer display," SIGGRAPH 1982
  12. > Proceedings, pp.297-307.
  13. > 
  14.  
  15. OK, since everybody else wants to show how they do it, I'll submit
  16. my own work.
  17.  
  18. This is my implementation of the median cut algorithm.  It does all of
  19. the optimizations suggested by Heckbert in his paper.  It will dither
  20. a 1152x900x24  image in 1.5 minutes.  Floyd-Steinberg dithering is
  21. available as an option.   You can choose the size of your target color-table
  22. to be anything you want.
  23.  
  24. This should not be too hard to port to other image formats besides the
  25. sun rasterfile format.
  26.  
  27. I have made a couple of improvements over the original algorithm.  In this
  28. version, the subdivision is always performed based on the largest box
  29. rather than arbitrarily dividing them all down.
  30.  
  31. Issues:  these are some improvements I've been meaning to make, but haven't
  32. gotten around to it:
  33.  
  34. Currently, the largest box is chosen based purely on its
  35. volume r*g*b.  This should be biased according to perceptual rules.
  36.  
  37.     volume = r*.3 + g*.59 + b*.11
  38.  
  39. this may produce a better image.
  40.  
  41. Also, gamma-correction should be applied on the output, and inversely on the
  42. input.
  43.  
  44. Dithering remains a problem.  Dithering only works if you have two values
  45. to dither between.  For instance, if you're dithering a blue field, and you
  46. only have one shade of blue in the color table, and the nearest other color
  47. is white, then a slightly lighter area of blue will produce blue with
  48. undesirable white specs in the middle of it.  We need a method of either
  49. making sure the color table contains reasonable bracketing or of making
  50. sure that the accumulated error in F-S dithering does not go beyond a
  51. certain limit.
  52.  
  53. Many images look better *without* dithering, try it both ways.
  54.  
  55.     -ed falk, falk@sun
  56.  
  57.  
  58. ==================================
  59. /* quantize stuff:
  60.  *
  61.  * median [-csize n] [-fsd] [ifile [ofile]]
  62.  *
  63.  * -csize n    - set colortable size.  Default is 256.
  64.  * -fsd        - use Floyd-Steinberg dithering.
  65.  *
  66.  */
  67.  
  68.  
  69. #include <stdio.h>
  70. #include <sys/types.h>
  71. #include <pixrect/pixrect_hs.h>
  72.  
  73. #define    MAX_CMAP_SIZE    256
  74.  
  75. #define    COLOR_DEPTH    8
  76. #define    MAX_COLOR    256
  77.  
  78. #define    B_DEPTH        5        /* # bits/pixel to use */
  79. #define    B_LEN        (1<<B_DEPTH)
  80.  
  81. #define    C_DEPTH        2
  82. #define    C_LEN        (1<<C_DEPTH)    /* # cells/color to use */
  83.  
  84.  
  85.  
  86. typedef    struct colorbox {
  87.       struct colorbox    *next, *prev ;
  88.       int            rmin,rmax, gmin,gmax, bmin,bmax ;
  89.       int            total ;
  90.     } Colorbox ;
  91.  
  92. typedef struct {
  93.       int    num_ents ;
  94.       int    entries[MAX_CMAP_SIZE][2] ;
  95.     } C_cell ;
  96.  
  97. static    struct rasterfile    rh;
  98. static    colormap_t        colormap, ocmap ;
  99. static    unsigned char        rm[MAX_CMAP_SIZE],
  100.                 gm[MAX_CMAP_SIZE],
  101.                 bm[MAX_CMAP_SIZE] ;
  102. static    int            bytes_per_pixel ;
  103. static    int            num_colors ;
  104.  
  105. static    int            histogram[B_LEN][B_LEN][B_LEN] ;
  106.  
  107. static    Colorbox        *freeboxes ;
  108. static    Colorbox        *usedboxes ;
  109. static    Colorbox        *largest_box() ;
  110.  
  111. static    C_cell            **ColorCells ;
  112.  
  113.  
  114. static    char    *usage =
  115. "usage: median [-csize n] [-fsd] [ifile [ofile]]\n" ;
  116.  
  117.  
  118. main(argc, argv)
  119.     int argc;
  120.     char **argv;
  121. {
  122. static    struct rasterfile    outheader ;
  123.     int            i, j ;
  124.     int            dither = 0 ;
  125.     char            *ifile_name = NULL, *ofile_name = NULL ;
  126.     Colorbox        *box_list, *ptr ;
  127.  
  128.     num_colors = MAX_CMAP_SIZE ;
  129.  
  130.     while(--argc)
  131.     {
  132.       ++argv ;
  133.       if(strcmp(*argv,"-fsd") == 0)
  134.         dither = 1 ;
  135.  
  136.       else if(strcmp(*argv,"-csize") == 0)
  137.       {
  138.         if(argc < 2)
  139.         {
  140.           fprintf(stderr,"-csize requires 2 arguments, %s",usage) ;
  141.           exit(1) ;
  142.         }
  143.         else
  144.         {
  145.           argc -= 1 ;
  146.           num_colors = atoi(*++argv) ;
  147.           if(num_colors > MAX_CMAP_SIZE)
  148.           {
  149.         fprintf(stderr,"-csize specifies colormap greater than %d\n%s",
  150.             MAX_CMAP_SIZE,usage) ;
  151.         exit(1) ;
  152.           }
  153.         }
  154.       }
  155.  
  156.       else if(**argv == '-')
  157.       {
  158.         fprintf(stderr,"unknown argument '%s'\n%s",*argv,usage) ;
  159.         exit(1) ;
  160.       }
  161.  
  162.       else if(!ifile_name)
  163.       {
  164.         ifile_name = *argv ;
  165.       }
  166.  
  167.       else if(!ofile_name)
  168.       {
  169.         ofile_name = *argv ;
  170.       }
  171.  
  172.       else
  173.       {
  174.         fprintf(stderr, "Too many file names\n%s",usage) ;
  175.         exit(1) ;
  176.       }
  177.     }
  178.  
  179.     if(ifile_name)
  180.       if (freopen(ifile_name, "r", stdin) == NULL)
  181.       {
  182.         fprintf(stderr,"Cannot open input file %s\n%s", ifile_name,usage);
  183.         exit(1);
  184.       }
  185.  
  186.     if(ofile_name)
  187.       if (freopen(ofile_name, "w", stdout) == NULL)
  188.       {
  189.         fprintf(stderr,"Cannot open output file %s\n%s", ofile_name,usage);
  190.         exit(1);
  191.       }
  192.  
  193.  
  194.     if (pr_load_header(stdin, &rh) != 0  ||
  195.         pr_load_colormap(stdin, &rh, &colormap) != 0 )
  196.     {
  197.       fprintf(stderr,"Error reading rasterfile header.\n");
  198.       exit(1);
  199.     }
  200.  
  201.  
  202.  
  203.     switch( rh.ras_depth )
  204.     {
  205.       case 24: bytes_per_pixel = 3 ; break ;
  206.       case 32: bytes_per_pixel = 4 ; break ;
  207.       default:
  208.         fprintf(stderr,
  209.             "Image is %d bits deep, median only works with 24 or 32\n");
  210.         exit(1) ;
  211.     }
  212.  
  213.  
  214.     /*
  215.      * STEP 1:  create empty boxes
  216.      */
  217.  
  218.     usedboxes = NULL ;
  219.     box_list = freeboxes =
  220.         (Colorbox *) malloc( num_colors * sizeof(Colorbox) ) ;
  221.  
  222.     freeboxes[0].next = &freeboxes[1] ;
  223.     freeboxes[0].prev = NULL ;
  224.     for(i=1; i<num_colors-1; ++i)
  225.     {
  226.       freeboxes[i].next = &freeboxes[i+1] ;
  227.       freeboxes[i].prev = &freeboxes[i-1] ;
  228.     }
  229.     freeboxes[num_colors-1].next = NULL ;
  230.     freeboxes[num_colors-1].prev = &freeboxes[num_colors-2] ;
  231.  
  232.  
  233.  
  234.     /*
  235.      * STEP 2: get histogram, initialize first box
  236.      */
  237.  
  238.     ptr = freeboxes ;
  239.     freeboxes = ptr->next ;
  240.     if(freeboxes) freeboxes->prev = NULL ;
  241.  
  242.     ptr->next = usedboxes ;
  243.     usedboxes = ptr ;
  244.     if(ptr->next) ptr->next->prev = ptr ;
  245.     
  246.     get_histogram(ptr) ;
  247.  
  248.     if( fseek(stdin, 0L, 0) != 0)
  249.     {
  250.       fprintf(stderr, "median: input must be from file\n") ;
  251.       exit(1) ;
  252.     }
  253.  
  254.  
  255.  
  256.     /*
  257.      * STEP 3: continually subdivide boxes until no more free
  258.      * boxes remain
  259.      */
  260.  
  261.     while(freeboxes != NULL)
  262.     {
  263.       ptr = largest_box() ;
  264.       splitbox(ptr) ;
  265.     }
  266.  
  267.  
  268.  
  269.  
  270.  
  271.  
  272.     /*
  273.      * STEP 4: assign colors to all boxes
  274.      */
  275.  
  276.     for(i=0, ptr=usedboxes; i<num_colors; ++i, ptr=ptr->next)
  277.       assign_color(ptr,&rm[i],&gm[i],&bm[i]) ;
  278.  
  279.     /* We're done with the boxes now */
  280.  
  281.     free(box_list) ;
  282.     box_list = freeboxes = usedboxes = NULL ;
  283.  
  284.  
  285.  
  286.  
  287.     /*
  288.      * STEP 5: scan histogram and map all values to closest color
  289.      */
  290.  
  291.     /* 5a: create cell list as described in Heckbert[2] */
  292.  
  293.     {
  294.       register C_cell    **ptr ;
  295.       register int        n = C_LEN*C_LEN*C_LEN ;
  296.  
  297.       ptr = ColorCells =
  298.         (C_cell **) malloc(C_LEN*C_LEN*C_LEN * sizeof(C_cell *));
  299.  
  300.       while( n-- > 0 )
  301.         *ptr++ = NULL ;
  302.     }
  303.  
  304.  
  305.  
  306.  
  307.     /* 5b: create mapping from truncated pixel space to color
  308.        table entries */
  309.  
  310.     map_colortable() ;
  311.  
  312.  
  313.  
  314.  
  315.  
  316.     /*
  317.      * STEP 6: scan image, match input values to table entries
  318.      */
  319.  
  320.     pr_load_header(stdin, &rh) ;
  321.     pr_load_colormap(stdin, &rh, &colormap) ;
  322.     pr_read_init(&rh) ;
  323.  
  324.     outheader.ras_magic = RAS_MAGIC ;
  325.     outheader.ras_width = rh.ras_width ;
  326.     outheader.ras_height = rh.ras_height ;
  327.     outheader.ras_depth = COLOR_DEPTH ;
  328.     outheader.ras_length =
  329.         ((outheader.ras_width+1) & ~1) * outheader.ras_height ;
  330.     outheader.ras_type = RT_STANDARD ;
  331.     outheader.ras_maptype = RMT_EQUAL_RGB ;
  332.     outheader.ras_maplength = num_colors*3 ;
  333.  
  334.     ocmap.type = RMT_EQUAL_RGB ;
  335.     ocmap.length = num_colors ;
  336.     ocmap.map[0] = rm ;
  337.     ocmap.map[1] = gm ;
  338.     ocmap.map[2] = bm ;
  339.  
  340.     pr_dump_header(stdout, &outheader, &ocmap) ;
  341.  
  342.     if(dither)
  343.       quant_fsdither() ;
  344.     else
  345.       quant() ;
  346.  
  347.     fclose(stdout) ;
  348. }
  349.  
  350.  
  351.  
  352.  
  353.  
  354. static
  355. get_histogram(box)
  356. register Colorbox    *box ;
  357. {
  358.     unsigned char    *inline ;
  359. register unsigned char    *inptr ;
  360.     int        i ;
  361. register int        j ;
  362.  
  363.  
  364.     pr_read_init(&rh) ;
  365.  
  366.     inline = (unsigned char *) alloca( rh.ras_width * bytes_per_pixel ) ;
  367.  
  368.     box->rmin = box->gmin = box->bmin = 999 ;
  369.     box->rmax = box->gmax = box->bmax = -1 ;
  370.     box->total = rh.ras_height*rh.ras_width ;
  371.  
  372.     {
  373.       register int *ptr = &histogram[0][0][0] ;
  374.       for(i=B_LEN*B_LEN*B_LEN; --i >= 0;)
  375.         *ptr++ = 0 ;
  376.     }
  377.  
  378.  
  379.     for(i=rh.ras_height; --i >= 0;)
  380.     {
  381.       pr_get_bytes(stdin, &rh, rh.ras_width * bytes_per_pixel, inline) ;
  382.       inptr = inline ;
  383.       for(j=rh.ras_width; --j >= 0;)
  384.       {
  385.         register int    red,green,blue ;
  386.  
  387.         if(rh.ras_depth == 32)
  388.           ++inptr ;
  389.         blue = *inptr++ ;
  390.         green = *inptr++ ;
  391.         red = *inptr++ ;
  392.  
  393.         if (rh.ras_maplength > 0)
  394.         {
  395.           red = *(colormap.map[0]+red) ;
  396.           green = *(colormap.map[1]+green) ;
  397.           blue = *(colormap.map[2]+blue) ;
  398.         }
  399.  
  400.         red >>= (COLOR_DEPTH-B_DEPTH) ;
  401.         green >>= (COLOR_DEPTH-B_DEPTH) ;
  402.         blue >>= (COLOR_DEPTH-B_DEPTH) ;
  403.  
  404.         if(red < box->rmin) box->rmin = red ;
  405.         if(red > box->rmax) box->rmax = red ;
  406.         if(green < box->gmin) box->gmin = green ;
  407.         if(green > box->gmax) box->gmax = green ;
  408.         if(blue < box->bmin) box->bmin = blue ;
  409.         if(blue > box->bmax) box->bmax = blue ;
  410.  
  411.         ++histogram[red][green][blue] ;
  412.       }
  413.     }
  414. }
  415.  
  416.  
  417.  
  418. static Colorbox *
  419. largest_box()
  420. {
  421. register Colorbox    *tmp = usedboxes, *ptr ;
  422. register int         size = -1 ;
  423.  
  424.     while(tmp != NULL)
  425.     {
  426.       if( (tmp->rmax > tmp->rmin  ||
  427.            tmp->gmax > tmp->gmin  ||
  428.            tmp->bmax > tmp->bmin)  &&  tmp->total > size )
  429.       {
  430.         ptr = tmp ;
  431.         size = tmp->total ;
  432.       }
  433.       tmp = tmp->next ;
  434.     }
  435.     return(ptr) ;
  436. }
  437.  
  438.  
  439.  
  440. static
  441. splitbox(ptr)
  442. register Colorbox    *ptr ;
  443. {
  444.     int        hist2[B_LEN] ;
  445.     int        first, last ;
  446. register Colorbox    *new ;
  447. register int        *iptr, *histp ;
  448. register int        i, j ;
  449.     enum        {RED,GREEN,BLUE} which ;
  450.  
  451.     /*
  452.      * see which axis is the largest, do a histogram along that
  453.      * axis.  Split at median point.  Contract both new boxes to
  454.      * fit points and return
  455.      */
  456.  
  457.     if( (i = ptr->rmax - ptr->rmin) >= (ptr->gmax - ptr->gmin)  &&
  458.         i >= (ptr->bmax - ptr->bmin) )
  459.       which = RED ;
  460.     else if( (ptr->gmax - ptr->gmin) >= (ptr->bmax - ptr->bmin) )
  461.       which = GREEN ;
  462.     else
  463.       which = BLUE ;
  464.  
  465.  
  466.  
  467.     /* get histogram along longest axis */
  468.     {
  469.       register int    ir,ig,ib ;
  470.       switch(which)
  471.       {
  472.         case RED:
  473.           histp = &hist2[ptr->rmin] ;
  474.           for(ir = ptr->rmin; ir <= ptr->rmax; ++ir)
  475.           {
  476.         *histp = 0 ;
  477.         for(ig = ptr->gmin; ig <= ptr->gmax; ++ig)
  478.         {
  479.           iptr = &histogram[ir][ig][ptr->bmin] ;
  480.           for(ib = ptr->bmin; ib <= ptr->bmax; ++ib)
  481.           {
  482.             *histp += *iptr ;
  483.             ++iptr ;
  484.           }
  485.         }
  486.         ++histp ;
  487.           }
  488.           first = ptr->rmin;  last = ptr->rmax ;
  489.           break ;
  490.  
  491.         case GREEN:
  492.           histp = &hist2[ptr->gmin] ;
  493.           for(ig = ptr->gmin; ig <= ptr->gmax; ++ig)
  494.           {
  495.         *histp = 0 ;
  496.         for(ir = ptr->rmin; ir <= ptr->rmax; ++ir)
  497.         {
  498.           iptr = &histogram[ir][ig][ptr->bmin] ;
  499.           for(ib = ptr->bmin; ib <= ptr->bmax; ++ib)
  500.           {
  501.             *histp += *iptr ;
  502.             ++iptr ;
  503.           }
  504.         }
  505.         ++histp ;
  506.           }
  507.           first = ptr->gmin;  last = ptr->gmax ;
  508.           break ;
  509.  
  510.         case BLUE:
  511.           histp = &hist2[ptr->bmin] ;
  512.           for(ib = ptr->bmin; ib <= ptr->bmax; ++ib)
  513.           {
  514.         *histp = 0 ;
  515.         for(ir = ptr->rmin; ir <= ptr->rmax; ++ir)
  516.         {
  517.           iptr = &histogram[ir][ptr->gmin][ib] ;
  518.           for(ig = ptr->gmin; ig <= ptr->gmax; ++ig)
  519.           {
  520.             *histp += *iptr ;
  521.             iptr += B_LEN ;
  522.           }
  523.         }
  524.         ++histp ;
  525.           }
  526.           first = ptr->bmin;  last = ptr->bmax ;
  527.           break ;
  528.       }
  529.     }
  530.  
  531.     /* find median point */
  532.     {
  533.       register int sum, sum2 ;
  534.       histp = &hist2[first] ;
  535.  
  536.       sum2 = ptr->total/2 ;
  537.       histp = &hist2[first] ;
  538.       sum = 0 ;
  539.       for( i=first; i<=last && (sum += *histp++) < sum2; ++i) ;
  540.       if( i == first ) ++i ;
  541.     }
  542.  
  543.     /* Create new box, re-allocate points */
  544.  
  545.     new = freeboxes ;
  546.     freeboxes = new->next ;
  547.     if(freeboxes) freeboxes->prev = NULL ;
  548.  
  549.     if(usedboxes) usedboxes->prev = new ;
  550.     new->next = usedboxes ;
  551.     usedboxes = new ;
  552.  
  553.     {
  554.       register int    sum1,sum2,j ;
  555.       int        total ;
  556.       total = ptr->total ;
  557.       histp = &hist2[first] ;
  558.       sum1 = 0 ;
  559.       for( j = first; j < i; ++j)
  560.         sum1 += *histp++ ;
  561.       sum2 = 0 ;
  562.       for( j = i; j <= last; ++j)
  563.         sum2 += *histp++ ;
  564. #ifdef    DEBUG
  565.       if( sum1+sum2 != total  ||  sum1 == 0  ||  sum2 == 0 )
  566.         fprintf(stderr, "??? splitbox: sum1=%d, sum2=%d, total=%d\n",
  567.         sum1,sum2,total) ;
  568. #endif    DEBUG
  569.       new->total = sum1 ;
  570.       ptr->total = sum2 ;
  571.     }
  572.  
  573.  
  574.     new->rmin = ptr->rmin ;
  575.     new->rmax = ptr->rmax ;
  576.     new->gmin = ptr->gmin ;
  577.     new->gmax = ptr->gmax ;
  578.     new->bmin = ptr->bmin ;
  579.     new->bmax = ptr->bmax ;
  580.     switch(which)
  581.     {
  582.       case RED:
  583.         new->rmax = i-1 ;
  584.         ptr->rmin = i ;
  585.         break ;
  586.  
  587.       case GREEN:
  588.         new->gmax = i-1 ;
  589.         ptr->gmin = i ;
  590.       break ;
  591.  
  592.       case BLUE:
  593.         new->bmax = i-1 ;
  594.         ptr->bmin = i ;
  595.       break ;
  596.     }
  597.     shrinkbox(new) ;
  598.     shrinkbox(ptr) ;
  599. }
  600.  
  601.  
  602.  
  603.  
  604. static
  605. shrinkbox(box)
  606. register Colorbox    *box ;
  607. {
  608. register int        *histp ;
  609. register int        ir,ig,ib ;
  610.  
  611.  
  612.     if(box->rmax > box->rmin)
  613.     {
  614.       for(ir = box->rmin; ir <= box->rmax; ++ir)
  615.         for(ig = box->gmin; ig <= box->gmax; ++ig)
  616.         {
  617.           histp = &histogram[ir][ig][box->bmin] ;
  618.           for(ib = box->bmin; ib <= box->bmax; ++ib)
  619.         if( *histp++ != 0 )
  620.         {
  621.           box->rmin = ir ;
  622.           goto have_rmin ;
  623.         }
  624.         }
  625. have_rmin:
  626.  
  627.       if(box->rmax > box->rmin)
  628.         for(ir = box->rmax; ir >= box->rmin; --ir)
  629.           for(ig = box->gmin; ig <= box->gmax; ++ig)
  630.           {
  631.         histp = &histogram[ir][ig][box->bmin] ;
  632.         for(ib = box->bmin; ib <= box->bmax; ++ib)
  633.           if( *histp++ != 0 )
  634.           {
  635.             box->rmax = ir ;
  636.             goto have_rmax ;
  637.           }
  638.           }
  639.     }
  640. have_rmax:
  641.  
  642.     if( box->gmax > box->gmin )
  643.     {
  644.       for(ig = box->gmin; ig <= box->gmax; ++ig)
  645.         for(ir = box->rmin; ir <= box->rmax; ++ir)
  646.         {
  647.           histp = &histogram[ir][ig][box->bmin] ;
  648.           for(ib = box->bmin; ib <= box->bmax; ++ib)
  649.         if( *histp++ != 0 )
  650.         {
  651.           box->gmin = ig ;
  652.           goto have_gmin ;
  653.         }
  654.         }
  655. have_gmin:
  656.  
  657.       if( box->gmax > box->gmin )
  658.         for(ig = box->gmax; ig >= box->gmin; --ig)
  659.           for(ir = box->rmin; ir <= box->rmax; ++ir)
  660.           {
  661.         histp = &histogram[ir][ig][box->bmin] ;
  662.         for(ib = box->bmin; ib <= box->bmax; ++ib)
  663.           if( *histp++ != 0 )
  664.           {
  665.             box->gmax = ig ;
  666.             goto have_gmax ;
  667.           }
  668.           }
  669.     }
  670. have_gmax:
  671.  
  672.     if( box->bmax > box->bmin )
  673.     {
  674.       for(ib = box->bmin; ib <= box->bmax; ++ib)
  675.         for(ir = box->rmin; ir <= box->rmax; ++ir)
  676.         {
  677.           histp = &histogram[ir][box->gmin][ib] ;
  678.           for(ig = box->gmin; ig <= box->gmax; ++ig)
  679.           {
  680.         if( *histp != 0 )
  681.         {
  682.           box->bmin = ib ;
  683.           goto have_bmin ;
  684.         }
  685.         histp += B_LEN ;
  686.           }
  687.         }
  688. have_bmin:
  689.  
  690.       if( box->bmax > box->bmin )
  691.         for(ib = box->bmax; ib >= box->bmin; --ib)
  692.           for(ir = box->rmin; ir <= box->rmax; ++ir)
  693.           {
  694.         histp = &histogram[ir][box->gmin][ib] ;
  695.         for(ig = box->gmin; ig <= box->gmax; ++ig)
  696.         {
  697.           if( *histp != 0 )
  698.           {
  699.             box->bmax = ib ;
  700.             goto have_bmax ;
  701.           }
  702.           histp += B_LEN ;
  703.         }
  704.           }
  705.     }
  706. have_bmax: return ;
  707. }
  708.  
  709.  
  710.  
  711.  
  712. static
  713. assign_color(ptr,rp,gp,bp)
  714. register Colorbox    *ptr ;
  715.     unsigned char    *rp,*gp,*bp ;
  716. {
  717. register int        *histp ;
  718. register int        ir,ig,ib ;
  719. register long        red = 0, green = 0, blue = 0 ;
  720.  
  721.     *rp = ((ptr->rmin + ptr->rmax) << (COLOR_DEPTH - B_DEPTH)) / 2 ;
  722.     *gp = ((ptr->gmin + ptr->gmax) << (COLOR_DEPTH - B_DEPTH)) / 2 ;
  723.     *bp = ((ptr->bmin + ptr->bmax) << (COLOR_DEPTH - B_DEPTH)) / 2 ;
  724.  
  725. #ifdef    COMMENT
  726. #ifdef    DEBUG
  727.     if( ptr->total <= 0 )
  728.     {
  729.       fprintf(stderr,"??? assign_color: total = %d\n",ptr->total) ;
  730.       return ;
  731.     }
  732. #endif    DEBUG
  733.  
  734.     for( ir = ptr->rmin; ir <= ptr->rmax; ++ir)
  735.       for( ig = ptr->gmin; ig <= ptr->gmax; ++ig)
  736.       {
  737.         histp = &histogram[ir][ig][ptr->bmin] ;
  738.         for( ib = ptr->bmin; ib <= ptr->bmax; ++ib)
  739.         {
  740.           if( *histp != 0 )
  741.           {
  742.         red += ir * *histp ;
  743.         green += ig * *histp ;
  744.         blue += ib * *histp ;
  745.           }
  746.           ++histp ;
  747.         }
  748.       }
  749.  
  750.     *rp = (red << (COLOR_DEPTH - B_DEPTH)) / ptr->total ;
  751.     *gp = (green << (COLOR_DEPTH - B_DEPTH)) / ptr->total ;
  752.     *bp = (blue << (COLOR_DEPTH - B_DEPTH)) / ptr->total ;
  753. #endif    COMMENT
  754. }
  755.  
  756.  
  757.  
  758. static    C_cell *
  759. create_colorcell(red,green,blue)
  760.     int        red,green,blue ;
  761. {
  762. register int        ir,ig,ib, i ;
  763.     int        mindist ;
  764. register C_cell        *ptr ;
  765.  
  766.     ir = red >> (COLOR_DEPTH-C_DEPTH) ;
  767.     ig = green >> (COLOR_DEPTH-C_DEPTH) ;
  768.     ib = blue >> (COLOR_DEPTH-C_DEPTH) ;
  769.  
  770.     red &= ~1 << (COLOR_DEPTH-C_DEPTH) ;
  771.     green &= ~1 << (COLOR_DEPTH-C_DEPTH) ;
  772.     blue &= ~1 << (COLOR_DEPTH-C_DEPTH) ;
  773.  
  774.     ptr = (C_cell *) malloc( sizeof(C_cell) ) ;
  775.     *( ColorCells + ir*C_LEN*C_LEN + ig*C_LEN + ib ) = ptr ;
  776.  
  777.     ptr->num_ents = 0 ;
  778.  
  779.  
  780.     /* step 1: find all colors inside this cell, while we're at
  781.           it, find distance of centermost point to furthest
  782.           corner */
  783.  
  784.     mindist = 99999999 ;
  785.     for(i=0; i<num_colors; ++i)
  786.       if( rm[i]>>(COLOR_DEPTH-C_DEPTH) == ir  &&
  787.           gm[i]>>(COLOR_DEPTH-C_DEPTH) == ig  &&
  788.           bm[i]>>(COLOR_DEPTH-C_DEPTH) == ib)
  789.       {
  790.         register int    tmp, dist ;
  791.  
  792.         ptr->entries[ptr->num_ents][0] = i ;
  793.         ptr->entries[ptr->num_ents][1] = 0 ;
  794.         ++ptr->num_ents ;
  795.  
  796.         tmp = rm[i] - red ;
  797.         if( tmp < (MAX_COLOR/C_LEN/2) )
  798.           tmp = MAX_COLOR/C_LEN-1 - tmp ;
  799.         dist = tmp*tmp ;
  800.  
  801.         tmp = gm[i] - green ;
  802.         if( tmp < (MAX_COLOR/C_LEN/2) )
  803.           tmp = MAX_COLOR/C_LEN-1 - tmp ;
  804.         dist += tmp*tmp ;
  805.  
  806.         tmp = bm[i] - blue ;
  807.         if( tmp < (MAX_COLOR/C_LEN/2) )
  808.           tmp = MAX_COLOR/C_LEN-1 - tmp ;
  809.         dist += tmp*tmp ;
  810.  
  811.         if( dist < mindist )
  812.           mindist = dist ;
  813.       }
  814.  
  815.  
  816.     /* step 3: find all points within that distance to box */
  817.  
  818.     for( i=0; i<num_colors; ++i )
  819.     {
  820.       register int    dist, tmp ;
  821.  
  822.       if( rm[i] >> (COLOR_DEPTH-C_DEPTH) != ir  ||
  823.           gm[i] >> (COLOR_DEPTH-C_DEPTH) != ig  ||
  824.           bm[i] >> (COLOR_DEPTH-C_DEPTH) != ib)
  825.       {
  826.         dist = 0 ;
  827.  
  828.         if( (tmp = red - rm[i]) > 0  ||
  829.         (tmp = rm[i] - (red + MAX_COLOR/C_LEN-1)) > 0 )
  830.           dist += tmp*tmp ;
  831.  
  832.         if( (tmp = green - gm[i]) > 0  ||
  833.         (tmp = gm[i] - (green + MAX_COLOR/C_LEN-1)) > 0 )
  834.           dist += tmp*tmp ;
  835.  
  836.         if( (tmp = blue - bm[i]) > 0  ||
  837.         (tmp = bm[i] - (blue + MAX_COLOR/C_LEN-1)) > 0 )
  838.           dist += tmp*tmp ;
  839.  
  840.         if( dist < mindist )
  841.         {
  842.           ptr->entries[ptr->num_ents][0] = i ;
  843.           ptr->entries[ptr->num_ents][1] = dist ;
  844.           ++ptr->num_ents ;
  845.         }
  846.       }
  847.     }
  848.  
  849.     /* sort color cells by distance, use cheap exchange sort */
  850.     {
  851.       register int    n, tmp, i ;
  852.       int        next_n ;
  853.  
  854.       n = ptr->num_ents - 1 ;
  855.       while( n > 0 )
  856.       {
  857.         next_n = 0 ;
  858.         for( i=0; i<n; ++i)
  859.         {
  860.           if( ptr->entries[i][1] > ptr->entries[i+1][1] )
  861.           {
  862.         tmp = ptr->entries[i][0] ;
  863.         ptr->entries[i][0] = ptr->entries[i+1][0] ;
  864.         ptr->entries[i+1][0] = tmp ;
  865.         tmp = ptr->entries[i][1] ;
  866.         ptr->entries[i][1] = ptr->entries[i+1][1] ;
  867.         ptr->entries[i+1][1] = tmp ;
  868.         next_n = i ;
  869.           }
  870.         }
  871.         n = next_n ;
  872.       }
  873.     }
  874.     return( ptr ) ;
  875. }
  876.  
  877.  
  878.  
  879.  
  880. static
  881. map_colortable()
  882. {
  883.     int        ir,ig,ib ;
  884. register int        *histp = &histogram[0][0][0] ;
  885. register C_cell        *cell ;
  886.  
  887.     for(ir = 0; ir < B_LEN; ++ir)
  888.       for(ig = 0; ig < B_LEN; ++ig)
  889.         for(ib = 0; ib < B_LEN; ++ib)
  890.         {
  891.           if( *histp == 0 )
  892.         *histp = -1 ;
  893.           else
  894.           {
  895.         int        i ;
  896.         register int    j, tmp, d2, dist ;
  897.  
  898.         cell = *( ColorCells + ( ((ir>>(B_DEPTH-C_DEPTH)) << C_DEPTH*2)
  899.                   + ((ig>>(B_DEPTH-C_DEPTH)) << C_DEPTH )
  900.               + (ib>>(B_DEPTH-C_DEPTH)) ) ) ;
  901.         
  902.         if(cell == NULL )
  903.           cell = create_colorcell( ir<<(COLOR_DEPTH-B_DEPTH),
  904.                        ig<<(COLOR_DEPTH-B_DEPTH),
  905.                        ib<<(COLOR_DEPTH-B_DEPTH) ) ;
  906.  
  907.         dist = 9999999 ;
  908.         for(i = 0;
  909.             i < cell->num_ents  &&  dist > cell->entries[i][1] ;
  910.             ++i)
  911.         {
  912.           j = cell->entries[i][0] ;
  913.           d2 = rm[j] - (ir << (COLOR_DEPTH-B_DEPTH)) ;
  914.           d2 *= d2 ;
  915.           tmp = gm[j] - (ig << (COLOR_DEPTH-B_DEPTH)) ;
  916.           d2 += tmp*tmp ;
  917.           tmp = bm[j] - (ib << (COLOR_DEPTH-B_DEPTH)) ;
  918.           d2 += tmp*tmp ;
  919.           if( d2 < dist )
  920.           {
  921.             dist = d2 ;
  922.             *histp = j ;
  923.           }
  924.         }
  925. #ifdef    DEBUG
  926.         if( dist == 9999999 )
  927.           fprintf(stderr,"??? map_colortable: dist=9999999\n") ;
  928. #endif    DEBUG
  929.           }
  930.           ++histp ;
  931.         }
  932. }
  933.  
  934.  
  935.  
  936.  
  937.  
  938.  
  939.  
  940.  
  941.  
  942.  
  943. /*
  944.  * straight quantization.  Each pixel is mapped to the colors
  945.  * closest to it.  Color values are rounded to the nearest color
  946.  * table entry.
  947.  */
  948.  
  949. static
  950. quant()
  951. {
  952.     unsigned char    *outline ;
  953.     unsigned char    *inline ;
  954. register unsigned char    *outptr ;
  955. register unsigned char    *inptr ;
  956. register int        i, j ;
  957.  
  958.  
  959.  
  960.  
  961.     inline = (unsigned char *) alloca( rh.ras_width * bytes_per_pixel ) ;
  962.  
  963.     outline = (unsigned char *) alloca( rh.ras_width ) ;
  964.  
  965.     for(i=0; i < rh.ras_height; ++i)
  966.     {
  967.       pr_get_bytes(stdin, &rh, rh.ras_width * bytes_per_pixel, inline) ;
  968.       inptr = inline ;
  969.       outptr = outline;
  970.       for(j=0; j < rh.ras_width; ++j)
  971.       {
  972.         register int    tmp,red,green,blue ;
  973.  
  974.         if(rh.ras_depth == 32)
  975.           ++inptr ;
  976.         blue = *inptr++ ;
  977.         green = *inptr++ ;
  978.         red = *inptr++ ;
  979.  
  980.         if (rh.ras_maplength > 0)
  981.         {
  982.           red = *(colormap.map[0]+red) ;
  983.           green = *(colormap.map[1]+green) ;
  984.           blue = *(colormap.map[2]+blue) ;
  985.         }
  986.  
  987.         red >>= (COLOR_DEPTH-B_DEPTH) ;
  988.         green >>= (COLOR_DEPTH-B_DEPTH) ;
  989.         blue >>= (COLOR_DEPTH-B_DEPTH) ;
  990.  
  991.         *outptr++ = histogram[red][green][blue] ;
  992.       }
  993.       fwrite(outline, 1, rh.ras_width, stdout) ;
  994.     }
  995. }
  996.  
  997.  
  998.  
  999.  
  1000.  
  1001. static
  1002. quant_fsdither()
  1003. {
  1004.     unsigned char    *outline, *inline ;
  1005.     short        *thisline, *nextline, *tmpptr ;
  1006.     unsigned char    *inptr ;
  1007. register unsigned char    *outptr ;
  1008. register short        *thisptr, *nextptr ;
  1009. register int        i, j ;
  1010.     int        imax, jmax ;
  1011.     int        lastline, lastpixel ;
  1012.  
  1013.  
  1014.     imax = rh.ras_height - 1 ;
  1015.     jmax = rh.ras_width - 1 ;
  1016.  
  1017.     inline = (unsigned char *) alloca( rh.ras_width * bytes_per_pixel ) ;
  1018.     thisline = (short *) alloca( rh.ras_width * 3 * sizeof(short)) ;
  1019.     nextline = (short *) alloca( rh.ras_width * 3 * sizeof(short)) ;
  1020.  
  1021.     outline = (unsigned char *) alloca( rh.ras_width ) ;
  1022.  
  1023.     /*
  1024.      * get first line
  1025.      */
  1026.     pr_get_bytes(stdin, &rh, rh.ras_width * bytes_per_pixel, inline) ;
  1027.     {
  1028.       register int tmp ;
  1029.       inptr = inline ;
  1030.       nextptr = nextline ;
  1031.       for(j=0; j < rh.ras_width; ++j)
  1032.       {
  1033.         if(rh.ras_depth == 32)
  1034.           ++inptr ;
  1035.         if( rh.ras_maplength > 0)
  1036.         {
  1037.           *nextptr++ = *(colormap.map[0]+*inptr++) ;
  1038.           *nextptr++ = *(colormap.map[1]+*inptr++) ;
  1039.           *nextptr++ = *(colormap.map[2]+*inptr++) ;
  1040.         }
  1041.         else
  1042.         {
  1043.           *nextptr++ = *inptr++ ;
  1044.           *nextptr++ = *inptr++ ;
  1045.           *nextptr++ = *inptr++ ;
  1046.         }
  1047.       }
  1048.     }
  1049.  
  1050.     for(i=0; i < rh.ras_height; ++i)
  1051.     {
  1052.       tmpptr = thisline ;
  1053.       thisline = nextline ;
  1054.       nextline = tmpptr ;
  1055.       lastline = (i == imax) ;
  1056.       pr_get_bytes(stdin, &rh, rh.ras_width * bytes_per_pixel, inline) ;
  1057.       {
  1058.         register int tmp ;
  1059.         inptr = inline ;
  1060.         nextptr = nextline ;
  1061.         for(j=0; j < rh.ras_width; ++j)
  1062.         {
  1063.           if(rh.ras_depth == 32)
  1064.         ++inptr ;
  1065.           if( rh.ras_maplength > 0)
  1066.           {
  1067.         *nextptr++ = *(colormap.map[0]+*inptr++) ;
  1068.         *nextptr++ = *(colormap.map[1]+*inptr++) ;
  1069.         *nextptr++ = *(colormap.map[2]+*inptr++) ;
  1070.           }
  1071.           else
  1072.           {
  1073.         *nextptr++ = *inptr++ ;
  1074.         *nextptr++ = *inptr++ ;
  1075.         *nextptr++ = *inptr++ ;
  1076.           }
  1077.         }
  1078.       }
  1079.  
  1080.       thisptr = thisline ;
  1081.       nextptr = nextline ;
  1082.       outptr = outline;
  1083.  
  1084.       for(j=0; j < rh.ras_width ; ++j)
  1085.       {
  1086.         int            red,green,blue ;
  1087.         register int    oval, r2,g2,b2 ;
  1088.  
  1089.         lastpixel = (j == jmax) ;
  1090.  
  1091.         b2 = *thisptr++ ;
  1092.         g2 = *thisptr++ ;
  1093.         r2 = *thisptr++ ;
  1094.  
  1095.  
  1096.         if( r2 < 0 ) r2 = 0 ;
  1097.         else if( r2 >= MAX_COLOR ) r2 = MAX_COLOR-1 ;
  1098.         if( g2 < 0 ) g2 = 0 ;
  1099.         else if( g2 >= MAX_COLOR ) g2 = MAX_COLOR-1 ;
  1100.         if( b2 < 0 ) b2 = 0 ;
  1101.         else if( b2 >= MAX_COLOR ) b2 = MAX_COLOR-1 ;
  1102.  
  1103.         red = r2 ;
  1104.         green = g2 ;
  1105.         blue = b2 ;
  1106.  
  1107.         r2 >>= (COLOR_DEPTH-B_DEPTH) ;
  1108.         g2 >>= (COLOR_DEPTH-B_DEPTH) ;
  1109.         b2 >>= (COLOR_DEPTH-B_DEPTH) ;
  1110.  
  1111.         if( (oval = histogram[r2][g2][b2]) == -1)
  1112.         {
  1113.           int        ci ;
  1114.           register int    cj, tmp, d2, dist ;
  1115.           register C_cell    *cell ;
  1116.  
  1117.           cell = *( ColorCells + ( ((r2>>(B_DEPTH-C_DEPTH)) << C_DEPTH*2)
  1118.             + ((g2>>(B_DEPTH-C_DEPTH)) << C_DEPTH )
  1119.             + (b2>>(B_DEPTH-C_DEPTH)) ) ) ;
  1120.           
  1121.           if(cell == NULL )
  1122.         cell = create_colorcell( red, green, blue ) ;
  1123.  
  1124.  
  1125.           dist = 9999999 ;
  1126.           for(ci = 0;
  1127.           ci < cell->num_ents  &&  dist > cell->entries[ci][1] ;
  1128.           ++ci)
  1129.           {
  1130.         cj = cell->entries[ci][0] ;
  1131.         d2 = (rm[cj] >> (COLOR_DEPTH-B_DEPTH)) - r2 ;
  1132.         d2 *= d2 ;
  1133.         tmp = (gm[cj] >> (COLOR_DEPTH-B_DEPTH)) - g2 ;
  1134.         d2 += tmp*tmp ;
  1135.         tmp = (bm[cj] >> (COLOR_DEPTH-B_DEPTH)) - b2 ;
  1136.         d2 += tmp*tmp ;
  1137.         if( d2 < dist )
  1138.         {
  1139.           dist = d2 ;
  1140.           oval = cj ;
  1141.         }
  1142.           }
  1143.           histogram[r2][g2][b2] = oval ;
  1144. #ifdef    DEBUG
  1145.           if( dist == 9999999 )
  1146.         fprintf(stderr,"??? quant: dist=9999999\n") ;
  1147. #endif    DEBUG
  1148.         }
  1149.  
  1150.         *outptr++ = oval ;
  1151.  
  1152.         red -= rm[oval] ;
  1153.         green -= gm[oval] ;
  1154.         blue -= bm[oval] ;
  1155.  
  1156.  
  1157.         if( !lastpixel )
  1158.         {
  1159.           thisptr[0] += blue * 7 / 16 ;
  1160.           thisptr[1] += green * 7 / 16 ;
  1161.           thisptr[2] += red * 7 / 16 ;
  1162.         }
  1163.         if( !lastline )
  1164.         {
  1165.           if( j != 0 )
  1166.           {
  1167.         nextptr[-3] += blue * 3 / 16 ;
  1168.         nextptr[-2] += green * 3 / 16 ;
  1169.         nextptr[-1] += red * 3 / 16 ;
  1170.           }
  1171.           nextptr[0] += blue * 5 / 16 ;
  1172.           nextptr[1] += green * 5 / 16 ;
  1173.           nextptr[2] += red * 5 / 16 ;
  1174.           if( !lastpixel )
  1175.           {
  1176.         nextptr[3] += blue / 16 ;
  1177.         nextptr[4] += green / 16 ;
  1178.         nextptr[5] += red / 16 ;
  1179.           }
  1180.           nextptr += 3 ;
  1181.         }
  1182.       }
  1183.       fwrite(outline, 1, rh.ras_width, stdout) ;
  1184.     }
  1185. }
  1186.  
  1187.  
  1188.  
  1189.  
  1190. #ifdef    DEBUG
  1191. showboxes()
  1192. {
  1193.     Colorbox    *ptr ;
  1194.  
  1195.     printf("boxes:\n") ;
  1196.     ptr = usedboxes ;
  1197.     while(ptr != NULL)
  1198.     {
  1199.       printf("  %d<r<%d, %d<g<%d, %d<b<%d\n", ptr->rmin,ptr->rmax,
  1200.           ptr->gmin,ptr->gmax, ptr->bmin,ptr->bmax) ;
  1201.       ptr = ptr->next ;
  1202.     }
  1203. }
  1204.  
  1205.  
  1206.  
  1207. showcell(ptr)
  1208. register C_cell    *ptr ;
  1209. {
  1210.     int    i ;
  1211.     int    j ;
  1212.  
  1213.     fprintf(stderr, "n=%d\n", ptr->num_ents) ;
  1214.     for(i=0; i<ptr->num_ents; ++i)
  1215.     {
  1216.       j = ptr->entries[i][0] ;
  1217.       fprintf(stderr, "(%d,%d): (%d,%d,%d)\n",
  1218.           j, ptr->entries[i][1], rm[j],gm[j],bm[j]) ;
  1219.     }
  1220. }
  1221.  
  1222.  
  1223.  
  1224. int
  1225. testbox(box)
  1226. register Colorbox    *box ;
  1227. {
  1228. register int        ir,ig,ib, total ;
  1229.  
  1230.     total = 0 ;
  1231.     for(ir=box->rmin; ir<=box->rmax; ++ir)
  1232.       for(ig=box->gmin; ig<=box->gmax; ++ig)
  1233.         for(ib=box->bmin; ib<=box->bmax; ++ib)
  1234.           total += histogram[ir][ig][ib] ;
  1235.  
  1236.     return total ;
  1237. }
  1238. #endif    DEBUG
  1239.  
  1240.  
  1241.  
  1242.  
  1243.  
  1244.  
  1245. /* pr_stream routines:  allow images to be read in stream
  1246.  * fasion rather than loading them into memory.
  1247.  *
  1248.  * These routines are used in conjunction with the pr_io routines.  It
  1249.  * is the application's responsibility to read the header and
  1250.  * colormaps
  1251.  *
  1252.  *
  1253.  * routines:
  1254.  *
  1255.  * pr_read_init(header)
  1256.  *    struct rasterfile *header ;
  1257.  *
  1258.  *    sets up to read an image
  1259.  *
  1260.  *
  1261.  * pr_get_bytes(file, header, len, buffer)
  1262.  *    FILE        *file ;
  1263.  *    struct rasterfile *header ;
  1264.  *    int        len ;
  1265.  *    char        *buffer ;
  1266.  *
  1267.  *    read data from input
  1268.  *
  1269.  *
  1270.  * int
  1271.  * pr_get_byte(file, header)
  1272.  *    FILE        *file ;
  1273.  *    struct rasterfile *header ;
  1274.  *
  1275.  *    read and return one byte from input
  1276.  *
  1277.  */
  1278.  
  1279. #include <stdio.h>
  1280. #include <sys/types.h>
  1281. #include <pixrect/pixrect.h>
  1282. #include <pixrect/memvar.h>
  1283. #include <pixrect/pr_io.h>
  1284. #include <rasterfile.h>
  1285.  
  1286. #define    ESCAPE    128
  1287.  
  1288.  
  1289.  
  1290. static    unsigned char    icvalue ;
  1291. static    int        icount ;
  1292.  
  1293.  
  1294.  
  1295. int
  1296. pr_read_init(header)
  1297.     struct rasterfile    *header ;
  1298. {
  1299.     icount = 0 ;
  1300. }
  1301.  
  1302.  
  1303.  
  1304. int
  1305. pr_get_bytes(file, header, len, buffer)
  1306. register FILE            *file ;
  1307. register struct rasterfile    *header ;
  1308. register int            len ;
  1309. register unsigned char        *buffer ;
  1310. {
  1311.     if(header->ras_type != RT_BYTE_ENCODED)
  1312.       return( fread(buffer, 1, len, file) ) ;
  1313.  
  1314.     else while(len-- > 0)
  1315.     {
  1316.       if(icount > 0)
  1317.       {
  1318.         --icount ;
  1319.         *buffer++ = icvalue ;
  1320.       }
  1321.       else
  1322.       {
  1323.         if( (icvalue = getc(file)) != ESCAPE )
  1324.           *buffer++ = icvalue ;
  1325.         else
  1326.         {
  1327.           if( (icount = getc(file)) == 0 )
  1328.         *buffer++ = ESCAPE ;
  1329.           else
  1330.           {
  1331.         *buffer++ = icvalue = getc(file) ;
  1332.           }
  1333.         }
  1334.       }
  1335.     }
  1336.     return 0 ;
  1337. }
  1338.  
  1339.  
  1340.  
  1341. unsigned char
  1342. pr_get_byte(file, header)
  1343.     FILE            *file ;
  1344.     struct rasterfile    *header ;
  1345. {
  1346.     unsigned char        rval ;
  1347.  
  1348.     pr_get_bytes(file, header, 1, &rval) ;
  1349.     return( rval ) ;
  1350. }
  1351. -- 
  1352.         -ed falk, sun microsystems
  1353.          sun!falk, falk@sun.com
  1354. terrorist, cryptography, DES, drugs, cipher, secret, decode, NSA, CIA, NRO.
  1355.