home *** CD-ROM | disk | FTP | other *** search
/ Magazyn Amiga Shareware Floppies / ma01.dms / ma01.adf / wasp / src / wriffdistr.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-01-01  |  28.8 KB  |  1,415 lines
  1. /* wasp - Copyright 1991 by Steven Reiz
  2.  * see COPYING and wasp.c for further info
  3.  * wriffdistr.c, 4/12/90 - 30/1/91, 2/6/91, 23/6/91,
  4.  * 3/7/91 - 10/7/91, 8/12/91, 1/1/92
  5.  * features: SORT_SLICED, SHAM_BLACK
  6.  */
  7.  
  8. struct hs_t {
  9.     char dist;
  10.     char pad;
  11.     short color;
  12. };
  13.  
  14. static char *sourcefile=__FILE__;
  15.  
  16. #include "wasp.h"
  17. #include "wriff.h"
  18. #define SORT_SLICED
  19. /* #define SHAM_BLACK */
  20.  
  21. static char regind;
  22. static long squares[]={
  23.     0, 1, 4, 9, 16, 25, 36, 49, 64, 81, 100, 121, 144, 169, 196, 225
  24. };
  25.  
  26.  
  27. void
  28. compute_distr(int firstrow, int lastrow)
  29. {
  30.     assert(counts!=NULL1);
  31.     if (icolors<=nregs)
  32.         no_distr();
  33.     else {
  34.         switch (distrmeth) {
  35.         case DISTRMETH_MOSTUSED:
  36.             mostused_distr();
  37.         break;
  38.         case DISTRMETH_WORSTFIRST:
  39.             worstfirst_distr();
  40.         break;
  41.         case DISTRMETH_EHB:
  42.             ehb_distr();
  43.         break;
  44.         case DISTRMETH_MUE:
  45.             mue_distr();
  46.         break;
  47.         case DISTRMETH_HAMSHARP:
  48.             hamsharp_distr();
  49.         break;
  50.         case DISTRMETH_CONTRACTION:
  51.             contraction_distr();
  52.         break;
  53.         default:
  54.         error0(E0_FATAL, E1_IFF, E2_OPTION, E3_DISTRMETH);
  55.         break;
  56.         }
  57.     }
  58.     if (xmode!=EHB) {
  59. #ifdef SORT_SLICED
  60.         sort_cm();
  61. #else
  62.         if (!slicedo)
  63.             sort_cm();
  64. #endif
  65.     }
  66.     if (slicedo) {
  67.     if (sliced==SLICED_PCHG)
  68.         restrict_changes=PCHG_NCHANGES;
  69.     if (restrict_changes!= -1)
  70.         restrict_cm();
  71.     }
  72.     count_pixels(firstrow, lastrow, 1);
  73.     fill_newcol();
  74. }
  75.  
  76.  
  77. void
  78. no_distr(void)
  79. {
  80.     short i;
  81.     u_long th, *p;
  82.     u_short *q;
  83.  
  84.     p=counts+NRGB;
  85.     th=threshold;
  86.     q=curcm;
  87.     i=NRGB-1;
  88.     if (th<=1) {
  89.         do {
  90.             if (*--p)
  91.                 *q++ =i;
  92.         } while (--i>=0);
  93.     } else {
  94.         do {
  95.             if (*--p >=th)
  96.                 *q++ =i;
  97.         } while (--i>=0);
  98.     }
  99.     assert(q==curcm+icolors);
  100.     while (q<curcm+nregs)
  101.         *q++ =0;
  102. }
  103.  
  104.  
  105. void
  106. mostused_distr(void)
  107. {
  108.     static u_long regcounts[MAXNREGS];
  109.     u_long *countp, *countp2, th;
  110.     static u_long *savecountp;
  111.     u_short *cmp;
  112.     short i;
  113.     
  114.     countp=regcounts;
  115.     i=nregs-1;
  116.     do {
  117.         *countp++ =0;
  118.     } while (--i>=0);
  119.     th=0;
  120.     countp=counts+NRGB;
  121.     i=NRGB-1;
  122.     do {
  123.         if (*--countp>th) {
  124.             th= *countp;
  125.             savecountp=countp;
  126.             countp=regcounts;
  127.             do {
  128.                 /* tight loop! */
  129.             } while (*countp++ >=th);
  130.             countp2=regcounts+nregs;
  131.             cmp=curcm+nregs;
  132.             if (countp2>countp) {
  133.                 do {
  134.                     *--countp2= *(countp2-2);
  135.                     *--cmp    = *(cmp-2);
  136.                 } while (countp2>countp);
  137.             }
  138.             *--countp2=th;
  139.             *--cmp=i;
  140.             th=regcounts[nregs-1];
  141.             countp=savecountp;
  142.         }
  143.     } while (--i>=0);
  144. #ifdef SHAM_BLACK
  145.     if (slicedo) {
  146.         for (i=0; i<nregs; ++i)
  147.             if (curcm[i]==0x000)
  148.                 break;
  149.         if (i>=nregs) { /* black isn't amongst the colors, add it */
  150.             curcm[nregs-1]=0x000;
  151.         }
  152.     }
  153. #endif
  154. }
  155.  
  156.  
  157. void
  158. worstfirst_distr(void)
  159. {
  160.     static int sicolors= -1;
  161.     static int startreg;
  162.     int i;
  163.     short rgb1, rgb2;
  164.  
  165.     if (sicolors<icolors) {
  166.     if (sicolors== -1) {
  167.             if (icolors<128)
  168.                 sicolors=128;
  169.             else
  170.                 sicolors=icolors;
  171.         } else {
  172.             free(wfmindist);
  173.             free(wfminind);
  174.             free(wfrgb);
  175.             if (icolors<256)
  176.                 sicolors=256;
  177.             else
  178.                 sicolors=icolors;
  179.         }
  180.         wfmindist=Malloc(sicolors);
  181.         wfminind=Malloc(sicolors);
  182.         wfrgb=Malloc(sicolors*3);
  183.     }
  184.     fill_wf_rgb();
  185.     startreg=2;
  186. #ifdef SHAM_BLACK
  187.     if (slicedo) {
  188.     startreg=1;
  189.     rgb1=0x000; regind=0;
  190.     fill_wf_dist_1(rgb1);
  191.     } else {
  192. #endif
  193.     find_2_furthest(&rgb1, &rgb2);
  194.     curcm[0]=rgb1; regind=0;
  195.     fill_wf_dist_1(rgb1);
  196.     curcm[1]=rgb2; regind=1;
  197.     fill_wf_dist_2(rgb2);
  198. #ifdef SHAM_BLACK
  199.     }
  200. #endif
  201.     for (i=startreg; i<nregs; ++i) {
  202.         find_1_furthest(&rgb1);
  203.         curcm[i]=rgb1; regind=i;
  204.         fill_wf_dist_2(rgb1);
  205.     }
  206.     fill_wf2rgbweight();
  207.     wf2_redo_curcm();
  208. }
  209.  
  210.  
  211. void
  212. fill_wf_rgb(void)
  213. {
  214.     long i;
  215.     u_long th, *p;
  216.     char *q;
  217.  
  218.     p=counts+NRGB;
  219.     th=threshold;
  220.     q=wfrgb;
  221.     i=NRGB-1;
  222.     if (th<=1) {
  223.         do {
  224.             if (*--p) {
  225.                 *q++ =i>>8;
  226.                 *q++ =(i>>4)&0x0f;
  227.                 *q++ =i&0x0f;
  228.             }
  229.         } while (--i>=0);
  230.     } else {
  231.         do {
  232.             if (*--p >=th) {
  233.                 *q++ =i>>8;
  234.                 *q++ =(i>>4)&0x0f;
  235.                 *q++ =i&0x0f;
  236.             }
  237.         } while (--i>=0);
  238.     }
  239.     assert(q==wfrgb+3*icolors);
  240. }
  241.  
  242.  
  243. void
  244. find_2_furthest(short *rgb1p, short *rgb2p)
  245. {
  246.     REG char *p, *q;
  247.     REG char r, g, b, t, dist, bigdist;
  248.     NON_REG char *bigp, *bigq;
  249.  
  250.     bigdist=0;
  251.     p=wfrgb+3*icolors;
  252.     while (p>wfrgb) {
  253.         b= *--p;
  254.         g= *--p;
  255.         r= *--p;
  256.         q=wfrgb;
  257.         do {
  258.             if ((dist=r- *q++)<0) dist= -dist;
  259.             if ((t=g- *q++)<0)    dist-=t; else dist+=t;
  260.             if ((t=b- *q++)<0)    dist-=t; else dist+=t;
  261.             if (dist>bigdist) {
  262.                 bigdist=dist;
  263.                 bigp=p;
  264.                 bigq=q;
  265.             }
  266.         } while (q<p);
  267.     }
  268.     *rgb1p=(*bigp<<8)    |(*(bigp+1)<<4)|(*(bigp+2));
  269.     *rgb2p=(*(bigq-3)<<8)|(*(bigq-2)<<4)|(*(bigq-1));
  270. }
  271.  
  272.  
  273. void
  274. fill_wf_dist_1(short color)
  275. {
  276.     short i;
  277.     char r, g, b, t, dist;
  278.     char *mdp, *mip, *rgbp;
  279.  
  280.     r=color>>8;
  281.     g=(color>>4)&0x0f;
  282.     b=color&0x0f;
  283.     mdp=wfmindist+icolors;
  284.     mip=wfminind+icolors;
  285.     rgbp=wfrgb+3*icolors;
  286.     i=icolors-1;
  287.     do {
  288.         if ((dist=b- *--rgbp)<0) dist= -dist;
  289.         if ((t=g- *--rgbp)<0)    dist-=t; else dist+=t;
  290.         if ((t=r- *--rgbp)<0)    dist-=t; else dist+=t;
  291.         *--mdp=dist;
  292.         *--mip=regind;
  293.     } while (--i>=0);
  294. }
  295.  
  296.  
  297. void
  298. fill_wf_dist_2(short color)
  299. {
  300.     short i;
  301.     char r, g, b, t, dist;
  302.     char *mdp, *mip, *rgbp;
  303.  
  304.     r=color>>8;
  305.     g=(color>>4)&0x0f;
  306.     b=color&0x0f;
  307.     mdp=wfmindist+icolors;
  308.     mip=wfminind+icolors;
  309.     rgbp=wfrgb+3*icolors;
  310.     i=icolors-1;
  311.     do {
  312.         if ((dist=b- *--rgbp)<0) dist= -dist;
  313.         if ((t=g- *--rgbp)<0)    dist-=t; else dist+=t;
  314.         if ((t=r- *--rgbp)<0)    dist-=t; else dist+=t;
  315.         if (dist< *--mdp) {
  316.             *mdp=dist;
  317.             *--mip=regind;
  318.         } else
  319.             --mip;
  320.     } while (--i>=0);
  321. }
  322.  
  323.  
  324. void
  325. find_1_furthest(short *rgbp)
  326. {
  327.     char *p, max;
  328.     short i, maxi;
  329.  
  330.     p=wfmindist+icolors;
  331.     max=0;
  332.     i=icolors-1;
  333.     do {
  334.         if (*--p>max) {
  335.             max= *p;
  336.             maxi=i;
  337.         }
  338.     } while (--i>=0);
  339.     p=wfrgb+3*maxi;
  340.     *rgbp=(*p<<8)|(*(p+1)<<4)|(*(p+2));
  341. }
  342.  
  343.  
  344. void
  345. fill_wf2rgbweight(void)
  346. {
  347.     char *rgbp, *mip;
  348.     u_long *weightp, weight;
  349.     short i;
  350.     char r, g, b;
  351.  
  352.     if (wf2rgbweight==NULL1)
  353.         wf2rgbweight=Malloc(nregs*4*sizeof(u_long));
  354.     weightp=wf2rgbweight;
  355.     i=nregs-1;
  356.     do {
  357.         *weightp++ =0;
  358.         *weightp++ =0;
  359.         *weightp++ =0;
  360.         *weightp++ =0;
  361.     } while (--i>=0);
  362.     rgbp=wfrgb;
  363.     mip=wfminind;
  364.     i=icolors-1;
  365.     do {
  366.         r= *rgbp++;
  367.         g= *rgbp++;
  368.         b= *rgbp++;
  369.         weight=counts[(r<<8)|(g<<4)|b];
  370.         weightp=wf2rgbweight+4* (*mip++);
  371.         *weightp++ += r*weight;
  372.         *weightp++ += g*weight;
  373.         *weightp++ += b*weight;
  374.         *weightp   += weight;
  375.     } while (--i>=0);
  376. }
  377.  
  378.  
  379. void
  380. wf2_redo_curcm(void)
  381. {
  382.     u_long *weightp, weight, r, g, b;
  383.     short i;
  384.     u_short *cmp;
  385.     
  386.     weightp=wf2rgbweight;
  387.     cmp=curcm;
  388.     i=nregs-1;
  389.     do {
  390.         r= *weightp++;
  391.         g= *weightp++;
  392.         b= *weightp++;
  393.         weight= *weightp++;
  394.         if (!weight)
  395.             weight=1;
  396.         r=(r+weight/2)/weight;
  397.         if (r>15) r=15;
  398.         g=(g+weight/2)/weight;
  399.         if (g>15) g=15;
  400.         b=(b+weight/2)/weight;
  401.         if (b>15) b=15;
  402.         *cmp++ =(r<<8)|(g<<4)|b;
  403.     } while (--i>=0);
  404. #ifdef SHAM_BLACK
  405.     if (slicedo)
  406.     black_darkest();
  407. #endif
  408. }
  409.  
  410.  
  411. void
  412. black_darkest(void)
  413. {
  414.     short i, besti, val, bestval;
  415.  
  416.     bestval=100;
  417.     for (i=0; i<nregs; ++i) {
  418.         val=(curcm[i]&0xf)+((curcm[i]>>4)&0xf)+(curcm[i]>>8);
  419.         if (val<bestval) {
  420.             bestval=val;
  421.             besti=i;
  422.         }
  423.     }
  424.     curcm[besti]=0x000;
  425. }
  426.  
  427.  
  428. void
  429. ehb_distr(void)
  430. {
  431.     static int sicolors= -1;
  432.     int i, nr;
  433.     short rgb1, rgb2;
  434.  
  435.     nr=nregs/2;
  436.     if (sicolors<icolors) {
  437.         if (sicolors== -1) {
  438.             if (icolors<128)
  439.                 sicolors=128;
  440.             else
  441.                 sicolors=icolors;
  442.         } else {
  443.             free(wfmindist);
  444.             free(wfminind);
  445.             free(wfrgb);
  446.             if (icolors<256)
  447.                 sicolors=256;
  448.             else
  449.                 sicolors=icolors;
  450.         }
  451.         wfmindist=Malloc(sicolors);
  452.         wfminind=Malloc(sicolors);
  453.         wfrgb=Malloc(sicolors*3);
  454.     }
  455.     fill_wf_rgb();
  456.     find_2_furthest(&rgb1, &rgb2);
  457.     regind=0;
  458.     fill_wf_dist_1(rgb1);
  459.     ehb_pair(rgb1, 0);
  460.     for (i=1; i<nr; ++i) {
  461.         find_1_furthest(&rgb1);
  462.         regind=i;
  463.         fill_wf_dist_2(rgb1);
  464.         ehb_pair(rgb1, i);
  465.     }
  466.     fill_wf2rgbweight();
  467.     ehb_redo_curcm();
  468. }
  469.  
  470.  
  471. void
  472. ehb_pair(short c1, int i)
  473. {
  474.     short c2, c3, c4;
  475.     short r, g, b;
  476.     int nr;
  477.  
  478.     nr=nregs/2;
  479.     r=c1>>8;
  480.     g=(c1>>4)&0x0f;
  481.     b=c1&0x0f;
  482.     c2=((r>>1)<<8)|((g>>1)<<4)|(b>>1);
  483.     if (r<8 && g<8 && b<8) {
  484.         c3=(r<<9)|(g<<5)|(b<<1);
  485.         ehb_worst(c2, c3, &c4);
  486.         if (c4==c2) {
  487.             curcm[i]=c1;
  488.             curcm[nr+i]=c2;
  489.             regind=nr+i;
  490.         } else {
  491.             curcm[i]=c4;
  492.             curcm[nr+i]=c1;
  493.             c2=c4;
  494.             ehb_minind_adjust((int)i, (int)(nr+i));
  495.             regind=i;
  496.         }
  497.     } else {
  498.         curcm[i]=c1;
  499.         curcm[nr+i]=c2;
  500.         regind=nr+i;
  501.     }
  502.     fill_wf_dist_2(c2);
  503. }
  504.  
  505.  
  506. void
  507. ehb_worst(short col1, short col8, short *worstp)
  508. {
  509.     short curcol;
  510.     int curdist, bestdist;
  511.     short i;
  512.  
  513.     bestdist=ehb_farcol(col1);
  514.     *worstp=col1;
  515.     i=7;
  516.     do {
  517.         curcol=col8;
  518.         if (i&1)
  519.             curcol|=0x001;
  520.         if (i&2)
  521.             curcol|=0x010;
  522.         if (i&4)
  523.             curcol|=0x100;
  524.         curdist=ehb_farcol(curcol);
  525.         if (curdist>bestdist) {
  526.             bestdist=curdist;
  527.             *worstp=curcol;
  528.         }
  529.     } while (--i>=0);
  530. }
  531.  
  532.  
  533. int
  534. ehb_farcol(short color)
  535. {
  536.     char r, g, b;
  537.     char *p, *q;
  538.     char t, dist, bigdist;
  539.     short i;
  540.  
  541.     r=color>>8;
  542.     g=(color>>4)&0x0f;
  543.     b=color&0x0f;
  544.     bigdist=0;
  545.     i=icolors-1;
  546.     p=wfmindist;
  547.     q=wfrgb;
  548.     do {
  549.         if ((dist=r- *q++)>0) dist= -dist;
  550.         if ((t=g- *q++)>0)    dist-=t; else dist+=t;
  551.         if ((t=b- *q++)>0)    dist-=t; else dist+=t;
  552.         dist+= *p++;
  553.         if (dist>bigdist)
  554.             bigdist=dist;
  555.     } while (--i>=0);
  556.     return (int)bigdist;
  557. }
  558.  
  559.  
  560. void
  561. ehb_minind_adjust(int from, int to)
  562. {
  563.     char *p;
  564.     char fromc, toc;
  565.     short i;
  566.  
  567.     p=wfminind+icolors;
  568.     i=icolors-1;
  569.     fromc=from;
  570.     toc=to;
  571.     do {
  572.         if (*--p ==fromc)
  573.             *p=toc;
  574.     } while (--i>=0);
  575. }
  576.  
  577.  
  578. void
  579. ehb_redo_curcm(void)
  580. {
  581.     short nr, i;
  582.     long r1, g1, b1, w1, r2, g2, b2, w2;
  583.     u_long *p1, *p2;
  584.  
  585.     nr=nregs/2;
  586.     p1=wf2rgbweight;
  587.     p2=wf2rgbweight+4*nr;
  588.     for (i=0; i<nr; ++i) {
  589.         r1= *p1++;     g1= *p1++;     b1= *p1++;     w1= *p1++;
  590.         r2=(*p2++)<<1; g2=(*p2++)<<1; b2=(*p2++)<<1; w2= *p2++;
  591.         w1=w1+w2;
  592.         assert(w1);
  593.         w2=w1/2;
  594.         r1=(r1+r2+w2)/w1;
  595.         if (r1>15) r1=15;
  596.         g1=(g1+g2+w2)/w1;
  597.         if (g1>15) g1=15;
  598.         b1=(b1+b2+w2)/w1;
  599.         if (b1>15) b1=15;
  600.         curcm[i]=(r1<<8)|(g1<<4)|b1;
  601.         r1>>=1; g1>>=1; b1>>=1;
  602.         curcm[nr+i]=(r1<<8)|(g1<<4)|b1;
  603.     }
  604. }
  605.  
  606.  
  607. void
  608. mue_distr(void)
  609. {
  610.     static u_long regcounts[MAXNREGS];
  611.     u_long *countp, *countp2, th;
  612.     static u_long *savecountp;
  613.     u_short *cmp;
  614.     short i, nr;
  615.     
  616.     nr=nregs/2;
  617.     countp=regcounts;
  618.     i=nr-1;
  619.     do {
  620.         *countp++ =0;
  621.     } while (--i>=0);
  622.     th=0;
  623.     countp=counts+NRGB;
  624.     i=NRGB-1;
  625.     do {
  626.         if (*--countp>th) {
  627.             th= *countp;
  628.             savecountp=countp;
  629.             countp=regcounts;
  630.             do {
  631.                 /* tight loop! */
  632.             } while (*countp++ >=th);
  633.             countp2=regcounts+nr;
  634.             cmp=curcm+nr;
  635.             if (countp2>countp) {
  636.                 do {
  637.                     *--countp2= *(countp2-2);
  638.                     *--cmp    = *(cmp-2);
  639.                 } while (countp2>countp);
  640.             }
  641.             *--countp2=th;
  642.             *--cmp=i;
  643.             th=regcounts[nr-1];
  644.             countp=savecountp;
  645.         }
  646.     } while (--i>=0);
  647.     mue_extend();
  648. }
  649.  
  650.  
  651. void
  652. mue_extend(void)
  653. {
  654.     short i;
  655.     short r, g, b, color;
  656.     u_short *p, *q;
  657.  
  658.     i=nregs/2;
  659.     p=curcm;
  660.     q=p+i;
  661.     --i;
  662.     do {
  663.         color= *p++;
  664.         r=color>>9;
  665.         g=(color>>5)&0x07;
  666.         b=(color>>1)&0x07;
  667.         *q++ = (r<<8)|(g<<4)|b;
  668.     } while (--i>=0);
  669. }
  670.  
  671.  
  672. void
  673. hamsharp_distr(void)
  674. {
  675.     static int sicolors= -1;
  676.     int i;
  677.  
  678.     if (sicolors<icolors) {
  679.         if (sicolors!= -1) {
  680.             for (i=0; i<sicolors; ++i)
  681.                 free(hshead[i]);
  682.             free(hshead);
  683.             free(hsmark);
  684.             free(wfminind);
  685.             free(wfrgb);
  686.         }
  687.         sicolors=icolors;
  688.         wfrgb=Malloc(sicolors*3);
  689.         wfminind=Malloc(sicolors);
  690.         hsmark=Malloc(sicolors);
  691.         hshead=Malloc(sicolors*sizeof(void *));
  692.         for (i=0; i<sicolors; ++i)
  693.             hshead[i]=Malloc(sicolors*sizeof(struct hs_t));
  694.     }
  695.     fill_wf_rgb();
  696.     fill_hsmark();
  697.     fill_hshead();
  698.     fill_hs_cm();
  699.     fill_wf2rgbweight();
  700.     wf2_redo_curcm();
  701. }
  702.  
  703.  
  704. void
  705. fill_hsmark(void)
  706. {
  707.     short i;
  708.     char *q;
  709.  
  710.     q=hsmark;
  711.     i=icolors-1;
  712.     do {
  713.         *q++ =1;
  714.     } while (--i>=0);
  715. }
  716.  
  717.  
  718. int
  719. cmphs(struct hs_t *p1, struct hs_t *p2)
  720. {
  721.     return p1->dist - p2->dist;
  722. }
  723.  
  724.  
  725. void
  726. fill_hshead(void)
  727. {
  728.     REG char *p, *q;
  729.     REG struct hs_t *hp;
  730.     NON_REG int ic;
  731.     REG short i;
  732.     REG char r, g, b, t, dist;
  733.  
  734.     p=wfrgb+3*icolors;
  735.     ic=icolors;
  736.     while (p>wfrgb) {
  737.         b= *--p;
  738.         g= *--p;
  739.         r= *--p;
  740.         q=wfrgb+3*icolors;
  741.         hp=hshead[--ic];
  742.         i=icolors-1;
  743.         do {
  744.             if ((dist=b- *--q)<0) dist= -dist;
  745.             if ((t=g- *--q)<0)    dist-=t; else dist+=t;
  746.             if ((t=r- *--q)<0)    dist-=t; else dist+=t;
  747.             hp->dist=dist;
  748.             hp->color=i;
  749.             ++hp;
  750.         } while (--i>=0);
  751.         qsort((char *)hshead[ic], (int)icolors, sizeof(struct hs_t),
  752.      (int (*)(void *, void *))cmphs);
  753.     }
  754. }
  755.  
  756.  
  757. void
  758. fill_hs_cm(void)
  759. {
  760.     NON_REG short j;
  761.     REG short i, k, bestk, colperreg;
  762.     REG long dist, bestdist;
  763.     REG struct hs_t *p;
  764.     REG char *q, *r;
  765.  
  766.     colperreg=icolors/nregs;
  767.     q=hsmark;
  768.     for (j=0; j<nregs; ++j) {
  769.         if (j==nregs-1)
  770.             colperreg+=icolors-nregs*colperreg;
  771.         bestdist=1000000000L;
  772.         for (k=0; k<icolors; ++k) {
  773.             dist=0;
  774.             p=hshead[k];
  775.             i=colperreg-1;
  776.             do {
  777.                 while (!q[p->color])
  778.                     ++p;
  779.                 dist+=p->dist;
  780.                 ++p;
  781.             } while (--i>=0);
  782.             if (dist<bestdist) {
  783.                 bestdist=dist;
  784.                 bestk=k;
  785.             }
  786.         }
  787.         p=hshead[bestk];
  788.         i=colperreg-1;
  789.         do {
  790.             while (!q[p->color])
  791.                 ++p;
  792.             q[p->color]=0;
  793.             wfminind[p->color]=j;
  794.             ++p;
  795.         } while (--i>=0);
  796.         r=wfrgb+3*bestk;
  797.         curcm[j]=(*r<<8)|(*(r+1)<<4)|*(r+2);
  798.     }
  799. }
  800.  
  801.  
  802. #define CONTRACT2
  803.  
  804. void
  805. contraction_distr(void)
  806. {
  807.     static int sicolors= -1;
  808.     int i;
  809.  
  810.     if (sicolors!= -1) {
  811.         for (i=0; i<sicolors; ++i)
  812.             if (cthead[i])
  813.                 free(cthead[i]);
  814.         free(cthead);
  815.         free(ctrgbw);
  816.     }
  817.     sicolors=icolors;
  818.     ctrgbw=Malloc(sicolors*4*sizeof(float));
  819.     cthead=Malloc(sicolors*sizeof(void *));
  820.     cthead[0]=NULL;
  821.     for (i=1; i<sicolors; ++i)
  822.         cthead[i]=Malloc(i*sizeof(float));
  823.     fill_ctrgbw();
  824.     fill_cthead();
  825.     do_contraction();
  826.     fill_ct_cm();
  827. }
  828.  
  829.  
  830. void
  831. fill_ctrgbw(void)
  832. {
  833.     short i;
  834.     u_long th, *p;
  835.     float *q;
  836.  
  837.     p=counts+NRGB;
  838.     th=threshold;
  839.     q=ctrgbw;
  840.     i=NRGB-1;
  841.     do {
  842.         if (*--p >=th) {
  843.             *q++ =(float)(i>>8);
  844.             *q++ =(float)((i>>4)&0x0f);
  845.             *q++ =(float)(i&0x0f);
  846.             *q++ =(float)*p;
  847.         }
  848.     } while (--i>=0);
  849.     assert(q==ctrgbw+4*icolors);
  850. }
  851.  
  852.  
  853. void
  854. fill_cthead(void)
  855. {
  856.     short i, j;
  857.     float *p, *q;
  858.     float r, g, b, dist, t;
  859.  
  860.     for (i=1; i<icolors; ++i) {
  861.         p=cthead[i];
  862.         q=ctrgbw+4*i;
  863.         r= *q++;
  864.         g= *q++;
  865.         b= *q;
  866.         q=ctrgbw;
  867. #ifdef CONTRACT2
  868.         for (j=0; j<i; ++j) {
  869.             t=r- *q++; dist=t*t;
  870.             t=g- *q++; dist+=t*t;
  871.             t=b- *q++; dist+=t*t;
  872.             ++q;
  873.             *p++ =dist;
  874.         }
  875. #else
  876.         for (j=0; j<i; ++j) {
  877.             if ((dist=r- *q++)<0) dist= -dist;
  878.             if ((t=g- *q++)<0)    dist-=t; else dist+=t;
  879.             if ((t=b- *q++)<0)    dist-=t; else dist+=t;
  880.             ++q;
  881.             *p++ =dist;
  882.         }
  883. #endif
  884.     }
  885. }
  886.  
  887.  
  888. void
  889. do_contraction(void)
  890. {
  891.     short ncolors, i, j, besti, bestj;
  892.     float *p, *q, bestdist;
  893.     float newr, newg, newb, neww, wi, wj;
  894.     float t, dist;
  895.  
  896.     ncolors=icolors;
  897.     while (ncolors>nregs) {
  898.         /* find minimum pair */
  899.         bestdist=2000.0;
  900.         for (i=0; i<icolors; ++i) {
  901.             p=cthead[i];
  902.             if (!p)
  903.                 continue;
  904.             j=0;
  905.             do {
  906.                 if (*p++ <bestdist) {
  907.                     besti=i;
  908.                     bestj=j;
  909.                     bestdist= *(p-1);
  910.                 }
  911.             } while (++j<i);
  912.         }
  913.         /* contract pair to new pair */
  914.         p=ctrgbw+4*besti;
  915.         q=ctrgbw+4*bestj;
  916.         wi=p[3];
  917.         wj=q[3];
  918.         neww=wi+wj;
  919.         newr=(wi*p[0]+wj*q[0])/neww;
  920.         newg=(wi*p[1]+wj*q[1])/neww;
  921.         newb=(wi*p[2]+wj*q[2])/neww;
  922.         /* delete 1 of pair */
  923.         if (bestj>besti) {
  924.             i=besti;
  925.             besti=bestj;
  926.             bestj=i;
  927.         }
  928.         free(cthead[besti]);
  929.         cthead[besti]=NULL;
  930.         ctrgbw[4*besti+3]= -1.0;
  931.         for (i=besti+1; i<icolors; ++i) {
  932.             p=cthead[i];
  933.             if (p)
  934.                 p[besti]=1000.0;
  935.         }
  936.         /* replace other of pair by new contracted color */
  937.         p=ctrgbw+4*bestj;
  938.         *p++ =newr;
  939.         *p++ =newg;
  940.         *p++ =newb;
  941.         *p   =neww;
  942.         for (i=bestj+1, q=ctrgbw+4*i; i<icolors; ++i) {
  943.             p=cthead[i];
  944.             if (p) {
  945. #ifdef CONTRACT2
  946.                 t=newr- *q++; dist=t*t;
  947.                 t=newg- *q++; dist+=t*t;
  948.                 t=newb- *q++; dist+=t*t;
  949. #else
  950.                 if ((dist=newr- *q++)<0) dist= -dist;
  951.                 if ((t=newg- *q++)<0)    dist-=t; else dist+=t;
  952.                 if ((t=newb- *q++)<0)    dist-=t; else dist+=t;
  953. #endif
  954.                 ++q;
  955.                 p[bestj]=dist;
  956.             } else {
  957.                 q+=4;
  958.             }
  959.         }
  960.         for (i=0, p=cthead[bestj]; i<bestj; ++i) {
  961.             q=ctrgbw+4*i;
  962.             if (q[3]< -0.5) {
  963.                 *p++ =1000.0;
  964.             } else {
  965. #ifdef CONTRACT2
  966.                 t=newr- *q++; dist=t*t;
  967.                 t=newg- *q++; dist+=t*t;
  968.                 t=newb- *q;   dist+=t*t;
  969. #else
  970.                 if ((dist=newr- *q++)<0) dist= -dist;
  971.                 if ((t=newg- *q++)<0)    dist-=t; else dist+=t;
  972.                 if ((t=newb- *q)<0)      dist-=t; else dist+=t;
  973. #endif
  974.                 *p++ =dist;
  975.             }
  976.         }
  977.         /* decrease number of colors */
  978.         --ncolors;
  979.     }
  980. }
  981.  
  982.  
  983. void
  984. fill_ct_cm(void)
  985. {
  986.     float *p;
  987.     u_short *q, color;
  988.     short i;
  989.  
  990.     p=ctrgbw;
  991.     q=curcm;
  992.     i=icolors-1;
  993.     do {
  994.         if (p[3]< -0.5) {
  995.             p+=4;
  996.         } else {
  997.             color =((int)(*p++ + 0.5)<<8)&0xf00;
  998.             color|=((int)(*p++ + 0.5)<<4)&0x0f0;
  999.             color|=((int)(*p++ + 0.5))   &0x00f;
  1000.             ++p;
  1001.             *q++ =color;
  1002.         }
  1003.     } while (--i>=0);
  1004.     assert(q==curcm+nregs);
  1005. #ifdef SHAM_BLACK
  1006.     if (slicedo)
  1007.     black_darkest();
  1008. #endif
  1009. }
  1010.  
  1011.  
  1012. void
  1013. fill_cmrgb(void)
  1014. {
  1015.     short i, color;
  1016.     u_short *p;
  1017.     char *q;
  1018.  
  1019.     if (cmrgb==NULL1)
  1020.     cmrgb=Malloc(3*nregs);
  1021.     p=curcm;
  1022.     q=cmrgb;
  1023.     i=nregs-1;
  1024.     do {
  1025.         color= *p++;
  1026.         *q++ = color>>8;
  1027.         *q++ = (color>>4)&0x0f;
  1028.         *q++ = color&0x0f;
  1029.     } while (--i>=0);
  1030. }
  1031.  
  1032.  
  1033. void
  1034. fill_newcol(void)
  1035. {
  1036.     short i, mini;
  1037.     char r, g, b, min, dist, t;
  1038.     NON_REG short savei;
  1039.     u_long *countp;
  1040.     char *rgbp;
  1041.  
  1042.     fill_cmrgb();
  1043.     if (newcol==NULL1) {
  1044.     newcol=Malloc(NRGB);
  1045.     error=Malloc(NRGB);
  1046.     error2=Malloc(NRGB*sizeof(u_long));
  1047.     }
  1048.     i=NRGB-1;
  1049.     countp=counts+NRGB;
  1050.     do {
  1051.         if (*--countp) {
  1052.             savei=i;
  1053.             r=i>>8;
  1054.             g=(i>>4)&0x0f;
  1055.             b=i&0x0f;
  1056.             min=100;
  1057.             rgbp=cmrgb+3*nregs;
  1058.             i=nregs-1;
  1059.             do {
  1060.                 if ((dist=b- *--rgbp)<0) dist= -dist;
  1061.                 if ((t=g- *--rgbp)<0)    dist-=t; else dist+=t;
  1062.                 if ((t=r- *--rgbp)<0)    dist-=t; else dist+=t;
  1063.                 if (dist<min) {
  1064.                     min=dist;
  1065.                     mini=i;
  1066.                 }
  1067.             } while (--i>=0);
  1068.             i=savei;
  1069.             newcol[i]=mini;
  1070.             error[i]=min;
  1071.         rgbp=cmrgb+3*mini;
  1072.         if ((t=r- *rgbp++)<0) t= -t; error2[i]=squares[t];
  1073.         if ((t=g- *rgbp++)<0) t= -t; error2[i]+=squares[t];
  1074.         if ((t=b- *rgbp  )<0) t= -t; error2[i]+=squares[t];
  1075.         }
  1076.     } while (--i>=0);
  1077. }
  1078.  
  1079.  
  1080. void
  1081. fill_newcol_count(void)
  1082. {
  1083.     short i, mini;
  1084.     char r, g, b, min, dist, t;
  1085.     NON_REG short savei;
  1086.     u_long *countp;
  1087.     char *rgbp;
  1088.  
  1089.     fill_cmrgb();
  1090.     if (newcol==NULL1) {
  1091.     newcol=Malloc(NRGB);
  1092.     error=Malloc(NRGB);
  1093.     error2=Malloc(NRGB*sizeof(u_long));
  1094.     }
  1095.     i=NRGB-1;
  1096.     countp=counts+NRGB;
  1097.     do {
  1098.         if (*--countp) {
  1099.             savei=i;
  1100.             r=i>>8;
  1101.             g=(i>>4)&0x0f;
  1102.             b=i&0x0f;
  1103.             min=100;
  1104.             rgbp=cmrgb+3*nregs;
  1105.             i=nregs-1;
  1106.             do {
  1107.                 if ((dist=b- *--rgbp)<0) dist= -dist;
  1108.                 if ((t=g- *--rgbp)<0)    dist-=t; else dist+=t;
  1109.                 if ((t=r- *--rgbp)<0)    dist-=t; else dist+=t;
  1110.                 if (dist<min) {
  1111.                     min=dist;
  1112.                     mini=i;
  1113.                 }
  1114.             } while (--i>=0);
  1115.             i=savei;
  1116.             newcol[i]=mini;
  1117.             error[i]=min;
  1118.         }
  1119.     } while (--i>=0);
  1120. }
  1121.  
  1122.  
  1123. void
  1124. restrict_cm()
  1125. {
  1126.     /* cm: array with all the colormaps, each nregs u_shorts
  1127.      * curcm: pointer into cm, pointing to the unrestricted
  1128.      *  colormap for the current line (or two lines in LACE)
  1129.      * curcm-nregs: pointer into cm, pointing to the colormap
  1130.      *  for the previous line
  1131.      * newmap: the restricted map for the current line that is
  1132.      *  produced by this function, curcm is overwritten with newmap
  1133.      * restrict_changes: max. nr. of registers may be different in
  1134.      *  curcm-nregs and the restricted curcm
  1135.      * reg_dist: reg_dist[oldn][newn] is the (manhattan) distance
  1136.      *  between (curcm-nregs)[oldn] and curcm[newn]
  1137.      * oldmin: oldmin[oldn] contains the minimum of reg_dist[oldn][x]
  1138.      * newmin: newmin[newn] contains the minimum of reg_dist[x][newn]
  1139.      * newrgb: newrgb[newn] contains r,g,b of curcm[newn]
  1140.      */
  1141.     u_short newmap[MAXNREGS];
  1142.     u_short *oldp, *newp, color;
  1143.     short i, j, ntodo, nregs1;
  1144.     struct reg_min {
  1145.         char reg;
  1146.         char dist;
  1147.     } oldmin[MAXNREGS], newmin[MAXNREGS], *minp;
  1148.     struct reg_rgb {
  1149.         char r, g, b;
  1150.     } newrgb[MAXNREGS], *rgbp;
  1151.     char *reg_dist[MAXNREGS], *charp;
  1152.  
  1153.     if (curcm==cm)    /* the first line is handled by CMAP */
  1154.         return;
  1155.     nregs1=nregs-1;
  1156.     /* mark all registers in newmap unoccupied */
  1157.     i=nregs1;
  1158.     newp=newmap;
  1159.     do {
  1160.         *newp++ =0xffff;
  1161.     } while (--i>=0);
  1162.     /* handle the colors in curcm that were already present in curcm-nregs */
  1163.     newp=curcm;
  1164.     ntodo=nregs;
  1165.     i=nregs1;
  1166.     do {
  1167.         color= *newp;
  1168.         oldp=curcm-nregs;
  1169.         j=nregs1;
  1170.         do {
  1171.             if ((*oldp & 0x0fff) == color) {
  1172.                 *oldp|=0x8000;
  1173.                 *newp|=0x8000;
  1174.                 newmap[nregs1-j]=color;
  1175.                 --ntodo;
  1176.                 break;
  1177.             }
  1178.             ++oldp;
  1179.         } while (--j>=0);
  1180.         ++newp;
  1181.     } while (--i>=0);
  1182.     /* handle at most restrict_changes in curcm that were not present in
  1183.      * curcm-nregs by searching min(ntodo, restrict_changes) times a pair
  1184.      * of regs oldreg,newreg in curcm-nregs and curcm such that oldreg
  1185.      * is not marked keep yet, newreg is not marked done yet, oldreg
  1186.      * has the largest minumum distance to the regs not done yet in curcm
  1187.      * (so we will overwrite the most useless register first),
  1188.      * newreg has the largest minimum distance to the regs marked keep
  1189.      * in curcm-nregs (so we will do the most interesting register first).
  1190.      */
  1191.     /* compute min(todo, restrict_changes) */
  1192.     if (restrict_changes<ntodo) {
  1193.         nrestricted+=ntodo-restrict_changes;
  1194.         ntodo=restrict_changes;
  1195.     }
  1196.     /* initialization of distance and min.distance tables */
  1197.     reg_dist[0]=NULL;
  1198.     if (ntodo) {
  1199.         /* allocate reg_dist */
  1200.         for (i=0; i<nregs; ++i)
  1201.             reg_dist[i]=Malloc(nregs);
  1202.         /* fill newrgb[newn] with the r,g,b values */
  1203.         newp=curcm;
  1204.         rgbp=newrgb;
  1205.         i=nregs1;
  1206.         do {
  1207.             color= *newp++;
  1208.             if (color & 0x8000)
  1209.                 rgbp->g= -1;
  1210.             else {
  1211.                 rgbp->r=color>>8;
  1212.                 rgbp->g=(color>>4)&0x0f;
  1213.                 rgbp->b=color&0x0f;
  1214.             }
  1215.             ++rgbp;
  1216.         } while (--i>=0);
  1217.         /* compute reg_dist[oldn][newn] */
  1218.         for (i=0; i<nregs; ++i) {
  1219.             char r, g, b, dist, t;
  1220.  
  1221.             color=(curcm-nregs)[i]&0x0fff;
  1222.             r=color>>8;
  1223.             g=(color>>4)&0x0f;
  1224.             b=color&0x0f;
  1225.             charp=reg_dist[i];
  1226.             rgbp=newrgb;
  1227.             j=nregs1;
  1228.             do {
  1229.                 if (rgbp->g == -1)
  1230.                     *charp++ =100;    /* not used */
  1231.                 else {
  1232.                     if ((dist=r-rgbp->r)<0) dist= -dist;
  1233.                     if ((t=g-rgbp->g)<0)    dist-=t; else dist+=t;
  1234.                     if ((t=b-rgbp->b)<0)    dist-=t; else dist+=t;
  1235.                     *charp++ =dist;
  1236.                 }
  1237.                 ++rgbp;
  1238.             } while (--j>=0);
  1239.         }
  1240.         /* fill oldmin[oldn] with minimum of reg_dist[oldn][x] */
  1241.         minp=oldmin;
  1242.         for (i=0; i<nregs; ++i) {
  1243.             if (newmap[i]==0xffff) {
  1244.             /* implies that (curcm-nregs)[i] isn't marked keep */
  1245.                 char minval=100;
  1246.                 int mini;
  1247.  
  1248.                 charp=reg_dist[i];
  1249.                 j=nregs1;
  1250.                 do {
  1251.                     if (*charp<minval) {
  1252.                         minval= *charp;
  1253.                         mini=nregs1-j;
  1254.                     }
  1255.                     ++charp;
  1256.                 } while (--j>=0);
  1257.                 minp->reg=mini;
  1258.                 minp->dist=minval;
  1259.             } else {
  1260.                 minp->dist= -1;
  1261.             }
  1262.             ++minp;
  1263.         }
  1264.         /* fill newmin[newn] with minimum of reg_dist[x][newn] */
  1265.         minp=newmin;
  1266.         rgbp=newrgb;
  1267.         for (i=0; i<nregs; ++i) {
  1268.             if (rgbp->g== -1) {
  1269.                 minp->dist= -1;
  1270.             } else {
  1271.                 char minval=100;
  1272.                 int mini;
  1273.  
  1274.                 for (j=0; j<nregs; ++j) {
  1275.                     if (newmap[j]!=0xffff
  1276.                      && reg_dist[j][i]<minval) {
  1277.                         minval=reg_dist[j][i];
  1278.                         mini=j;
  1279.                     }
  1280.                 }
  1281.                 minp->reg=mini;
  1282.                 minp->dist=minval;
  1283.             }
  1284.             ++minp;
  1285.             ++rgbp;
  1286.         }
  1287.     }
  1288.     /* one by one search for the ntodo register pairs */
  1289.     while (ntodo>0) {
  1290.         int oldreg= -1, newreg= -1;
  1291.         char maxval;
  1292.  
  1293.         --ntodo;
  1294.         /* find the oldreg */
  1295.         maxval= -1;
  1296.         minp=oldmin;
  1297.         i=nregs1;
  1298.         do {
  1299.             if (minp->dist>maxval) {
  1300.                 maxval=minp->dist;
  1301.                 oldreg=nregs1-i;
  1302.             }
  1303.             ++minp;
  1304.         } while (--i>=0);
  1305.         assert(oldreg>=0);
  1306.         /* find the newreg */
  1307.         maxval= -1;
  1308.         minp=newmin;
  1309.         i=nregs1;
  1310.         do {
  1311.             if (minp->dist>maxval) {
  1312.                 maxval=minp->dist;
  1313.                 newreg=nregs1-i;
  1314.             }
  1315.             ++minp;
  1316.         } while (--i>=0);
  1317.         assert(newreg>=0);
  1318.         /* mark newreg done, oldreg keep, and fill in value */
  1319.         color=curcm[newreg];
  1320.         assert(!(color & 0xf000));
  1321.         assert(newmap[oldreg]==0xffff);
  1322.         newmap[oldreg]=color;
  1323.         oldmin[oldreg].dist= -1;
  1324.         newmin[newreg].dist= -1;
  1325.         newrgb[newreg].g= -1;
  1326.         /* fix up oldmin */
  1327.         minp=oldmin;
  1328.         for (i=0; i<nregs; ++i) {
  1329.             if (newmap[i]==0xffff) {
  1330.                 charp=reg_dist[i];
  1331.                 charp[newreg]=100;
  1332.                 if (minp->reg==newreg) {
  1333.                     char minval=100;
  1334.                     int mini;
  1335.  
  1336.                     j=nregs1;
  1337.                     do {
  1338.                         if (*charp<minval) {
  1339.                             minval= *charp;
  1340.                             mini=nregs1-j;
  1341.                         }
  1342.                         ++charp;
  1343.                     } while (--j>=0);
  1344.                     minp->reg=mini;
  1345.                     minp->dist=minval;
  1346.                 }
  1347.             }
  1348.             ++minp;
  1349.         }
  1350.         /* fix up one row of reg_dist */
  1351.         {
  1352.             char r, g, b, dist, t;
  1353.  
  1354.             color=newmap[oldreg];
  1355.             r=color>>8;
  1356.             g=(color>>4)&0x0f;
  1357.             b=color&0x0f;
  1358.             charp=reg_dist[oldreg];
  1359.             rgbp=newrgb;
  1360.             j=nregs1;
  1361.             do {
  1362.                 if (rgbp->g == -1)
  1363.                     *charp++ =100;    /* not used */
  1364.                 else {
  1365.                     if ((dist=r-rgbp->r)<0) dist= -dist;
  1366.                     if ((t=g-rgbp->g)<0)    dist-=t; else dist+=t;
  1367.                     if ((t=b-rgbp->b)<0)    dist-=t; else dist+=t;
  1368.                     *charp++ =dist;
  1369.                 }
  1370.                 ++rgbp;
  1371.             } while (--j>=0);
  1372.         }
  1373.         /* fix up newmin */
  1374.         minp=newmin;
  1375.         rgbp=newrgb;
  1376.         for (i=0; i<nregs; ++i) {
  1377.             if (rgbp->g!= -1) {
  1378.                 if (reg_dist[oldreg][i]<minp->dist) {
  1379.                     minp->dist=reg_dist[oldreg][i];
  1380.                     minp->reg=oldreg;
  1381.                 }
  1382.             }
  1383.             ++minp;
  1384.             ++rgbp;
  1385.         }
  1386.     }
  1387.     /* clean up distance table */
  1388.     if (reg_dist[0]) {
  1389.         for (i=0; i<nregs; ++i)
  1390.             free(reg_dist[i]);
  1391.     }
  1392.     /* the still unoccupied regs in newmap will have to be filled with
  1393.      * the same value as in curcm-nregs, strip hi bits from curcm-nregs
  1394.      */
  1395.     i=nregs1;
  1396.     newp=newmap;
  1397.     oldp=curcm-nregs;
  1398.     do {
  1399.         *oldp&=0x0fff;
  1400.         if (*newp & 0xf000)
  1401.             *newp++ = *oldp++;
  1402.         else {
  1403.             ++newp;
  1404.             ++oldp;
  1405.         }
  1406.     } while (--i>=0);
  1407.     /* overwrite curcm with newmap */
  1408.     i=nregs1;
  1409.     newp=curcm;
  1410.     oldp=newmap;
  1411.     do {
  1412.         *newp++ = *oldp++;
  1413.     } while (--i>=0);
  1414. }
  1415.