home *** CD-ROM | disk | FTP | other *** search
/ vsiftp.vmssoftware.com / VSIPUBLIC@vsiftp.vmssoftware.com.tar / FREEWARE / FREEWARE40.ZIP / xpaint-247 / iprocess.c < prev    next >
C/C++ Source or Header  |  1997-01-03  |  36KB  |  1,524 lines
  1. /* +-------------------------------------------------------------------+ */
  2. /* | Copyright 1993, David Koblas (koblas@netcom.com)               | */
  3. /* | Copyright 1995, 1996 Torsten Martinsen (bullestock@dk-online.dk)  | */
  4. /* |                                       | */
  5. /* | Permission to use, copy, modify, and to distribute this software  | */
  6. /* | and its documentation for any purpose is hereby granted without   | */
  7. /* | fee, provided that the above copyright notice appear in all       | */
  8. /* | copies and that both that copyright notice and this permission    | */
  9. /* | notice appear in supporting documentation.     There is no           | */
  10. /* | representations about the suitability of this software for           | */
  11. /* | any purpose.  this software is provided "as is" without express   | */
  12. /* | or implied warranty.                           | */
  13. /* |                                       | */
  14. /* +-------------------------------------------------------------------+ */
  15.  
  16. /* $Id: iprocess.c,v 1.11 1996/05/09 07:12:23 torsten Exp $ */
  17.  
  18. #include <X11/Intrinsic.h>
  19. #include <X11/StringDefs.h>
  20. #include <X11/Xatom.h>
  21. #include <stdio.h>
  22. #include <stdlib.h>
  23. #include <math.h>
  24. #include "xpaint.h"
  25. #include "image.h"
  26. #include "misc.h"
  27. #include "protocol.h"
  28.  
  29. #ifndef M_PI
  30. #define M_PI            3.14159265358979323846
  31. #endif
  32.  
  33. /*
  34. **  Some real image processing
  35.  */
  36.  
  37. extern struct imageprocessinfo ImgProcessInfo;
  38.  
  39. #define CLAMP(low, value, high) \
  40.         if (value < low) value = low; else if (value > high) value = high
  41. #define DEG2RAD    (M_PI/180.0)
  42.  
  43. #define ConvMatrixSize    3
  44. typedef float ConvMatrix[ConvMatrixSize * ConvMatrixSize];
  45.  
  46. static int *histogram(Image * input);
  47.  
  48. static Image *
  49. convolve(Image * input, float *mat, int n,
  50.      unsigned char *basePixel, Boolean absFlag)
  51. {
  52.     int x, y, xx, yy, xv, yv;
  53.     float sum;
  54.     unsigned char *p;
  55.     unsigned char *op;
  56.     float r, g, b;
  57.     int ir, ig, ib;
  58.     Image *output;
  59.  
  60.     output = ImageNew(input->width, input->height);
  61.     op = output->data;
  62.  
  63.     for (y = 0; y < input->height; y++) {
  64.     for (x = 0; x < input->width; x++) {
  65.         r = g = b = 0;
  66.         sum = 0;
  67.         for (yy = 0; yy < n; yy++) {
  68.         for (xx = 0; xx < n; xx++) {
  69.             xv = x + xx - 1;
  70.             yv = y + yy - 1;
  71.             if (xv < 0 || yv < 0)
  72.             continue;
  73.             if (xv >= input->width || yv >= input->height)
  74.             continue;
  75.             p = ImagePixel(input, xv, yv);
  76.             r += (float) *p++ * mat[xx + n * yy];
  77.             g += (float) *p++ * mat[xx + n * yy];
  78.             b += (float) *p * mat[xx + n * yy];
  79.             sum += mat[xx + n * yy];
  80.         }
  81.         }
  82.         if (sum <= 0)
  83.         sum = 1;
  84.         if (absFlag) {
  85.         if (r < 0)
  86.             r = -r;
  87.         if (g < 0)
  88.             g = -g;
  89.         if (b < 0)
  90.             b = -b;
  91.         }
  92.         ir = r / sum;
  93.         ig = g / sum;
  94.         ib = b / sum;
  95.         if (basePixel) {
  96.         ir += basePixel[0];
  97.         ig += basePixel[1];
  98.         ib += basePixel[2];
  99.         }
  100.         CLAMP(0, ir, 255);
  101.         CLAMP(0, ig, 255);
  102.         CLAMP(0, ib, 255);
  103.  
  104.         *op++ = ir;
  105.         *op++ = ig;
  106.         *op++ = ib;
  107.     }
  108.  
  109.     if (y % 16 == 0)
  110.         StateTimeStep();
  111.     }
  112.  
  113.     return output;
  114. }
  115.  
  116. /*
  117. **  rescale values from 0..255
  118.  */
  119. static void 
  120. normalize(Image * image)
  121. {
  122.     int i, count;
  123.     unsigned char *sp;
  124.     unsigned char *ip;
  125.     int maxval = 0;
  126.  
  127.     if (image->cmapSize != 0) {
  128.     sp = image->cmapData;
  129.     count = image->cmapSize * 3;
  130.     } else {
  131.     sp = image->data;
  132.     count = image->width * image->height * 3;
  133.     }
  134.  
  135.     for (ip = sp, i = 0; i < count; i++, ip++)
  136.     if (*ip > maxval)
  137.         maxval = *ip;
  138.     if (maxval == 0)
  139.     return;
  140.     for (ip = sp, i = 0; i < count; i++, ip++)
  141.     *ip = ((int) *ip * 255) / maxval;
  142. }
  143.  
  144. /*
  145. **  Convert image into a monochrome image
  146.  */
  147. static void 
  148. monochrome(Image * image)
  149. {
  150.     int i, count;
  151.     unsigned char *sp;
  152.     unsigned char *ip;
  153.     int v;
  154.  
  155.     if (image->cmapSize != 0) {
  156.     sp = image->cmapData;
  157.     count = image->cmapSize;
  158.     } else {
  159.     sp = image->data;
  160.     count = image->width * image->height;
  161.     }
  162.  
  163.     for (ip = sp, i = 0; i < count; i++) {
  164.     v = (ip[0] * 11 + ip[1] * 16 + ip[2] * 5) >> 5;        /* pp = .33R+.5G+.17B */
  165.     *ip++ = v;
  166.     *ip++ = v;
  167.     *ip++ = v;
  168.     }
  169. }
  170.  
  171.  
  172. Image *
  173. ImageSmooth(Image * input)
  174. {
  175.     float *mat;
  176.     int n, i;
  177.     Image *output;
  178.  
  179.     /*
  180.      * Build n x n convolution matrix (all 1's)
  181.      */
  182.     n = ImgProcessInfo.smoothMaskSize;
  183.     mat = malloc(n * n * sizeof(float));
  184.     for (i = 0; i < n * n; ++i)
  185.     mat[i] = 1.0;
  186.  
  187.     output = convolve(input, mat, n, NULL, False);
  188.     free(mat);
  189.     return output;
  190. }
  191.  
  192. Image *
  193. ImageSharpen(Image * input)
  194. {
  195.     static ConvMatrix mat =
  196.     {
  197.     -1, -2, -1,
  198.     -2, 20, -2,
  199.     -1, -2, -1
  200.     };
  201.     return convolve(input, mat, ConvMatrixSize, NULL, False);
  202. }
  203.  
  204. Image *
  205. ImageEdge(Image * input)
  206. {
  207.     static ConvMatrix mat =
  208.     {
  209.     -1, -2, 0,
  210.     -2, 0, 2,
  211.     0, 2, 1
  212.     };
  213.     Image *image = convolve(input, mat, ConvMatrixSize, NULL, True);
  214.  
  215.     normalize(image);
  216.  
  217.     return image;
  218. }
  219.  
  220. Image *
  221. ImageEmbose(Image * input)
  222. {
  223.     static ConvMatrix mat =
  224.     {
  225.     -1, -2, 0,
  226.     -2, 0, 2,
  227.     0, 2, 1
  228.     };
  229.     static unsigned char base[3] =
  230.     {128, 128, 128};
  231.     Image *image = convolve(input, mat, ConvMatrixSize, base, False);
  232.  
  233.     monochrome(image);
  234.     normalize(image);
  235.  
  236.     return image;
  237. }
  238.  
  239. Image *
  240. ImageInvert(Image * input)
  241. {
  242.     Image *output = ImageNewCmap(input->width, input->height, input->cmapSize);
  243.     int i;
  244.     unsigned char *ip, *op;
  245.     int count;
  246.  
  247.     /*
  248.     **     If the input has a colormap, just invert that.
  249.      */
  250.     if (input->cmapSize != 0) {
  251.     ip = input->cmapData;
  252.     op = output->cmapData;
  253.     count = input->cmapSize;
  254.  
  255.     memcpy(output->data, input->data,
  256.          sizeof(char) * input->scale * input->width * input->height);
  257.     } else {
  258.     ip = input->data;
  259.     op = output->data;
  260.     count = input->width * input->height;
  261.     }
  262.  
  263.     for (i = 0; i < count; i++) {
  264.     *op++ = 255 - *ip++;
  265.     *op++ = 255 - *ip++;
  266.     *op++ = 255 - *ip++;
  267.     }
  268.  
  269.     return output;
  270. }
  271.  
  272. Image *
  273. ImageOilPaint(Image * input)
  274. {
  275.     Image *output = ImageNew(input->width, input->height);
  276.     unsigned char *op = output->data;
  277.     int x, y, xx, yy, i;
  278.     int rVal, gVal, bVal;
  279.     int rCnt, gCnt, bCnt;
  280.     int rHist[256], gHist[256], bHist[256];
  281.     int oilArea_2;
  282.  
  283.  
  284.     oilArea_2 = ImgProcessInfo.oilArea / 2;
  285.  
  286.     for (y = 0; y < input->height; y++) {
  287.     for (x = 0; x < input->width; x++) {
  288.         /*
  289.         **    compute histogram of (on-screen hunk of) n*n 
  290.         **      region centered plane
  291.          */
  292.  
  293.         rCnt = gCnt = bCnt = 0;
  294.         rVal = gVal = bVal = 0;
  295.         for (i = 0; i < XtNumber(rHist); i++)
  296.         rHist[i] = gHist[i] = bHist[i] = 0;
  297.  
  298.         for (yy = y - oilArea_2; yy < y + oilArea_2; yy++) {
  299.         if (yy < 0 || yy >= input->height)
  300.             continue;
  301.         for (xx = x - oilArea_2; xx < x + oilArea_2; xx++) {
  302.             int c, p;
  303.             unsigned char *rgb;
  304.  
  305.             if (xx < 0 || xx >= input->width)
  306.             continue;
  307.  
  308.             rgb = ImagePixel(input, xx, yy);
  309.  
  310.             if ((c = ++rHist[(p = rgb[0]) / 4]) > rCnt) {
  311.             rVal = p;
  312.             rCnt = c;
  313.             }
  314.             if ((c = ++gHist[(p = rgb[1]) / 4]) > gCnt) {
  315.             gVal = p;
  316.             gCnt = c;
  317.             }
  318.             if ((c = ++bHist[(p = rgb[2]) / 4]) > bCnt) {
  319.             bVal = p;
  320.             bCnt = c;
  321.             }
  322.         }
  323.         }
  324.  
  325.         *op++ = rVal;
  326.         *op++ = gVal;
  327.         *op++ = bVal;
  328.     }
  329.  
  330.     if (y % 16 == 0)
  331.         StateTimeStep();
  332.     }
  333.  
  334.     return output;
  335. }
  336.  
  337.  
  338. /*
  339.  * Add random noise to an image.
  340.  */
  341. Image *
  342. ImageAddNoise(Image * input)
  343. {
  344.     unsigned char *p;
  345.     int x, y, r, g, b, ddelta;
  346.     Image *output = ImageNew(input->width, input->height);
  347.     unsigned char *op = output->data;
  348.  
  349.  
  350.     ddelta = ImgProcessInfo.noiseDelta * 2;
  351.  
  352.     for (y = 0; y < input->height; y++) {
  353.     for (x = 0; x < input->width; x++) {
  354.         p = ImagePixel(input, x, y);
  355.         r = p[0] + ddelta * gauss();
  356.         CLAMP(0, r, 255);
  357.         g = p[1] + ddelta * gauss();
  358.         CLAMP(0, g, 255);
  359.         b = p[2] + ddelta * gauss();
  360.         CLAMP(0, b, 255);
  361.         *op++ = r;
  362.         *op++ = g;
  363.         *op++ = b;
  364.     }
  365.     }
  366.     return output;
  367. }
  368.  
  369. /*
  370.  * This function works in-place.
  371.  * Because it swaps pixels in the input image, which may be colour mapped
  372.  * or not, is has knowledge about the Image format that should probably
  373.  * have stayed in image.h. Too bad.
  374.  */
  375. Image *
  376. ImageSpread(Image * input)
  377. {
  378.     int w = input->width, h = input->height;
  379.     unsigned char *p, *np;
  380.     int x, y, minx, miny, maxx, maxy, xn, yn, dist;
  381.  
  382.  
  383.     dist = ImgProcessInfo.spreadDistance;
  384.     for (y = 0; y < h; y++) {
  385.     for (x = 0; x < w; x++) {
  386.         p = ImagePixel(input, x, y);
  387.         /* find random neighbour pixel within +- dist */
  388.         minx = x - dist;
  389.         if (minx < 0)
  390.         minx = 0;
  391.         maxx = x + dist;
  392.         if (maxx >= w)
  393.         maxx = w - 1;
  394.         xn = RANDOMI2(minx, maxx);
  395.  
  396.         miny = y - dist;
  397.         if (miny < 0)
  398.         miny = 0;
  399.         maxy = y + dist;
  400.         if (maxy >= h)
  401.         maxy = h - 1;
  402.         yn = RANDOMI2(miny, maxy);
  403.         /* swap pixel with neighbour */
  404.         np = ImagePixel(input, xn, yn);
  405.         if (input->cmapSize == 0) {
  406.         unsigned char rn, gn, bn;
  407.  
  408.         rn = np[0];
  409.         gn = np[1];
  410.         bn = np[2];
  411.         np[0] = p[0];
  412.         np[1] = p[1];
  413.         np[2] = p[2];
  414.         p[0] = rn;
  415.         p[1] = gn;
  416.         p[2] = bn;
  417.         } else if (input->cmapSize > 256) {
  418.         unsigned short i, *ps, *nps;
  419.  
  420.         ps = &((unsigned short *) input->data)[y * w + x];
  421.         nps = &((unsigned short *) input->data)[yn * w + xn];
  422.         i = *nps;
  423.         *nps = *ps;
  424.         *ps = i;
  425.         } else {        /* colour map of 256 or less */
  426.         unsigned char i;
  427.  
  428.         p = &input->data[y * w + x];
  429.         np = &input->data[yn * w + xn];
  430.         i = *np;
  431.         *np = *p;
  432.         *p = i;
  433.         }
  434.     }
  435.     }
  436.  
  437.     return input;
  438. }
  439.  
  440. Image *
  441. ImageBlend(Image * input)
  442. {
  443.     int w = input->width, h = input->height;
  444.     Image *output = ImageNew(w, h);
  445.     unsigned char *op = output->data;
  446.     unsigned char *rgb;
  447.     int x, y, n, ox, oy, px, py, dx, dy;
  448.     int rSum, gSum, bSum, r, g, b;
  449.     float diagonal, gradient, ap;
  450.  
  451.  
  452.     ox = w / 2;
  453.     oy = h / 2;
  454.     diagonal = ((double) h) / w;
  455.  
  456.     /*
  457.      * Compute average of pixels on edge of region. While we are
  458.      * at it, we copy the edge to the output as well.
  459.      */
  460.     rSum = gSum = bSum = 0;
  461.     for (x = 0; x < w; ++x) {
  462.     rgb = ImagePixel(input, x, 0);
  463.     r = rgb[0];
  464.     g = rgb[1];
  465.     b = rgb[2];
  466.     rSum += r;
  467.     gSum += g;
  468.     bSum += b;
  469.     *op++ = r;
  470.     *op++ = g;
  471.     *op++ = b;
  472.     }
  473.     op = ImagePixel(output, 0, h - 1);
  474.     for (x = 0; x < w; ++x) {
  475.     rgb = ImagePixel(input, x, h - 1);
  476.     r = rgb[0];
  477.     g = rgb[1];
  478.     b = rgb[2];
  479.     rSum += r;
  480.     gSum += g;
  481.     bSum += b;
  482.     *op++ = r;
  483.     *op++ = g;
  484.     *op++ = b;
  485.     }
  486.     for (y = 1; y < h - 1; ++y) {
  487.     rgb = ImagePixel(input, 0, y);
  488.     r = rgb[0];
  489.     g = rgb[1];
  490.     b = rgb[2];
  491.     rSum += r;
  492.     gSum += g;
  493.     bSum += b;
  494.     op = ImagePixel(output, 0, y);
  495.     *op++ = r;
  496.     *op++ = g;
  497.     *op++ = b;
  498.  
  499.     rgb = ImagePixel(input, w - 1, y);
  500.     r = rgb[0];
  501.     g = rgb[1];
  502.     b = rgb[2];
  503.     rSum += r;
  504.     gSum += g;
  505.     bSum += b;
  506.     op = ImagePixel(output, w - 1, y);
  507.     *op++ = r;
  508.     *op++ = g;
  509.     *op++ = b;
  510.     }
  511.  
  512.     n = 2 * (w + h - 2);    /* total # of pixels on edge */
  513.     rSum /= n;
  514.     gSum /= n;
  515.     bSum /= n;
  516.     op = ImagePixel(output, 1, 1);
  517.  
  518.     /* compute colour of each point in the interior */
  519.     for (y = 1; y < h - 1; y++) {
  520.     for (x = 1; x < w - 1; x++) {
  521.         dx = x - ox;
  522.         dy = y - oy;
  523.         if (dx == 0 && dy == 0) {    /* this is the centre point */
  524.         r = rSum;
  525.         g = gSum;
  526.         b = bSum;
  527.         } else {
  528.         if (dx == 0) {    /* special case 1 */
  529.             px = ox;
  530.             py = (dy > 0) ? h - 1 : 0;
  531.         } else if (dy == 0) {    /* special case 2 */
  532.             py = oy;
  533.             px = (dx > 0) ? w - 1 : 0;
  534.         } else {    /* general case */
  535.             gradient = ((double) dy) / dx;
  536.             if (fabs(gradient) < fabs(diagonal)) {
  537.             px = (dx > 0) ? w - 1 : 0;
  538.             py = oy + ((px - ox) * dy) / dx;
  539.             } else {
  540.             py = (dy > 0) ? h - 1 : 0;
  541.             px = ox + ((py - oy) * dx) / dy;
  542.             }
  543.         }
  544.  
  545.         /*
  546.          * given O(ox,oy), P(px,py), and A(x,y), compute
  547.          *   |AO|/|PO|
  548.          */
  549.         ap = sqrt((double) (dx * dx) + (double) (dy * dy)) /
  550.             sqrt((double) ((px - ox) * (px - ox)) +
  551.              (double) ((py - oy) * (py - oy)));
  552.  
  553.         rgb = ImagePixel(input, px, py);
  554.         r = (1 - ap) * rSum + ap * rgb[0];
  555.         g = (1 - ap) * gSum + ap * rgb[1];
  556.         b = (1 - ap) * bSum + ap * rgb[2];
  557.         }
  558.         *op++ = r;
  559.         *op++ = g;
  560.         *op++ = b;
  561.     }
  562.     op += 6;        /* skip last on this line and first on next */
  563.     }
  564.  
  565.     return output;
  566. }
  567.  
  568. Image *
  569. ImagePixelize(Image * input)
  570. {
  571.     int width = input->width, height = input->height;
  572.     Image *output = ImageNew(width, height);
  573.     unsigned char *op, *rgb;
  574.     int x, y, xx, yy, n;
  575.     int rSum, gSum, bSum;
  576.     int xsize, ysize;
  577.  
  578.  
  579.     xsize = ImgProcessInfo.pixelizeXSize;
  580.     ysize = ImgProcessInfo.pixelizeYSize;
  581.  
  582.     if (xsize > width)
  583.     xsize = width;
  584.     if (ysize > height)
  585.     ysize = height;
  586.     n = xsize * ysize;
  587.  
  588.     for (y = 0; y < height; y += ysize) {
  589.     for (x = 0; x < width; x += xsize) {
  590.         /*
  591.          *    compute average of pixels inside megapixel
  592.          */
  593.         rSum = gSum = bSum = 0;
  594.         for (yy = y; yy < y + ysize; yy++) {
  595.         if (yy >= height)
  596.             continue;
  597.         for (xx = x; xx < x + xsize; xx++) {
  598.             if (xx >= width)
  599.             continue;
  600.             rgb = ImagePixel(input, xx, yy);
  601.             rSum += rgb[0];
  602.             gSum += rgb[1];
  603.             bSum += rgb[2];
  604.         }
  605.         }
  606.         rSum /= n;
  607.         gSum /= n;
  608.         bSum /= n;
  609.         /*
  610.          * replace each pixel in megapixel with average
  611.          */
  612.         for (yy = y; yy < y + ysize; yy++) {
  613.         if (yy >= height)
  614.             continue;
  615.         for (xx = x; xx < x + xsize; xx++) {
  616.             if (xx >= width)
  617.             continue;
  618.             op = ImagePixel(output, xx, yy);
  619.             *op++ = rSum;
  620.             *op++ = gSum;
  621.             *op = bSum;
  622.         }
  623.         }
  624.     }
  625.  
  626.     if (y % 16 == 0)
  627.         StateTimeStep();
  628.     }
  629.  
  630.     /* if any incomplete megapixels are left, copy them to the output */
  631.     for (x = (width / xsize) * xsize; x < width; ++x)
  632.     for (y = 0; y < height; ++y) {
  633.         rgb = ImagePixel(input, x, y);
  634.         op = ImagePixel(output, x, y);
  635.         *op++ = *rgb++;
  636.         *op++ = *rgb++;
  637.         *op = *rgb;
  638.     }
  639.     for (y = (height / ysize) * ysize; y < height; ++y)
  640.     for (x = 0; x < width; ++x) {
  641.         rgb = ImagePixel(input, x, y);
  642.         op = ImagePixel(output, x, y);
  643.         *op++ = *rgb++;
  644.         *op++ = *rgb++;
  645.         *op = *rgb;
  646.     }
  647.  
  648.     return output;
  649. }
  650.  
  651.  
  652. #define SWAP(a, b)    { int t = (a); (a) = (b); (b) = t; }
  653.  
  654. Image *
  655. ImageDespeckle(Image * input)
  656. {
  657.     int width = input->width, height = input->height;
  658.     Image *output = ImageNew(width, height);
  659.     unsigned char *op, *rgb;
  660.     int x, y, xx, yy, mask, mask2, i, j, k, l;
  661.     int *ra, *ga, *ba;
  662.  
  663.     mask = ImgProcessInfo.despeckleMask;
  664.     if (mask > width)
  665.     mask = width;
  666.     if (mask > height)
  667.     mask = height;
  668.     mask2 = mask / 2;
  669.  
  670.     /* arrays for storing pixels inside mask */
  671.     ra = xmalloc(mask * mask * sizeof(int));
  672.     ga = xmalloc(mask * mask * sizeof(int));
  673.     ba = xmalloc(mask * mask * sizeof(int));
  674.  
  675.     op = output->data;
  676.     for (y = 0; y < height; ++y) {
  677.     for (x = 0; x < width; ++x) {
  678.         i = 0;
  679.         for (yy = MAX(0, y - mask2); yy < MIN(height, y + mask2); ++yy)
  680.         for (xx = MAX(0, x - mask2); xx < MIN(width, x + mask2); ++xx) {
  681.             rgb = ImagePixel(input, xx, yy);
  682.             ra[i] = *rgb++;
  683.             ga[i] = *rgb++;
  684.             ba[i] = *rgb;
  685.             ++i;
  686.         }
  687.         /*
  688.          * now find median by (shell-)sorting the arrays and
  689.          * picking the center value
  690.          */
  691.         for (j = i / 2; j > 0; j = j / 2)
  692.         for (k = j; k < i; k++) {
  693.             for (l = k - j; l >= 0 && ra[l] > ra[l + j]; l -= j)
  694.             SWAP(ra[l], ra[l + j]);
  695.             for (l = k - j; l >= 0 && ga[l] > ga[l + j]; l -= j)
  696.             SWAP(ga[l], ga[l + j]);
  697.             for (l = k - j; l >= 0 && ba[l] > ba[l + j]; l -= j)
  698.             SWAP(ba[l], ba[l + j]);
  699.         }
  700.         if (i & 1) {    /* uneven number of data points */
  701.         *op++ = ra[i / 2];
  702.         *op++ = ga[i / 2];
  703.         *op++ = ba[i / 2];
  704.         } else {        /* even, take average */
  705.         *op++ = (ra[i / 2 - 1] + ra[i / 2]) / 2;
  706.         *op++ = (ga[i / 2 - 1] + ga[i / 2]) / 2;
  707.         *op++ = (ba[i / 2 - 1] + ba[i / 2]) / 2;
  708.         }
  709.     }
  710.     if (y % 16 == 0)
  711.         StateTimeStep();
  712.     }
  713.  
  714.     free(ra);
  715.     free(ga);
  716.     free(ba);
  717.  
  718.     return output;
  719. }
  720.  
  721.  
  722. /*
  723.  * Normalize the contrast of an image.
  724.  */
  725. Image *
  726. ImageNormContrast(Image * input)
  727. {
  728.     unsigned char *p;
  729.     int *hist;
  730.     int x, y, w, h, i, max = 0, cumsum, limit;
  731.     Image *output = ImageNew(input->width, input->height);
  732.     unsigned char *op = output->data;
  733.     int blackpcnt, blackval, whitepcnt, whiteval;
  734.  
  735.  
  736.     blackpcnt = ImgProcessInfo.contrastB;
  737.     whitepcnt = 100 - ImgProcessInfo.contrastW;
  738.  
  739.     hist = histogram(input);
  740.  
  741.     w = input->width;
  742.     h = input->height;
  743.  
  744.     /* Find lower knee point in histogram to determine 'black' threshold. */
  745.     cumsum = 0;
  746.     limit = w * h * blackpcnt / 100;
  747.     for (i = 0; i < 256; ++i)
  748.     if ((cumsum += hist[i]) > limit)
  749.         break;
  750.     blackval = i;
  751.  
  752.     /* Likewise for 'white' threshold. */
  753.     cumsum = 0;
  754.     limit = w * h * whitepcnt / 100;
  755.     for (i = 256; i > 0; --i) {
  756.     if ((max == 0) && hist[i])
  757.         max = i + 1;    /* max is index of highest non-zero entry */
  758.     if ((cumsum += hist[i]) > limit)
  759.         break;
  760.     }
  761.     whiteval = i;
  762.  
  763.     /*
  764.      * Now create a histogram with a linear ramp
  765.      * from (blackval, 0) to (whiteval, max)
  766.      */
  767.     for (i = 0; i < blackval; ++i)
  768.     hist[i] = 0;
  769.     for (; i <= whiteval; ++i)
  770.     hist[i] = max * (i - blackval) / (whiteval - blackval);
  771.     for (; i < 256; ++i)
  772.     hist[i] = max;
  773.  
  774.     for (y = 0; y < h; y++) {
  775.     for (x = 0; x < w; x++) {
  776.         p = ImagePixel(input, x, y);
  777.         *op++ = hist[*p++];
  778.         *op++ = hist[*p++];
  779.         *op++ = hist[*p];
  780.     }
  781.     }
  782.     free(hist);
  783.  
  784.     return output;
  785. }
  786.  
  787. /*
  788.  * Build histogram data for the image.
  789.  */
  790. static int *
  791. histogram(Image * input)
  792. {
  793.     unsigned char *p;
  794.     int x, y, w, h, i, r, g, b, *hist;
  795.  
  796.     w = input->width;
  797.     h = input->height;
  798.  
  799.     /* Make a table of 256 zeros plus one sentinel */
  800.     hist = xmalloc(257 * sizeof(int));
  801.     for (i = 0; i <= 256; ++i)
  802.     hist[i] = 0;
  803.  
  804.     /*
  805.      * Build histogram of intensity values.
  806.      */
  807.     for (y = 0; y < h; y++) {
  808.     for (x = 0; x < w; x++) {
  809.         p = ImagePixel(input, x, y);
  810.         r = *p++;
  811.         g = *p++;
  812.         b = *p;
  813.         i = (r * 169 + g * 256 + b * 87) / 512;    /* i = .33R+.5G+.17B */
  814.         ++hist[i];
  815.     }
  816.     }
  817.     return hist;
  818. }
  819.  
  820. #if 0
  821.  
  822. #define HIST_H      120        /* Height of histogram */
  823. #define HIST_W      256        /* Width of histogram */
  824. #define HIST_B        2        /* Horizontal border width for histogram */
  825. #define HIST_R       12        /* Height of ruler below histogram */
  826. #define HIST_T        9        /* Height of tick marks on ruler */
  827. #define LIGHTGREY 200        /* Greyscale level used for bars */
  828. #define MIDGREY      255        /* Greyscale level used for ruler */
  829. #define DARKGREY  120        /* Greyscale level used for background */
  830.  
  831. /*
  832.  * Create an image (width 256, height HIST_H) depicting the histogram data.
  833.  */
  834. Image *
  835. HistogramImage(char *hist)
  836. {
  837.     unsigned char *p;
  838.     int x, y, i;
  839.     int max, e:
  840.     unsigned char *ihist;
  841.     Image *output = ImageNewGrey(HIST_W + 2 * HIST_B, HIST_H + HIST_R);;
  842.  
  843.     max = 0;
  844.     ihist = xmalloc(128 * sizeof(int));
  845.     for (i = 0; i < 128; ++i) {
  846.     e = ihist[i] = hist[i * 2] + hist[i * 2 + 1];
  847.     if (e > max)
  848.         max = e;
  849.     }
  850.     /*
  851.      * Now create histogram image.
  852.      * Only display 128 bins to make the histogram easier to read.
  853.      * Leave a border of a few pixels in each side.
  854.      */
  855.  
  856.     memset(output->data, DARKGREY, (HIST_W + 2 * HIST_B) * (HIST_H + HIST_R));
  857.     for (x = 0; x < 128; ++x) {
  858.     i = ihist[x] * (HIST_H - 1) / max;
  859.     for (y = HIST_H - 1 - i; y < HIST_H - 1; y++) {
  860.         p = output->data + y * (HIST_W + 2 * HIST_B) + x * 2 + HIST_B;
  861.         *p++ = LIGHTGREY;
  862.         *p = LIGHTGREY;
  863.     }
  864.     }
  865.     /* Ruler */
  866.     memset(output->data + HIST_H * (HIST_W + 2 * HIST_B) + HIST_B,
  867.        MIDGREY, HIST_W);
  868.     for (x = 0; x < 5; ++x) {
  869.     int tick[] =
  870.     {0, HIST_W / 4, HIST_W / 2, HIST_W * 3 / 4, HIST_W - 1};
  871.  
  872.     for (y = 0; y < HIST_T; y++) {
  873.         p = output->data + (y + HIST_H + 1) * (HIST_W + 2 * HIST_B)
  874.         + tick[x] + HIST_B;
  875.         *p = MIDGREY;
  876.     }
  877.     }
  878.     free(ihist);
  879.  
  880.     return output;
  881. }
  882. #endif
  883.  
  884. #ifdef FEATURE_TILT
  885. /*
  886.  * Tilt an image.
  887.  */
  888.  
  889. /*
  890.  * Compute the interpolated RGB values of a 1 x 1 pixel centered
  891.  * at (x,y) and store these values in op[0] trough op[2].
  892.  * The interpolation is based on weighting the RGB values proportionally
  893.  * to the overlapping areas.
  894.  */
  895. static void 
  896. interpol(float x, float y, Image * input, unsigned char *op)
  897. {
  898.     int xl, xh, yl, yh;
  899.     float a, b, c, d, A, B, C, D;
  900.     unsigned char *plh, *phh, *pll, *phl;
  901.  
  902.     xl = floor(x);
  903.     xh = ceil(x);
  904.     yl = floor(y);
  905.     yh = ceil(y);
  906.  
  907.     pll = ImagePixel(input, xl, yl);
  908.  
  909.     if (xh == xl) {
  910.     if (yh == yl) {
  911.         /* pixel coincides with one pixel */
  912.         *op++ = *pll++;
  913.         *op++ = *pll++;
  914.         *op++ = *pll++;
  915.     } else {
  916.         /* pixel overlaps with two pixels on top of each other */
  917.         plh = ImagePixel(input, xl, yh);
  918.         A = y - yl;
  919.         B = yh - y;
  920.         *op++ = A * *plh++ + B * *pll++;
  921.         *op++ = A * *plh++ + B * *pll++;
  922.         *op++ = A * *plh++ + B * *pll++;
  923.     }
  924.     } else if (yh == yl) {
  925.     /* pixel overlaps with two side-by-side pixels */
  926.     phl = ImagePixel(input, xh, yl);
  927.     A = x - xl;
  928.     B = xh - x;
  929.     *op++ = A * *phl++ + B * *pll++;
  930.     *op++ = A * *phl++ + B * *pll++;
  931.     *op++ = A * *phl++ + B * *pll++;
  932.     } else {
  933.     /* general case: pixel overlaps all four neighbour pixels */
  934.     a = xh - x;
  935.     b = y - yl;
  936.     c = x - xl;
  937.     d = yh - y;
  938.     A = a * b;
  939.     B = b * c;
  940.     C = a * d;
  941.     D = c * d;
  942.     plh = ImagePixel(input, xl, yh);
  943.     phl = ImagePixel(input, xh, yl);
  944.     phh = ImagePixel(input, xh, yh);
  945.     *op++ = A * *plh++ + B * *phh++ + C * *pll++ + D * *phl++;
  946.     *op++ = A * *plh++ + B * *phh++ + C * *pll++ + D * *phl++;
  947.     *op++ = A * *plh++ + B * *phh++ + C * *pll++ + D * *phl++;
  948.     }
  949. }
  950.  
  951.  
  952. Image *
  953. ImageTilt(Image * input)
  954. {
  955.     int x, y, br, bg, bb;
  956.     int width = input->width, height = input->height;
  957.     Image *output;
  958.     unsigned char *op;
  959.     float xt, yt;
  960.     float X1, Y1, X2, Y2;
  961.  
  962.  
  963.     X1 = width * ImgProcessInfo.tiltX1 / 100;
  964.     Y1 = height * ImgProcessInfo.tiltY1 / 100;
  965.     X2 = width * ImgProcessInfo.tiltX2 / 100;
  966.     Y2 = height * ImgProcessInfo.tiltY2 / 100;
  967.  
  968.     if (((Y1 >= 0) && (Y1 < height)) || ((X2 >= 0) && (X2 < width)))
  969.     return input;
  970.  
  971.     output = ImageNew(width, height);
  972.  
  973.     /* Get RGB values of background colour */
  974.     br = ImgProcessInfo.background->red / 256;
  975.     bg = ImgProcessInfo.background->green / 256;
  976.     bb = ImgProcessInfo.background->blue / 256;
  977.  
  978.     op = output->data;
  979.     for (x = 0; x < width * height; ++x) {
  980.     *op++ = br;
  981.     *op++ = bg;
  982.     *op++ = bb;
  983.     }
  984.  
  985.     for (y = 0; y < height; y++) {
  986.     for (x = 0; x < width; x++) {
  987.         op = ImagePixel(output, x, y);
  988.         /* find the input coords corresponding to output coords (x,y) */
  989.         xt = (X1 * y - Y1 * x) / (y - Y1);
  990.         yt = (-X2 * y + x * Y2) / (x - X2);
  991.         if ((xt >= 0) && (xt < width) && (yt >= 0) && (yt < height)) {
  992. #if 0
  993.         unsigned char *p;
  994.         /* rounding */
  995.         p = ImagePixel(input, (int) (xt + 0.5), (int) (yt + 0.5));
  996.         *op++ = *p++;
  997.         *op++ = *p++;
  998.         *op = *p;
  999. #else
  1000.         interpol(xt, yt, input, op);
  1001. #endif
  1002.         }
  1003.     }
  1004.     }
  1005.  
  1006.     return output;
  1007. }
  1008. #endif
  1009.  
  1010. /*
  1011.  * Produce the 'solarization' effect seen when exposing a 
  1012.  * photographic film to light during the development process.
  1013.  * Done here by inverting all pixels above threshold level.
  1014.  */
  1015. Image *
  1016. ImageSolarize(Image * input)
  1017. {
  1018.     Image *output = ImageNewCmap(input->width, input->height, input->cmapSize);
  1019.     int i;
  1020.     unsigned char *ip, *op;
  1021.     int count, limit;
  1022.  
  1023.     limit = ImgProcessInfo.solarizeThreshold * 255 / 100;
  1024.  
  1025.     /*
  1026.     **     If the input has a colormap, just tweak that.
  1027.      */
  1028.     if (input->cmapSize != 0) {
  1029.     ip = input->cmapData;
  1030.     op = output->cmapData;
  1031.     count = input->cmapSize;
  1032.  
  1033.     memcpy(output->data, input->data,
  1034.          sizeof(char) * input->scale * input->width * input->height);
  1035.     } else {
  1036.     ip = input->data;
  1037.     op = output->data;
  1038.     count = input->width * input->height;
  1039.     }
  1040.  
  1041.     for (i = 0; i < count; i++) {
  1042.     *op++ = (*ip > limit) ? 255 - *ip++ : *ip++;
  1043.     *op++ = (*ip > limit) ? 255 - *ip++ : *ip++;
  1044.     *op++ = (*ip > limit) ? 255 - *ip++ : *ip++;
  1045.     }
  1046.  
  1047.     return output;
  1048. }
  1049.  
  1050. /*
  1051.  * Colourmap quantization. Hacked from ppmqvga.c, which is part of Netpbm and
  1052.  * was originally written by Lyle Rains (lrains@netcom.com) and Bill Davidsen
  1053.  * (davidsen@crd.ge.com).
  1054.  * You are probably not supposed to understand this.
  1055.  */
  1056.  
  1057. #define RED_BITS   5
  1058. #define GREEN_BITS 6
  1059. #define BLUE_BITS  5
  1060.  
  1061. #define MAX_RED       (1 << RED_BITS)
  1062. #define MAX_GREEN  (1 << GREEN_BITS)
  1063. #define MAX_BLUE   (1 << BLUE_BITS)
  1064.  
  1065. #define MAXWEIGHT  128
  1066. #define STDWEIGHT_DIV  (2 << 8)
  1067. #define STDWEIGHT_MUL  (2 << 10)
  1068. #define GAIN       4
  1069.  
  1070. static int r, g, b, clutx, rep_threshold, rep_weight, dr, dg, db;
  1071. static int *color_cube, ncolors;
  1072. static unsigned char *clut;
  1073.  
  1074. #define CUBEINDEX(r,g,b)  (r)*(MAX_GREEN*MAX_BLUE) + (g)*MAX_BLUE + (b)
  1075.  
  1076. static void 
  1077. diffuse(void)
  1078. {
  1079.     int _7_32nds, _3_32nds, _1_16th;
  1080.  
  1081.     if (clutx < ncolors) {
  1082.     if (color_cube[CUBEINDEX(r, g, b)] > rep_threshold) {
  1083.         clut[clutx * 4 + 0] = ((2 * r + 1) * 256) / (2 * MAX_RED);
  1084.         clut[clutx * 4 + 1] = ((2 * g + 1) * 256) / (2 * MAX_GREEN);
  1085.         clut[clutx * 4 + 2] = ((2 * b + 1) * 256) / (2 * MAX_BLUE);
  1086.         ++clutx;
  1087.         color_cube[CUBEINDEX(r, g, b)] -= rep_weight;
  1088.     }
  1089.     _7_32nds = (7 * color_cube[CUBEINDEX(r, g, b)]) / 32;
  1090.     _3_32nds = (3 * color_cube[CUBEINDEX(r, g, b)]) / 32;
  1091.     _1_16th = color_cube[CUBEINDEX(r, g, b)] - 3 * (_7_32nds + _3_32nds);
  1092.     color_cube[CUBEINDEX(r, g, b)] = 0;
  1093.     /* spread error evenly in color space. */
  1094.     color_cube[CUBEINDEX(r, g, b + db)] += _7_32nds;
  1095.     color_cube[CUBEINDEX(r, g + dg, b)] += _7_32nds;
  1096.     color_cube[CUBEINDEX(r + dr, g, b)] += _7_32nds;
  1097.     color_cube[CUBEINDEX(r, g + dg, b + db)] += _3_32nds;
  1098.     color_cube[CUBEINDEX(r + dr, g, b + db)] += _3_32nds;
  1099.     color_cube[CUBEINDEX(r + dr, g + dg, b)] += _3_32nds;
  1100.     color_cube[CUBEINDEX(r + dr, g + dg, b + db)] += _1_16th;
  1101.  
  1102.     /*
  1103.      * Conserve the error at edges if possible
  1104.      * (which it is, except the last pixel)
  1105.      */
  1106.     if (color_cube[CUBEINDEX(r, g, b)] != 0) {
  1107.         if (dg != 0)
  1108.         color_cube[CUBEINDEX(r, g + dg, b)] +=
  1109.             color_cube[CUBEINDEX(r, g, b)];
  1110.         else if (dr != 0)
  1111.         color_cube[CUBEINDEX(r + dr, g, b)] +=
  1112.             color_cube[CUBEINDEX(r, g, b)];
  1113.         else if (db != 0)
  1114.         color_cube[CUBEINDEX(r, g, b) + db] +=
  1115.             color_cube[CUBEINDEX(r, g, b)];
  1116.         else
  1117.         fprintf(stderr, "lost error term\n");
  1118.     }
  1119.     }
  1120.     color_cube[CUBEINDEX(r, g, b)] = -1;
  1121. }
  1122.  
  1123. /*
  1124. ** Find representative color nearest to requested color.  Check color cube
  1125. ** for a cached color index.  If not cached, compute nearest and cache result.
  1126.  */
  1127. static int 
  1128. nearest_color(unsigned char *p)
  1129. {
  1130.     register unsigned char *test;
  1131.     register unsigned i;
  1132.     unsigned long min_dist_sqd, dist_sqd;
  1133.     int nearest = 0;
  1134.     int *cache;
  1135.     int r, g, b;
  1136.  
  1137.     r = *p++;
  1138.     g = *p++;
  1139.     b = *p;
  1140.     cache = &(color_cube[CUBEINDEX((r << RED_BITS) / 256,
  1141.                    (g << GREEN_BITS) / 256,
  1142.                    (b << BLUE_BITS) / 256)]);
  1143.     if (*cache >= 0)
  1144.     return *cache;
  1145.     min_dist_sqd = ~0;
  1146.     for (i = 0; i < ncolors; ++i) {
  1147.     test = &clut[i * 4];
  1148.     dist_sqd =
  1149.         3 * (r - test[0]) * (r - test[0]) +
  1150.         4 * (g - test[1]) * (g - test[1]) +
  1151.         2 * (b - test[2]) * (b - test[2]);
  1152.     if (dist_sqd < min_dist_sqd) {
  1153.         nearest = i;
  1154.         min_dist_sqd = dist_sqd;
  1155.     }
  1156.     }
  1157.     return (*cache = nearest);
  1158. }
  1159.  
  1160. Image *
  1161. ImageQuantize(Image * input)
  1162. {
  1163.     int w = input->width, h = input->height;
  1164.     Image *output = ImageNew(w, h);
  1165.     unsigned char *op = output->data;
  1166.     int *weight_convert;
  1167.     int k, x, y, i, j, nearest;
  1168.     int total_weight, cum_weight[MAX_GREEN];
  1169.     int *erropt;
  1170.     int *errP;
  1171.     unsigned char *p, *clutP;
  1172.  
  1173.  
  1174.     ncolors = ImgProcessInfo.quantizeColors;
  1175.  
  1176.     /* Make a color cube, initialized to zeros */
  1177.     color_cube = calloc(MAX_RED * MAX_GREEN * MAX_BLUE, sizeof(int));
  1178.     clut = calloc(ncolors * 4, sizeof(int));
  1179.     erropt = calloc(ncolors * 4, sizeof(int));
  1180.  
  1181.     clutx = 0;
  1182.  
  1183.     /* Count all occurrances of each color */
  1184.     for (y = 0; y < h; ++y)
  1185.     for (x = 0; x < w; ++x) {
  1186.         p = ImagePixel(input, x, y);
  1187.         r = p[0] / (256 / MAX_RED);
  1188.         g = p[1] / (256 / MAX_GREEN);
  1189.         b = p[2] / (256 / MAX_BLUE);
  1190.         ++color_cube[CUBEINDEX(r, g, b)];
  1191.     }
  1192.  
  1193.     /* Initialize logarithmic weighting table */
  1194.     weight_convert = malloc(MAXWEIGHT * sizeof(int));
  1195.     weight_convert[0] = 0;
  1196.     for (i = 1; i < MAXWEIGHT; ++i) {
  1197.     weight_convert[i] = (int) (100.0 * log((double) i));
  1198.     }
  1199.  
  1200.     k = w * h;
  1201.     if ((k /= STDWEIGHT_DIV) == 0)
  1202.     k = 1;
  1203.     total_weight = i = 0;
  1204.     for (g = 0; g < MAX_GREEN; ++g) {
  1205.     for (r = 0; r < MAX_RED; ++r) {
  1206.         for (b = 0; b < MAX_BLUE; ++b) {
  1207.         register int weight;
  1208.         /* Normalize the weights, independent of picture size. */
  1209.         weight = color_cube[CUBEINDEX(r, g, b)] * STDWEIGHT_MUL;
  1210.         weight /= k;
  1211.         if (weight)
  1212.             ++i;
  1213.         if (weight >= MAXWEIGHT)
  1214.             weight = MAXWEIGHT - 1;
  1215.         total_weight += (color_cube[CUBEINDEX(r, g, b)]
  1216.                  = weight_convert[weight]);
  1217.         }
  1218.     }
  1219.     cum_weight[g] = total_weight;
  1220.     }
  1221.     rep_weight = total_weight / ncolors;
  1222.  
  1223.     /* Magic foo-foo dust here.     What IS the correct way to select threshold? */
  1224.     rep_threshold = total_weight * (28 + 110000 / i) / 95000;
  1225.  
  1226.     /*
  1227.      * Do a 3-D error diffusion dither on the data in the color cube
  1228.      * to select the representative colors.  Do the dither back and forth in
  1229.      * such a manner that all the error is conserved (none lost at the edges).
  1230.      */
  1231.  
  1232.     dg = 1;
  1233.     for (g = 0; g < MAX_GREEN; ++g) {
  1234.     dr = 1;
  1235.     for (r = 0; r < MAX_RED; ++r) {
  1236.         db = 1;
  1237.         for (b = 0; b < MAX_BLUE - 1; ++b)
  1238.         diffuse();
  1239.         db = 0;
  1240.         diffuse();
  1241.         ++b;
  1242.         if (++r == MAX_RED - 1)
  1243.         dr = 0;
  1244.         db = -1;
  1245.         while (--b > 0)
  1246.         diffuse();
  1247.         db = 0;
  1248.         diffuse();
  1249.     }
  1250.     /* Modify threshold to keep rep points proportionally distributed */
  1251.     if ((j = clutx - (ncolors * cum_weight[g]) / total_weight) != 0)
  1252.         rep_threshold += j * GAIN;
  1253.  
  1254.     if (++g == MAX_GREEN - 1)
  1255.         dg = 0;
  1256.     dr = -1;
  1257.     while (r-- > 0) {
  1258.         db = 1;
  1259.         for (b = 0; b < MAX_BLUE - 1; ++b)
  1260.         diffuse();
  1261.         db = 0;
  1262.         diffuse();
  1263.         ++b;
  1264.         if (--r == 0)
  1265.         dr = 0;
  1266.         db = -1;
  1267.         while (--b > 0)
  1268.         diffuse();
  1269.         db = 0;
  1270.         diffuse();
  1271.     }
  1272.     /* Modify threshold to keep rep points proportionally distributed */
  1273.     if ((j = clutx - (ncolors * cum_weight[g]) / total_weight) != 0)
  1274.         rep_threshold += j * GAIN;
  1275.     }
  1276.  
  1277.     /*
  1278.      * Check the error associated with the use of each color, and
  1279.      * change the value of the color to minimize the error.
  1280.      */
  1281.     for (y = 0; y < h; ++y) {
  1282.     for (x = 0; x < w; ++x) {
  1283.         p = ImagePixel(input, x, y);
  1284.         nearest = nearest_color(p);
  1285.         errP = &erropt[nearest * 4];
  1286.         clutP = &clut[nearest * 4];
  1287.         errP[0] += *p++ - clutP[0];
  1288.         errP[1] += *p++ - clutP[1];
  1289.         errP[2] += *p - clutP[2];
  1290.         ++errP[3];
  1291.     }
  1292.     }
  1293.  
  1294.     for (i = 0; i < ncolors; ++i) {
  1295.     clutP = &clut[i * 4];
  1296.     errP = &erropt[i * 4];
  1297.     j = errP[3];
  1298.     if (j > 0)
  1299.         j *= 4;
  1300.     else if (j == 0)
  1301.         j = 1;
  1302.     clutP[0] += (errP[0] / j) * 4;
  1303.     clutP[1] += (errP[1] / j) * 4;
  1304.     clutP[2] += (errP[2] / j) * 4;
  1305.     }
  1306.  
  1307.     /* Reset the color cache. */
  1308.     for (i = 0; i < MAX_RED * MAX_GREEN * MAX_BLUE; ++i)
  1309.     color_cube[i] = -1;
  1310.  
  1311.     /*
  1312.      * Map the colors in the image to their closest match in the new colormap.
  1313.      */
  1314.     for (y = 0; y < h; ++y) {
  1315.     for (x = 0; x < w; ++x) {
  1316.         p = ImagePixel(input, x, y);
  1317.         nearest = nearest_color(p);
  1318.         clutP = &clut[nearest * 4];
  1319.         *op++ = *clutP++;
  1320.         *op++ = *clutP++;
  1321.         *op++ = *clutP;
  1322.     }
  1323.     }
  1324.  
  1325.     free(clut);
  1326.     free(erropt);
  1327.     free(weight_convert);
  1328.  
  1329.     return output;
  1330. }
  1331.  
  1332. /*
  1333.  * Convert an image to grey scale.
  1334.  */
  1335. Image *
  1336. ImageGrey(Image * input)
  1337. {
  1338.     Image *output = ImageNew(input->width, input->height);
  1339.     unsigned char *p, *op = output->data;
  1340.     int x, y, v;
  1341.  
  1342.  
  1343.     for (y = 0; y < input->height; y++) {
  1344.     for (x = 0; x < input->width; x++) {
  1345.         p = ImagePixel(input, x, y);
  1346.         v = (p[0] * 11 + p[1] * 16 + p[2] * 5) >> 5;    /* .33R + .5G + .17B */
  1347.         *op++ = v;
  1348.         *op++ = v;
  1349.         *op++ = v;
  1350.     }
  1351.     }
  1352.     return output;
  1353. }
  1354.  
  1355. #define GRAY(p)        (((p)[0]*169 + (p)[1]*256 + (p)[2]*87)/512)
  1356. #define SQ(x)        ((x)*(x))
  1357.  
  1358. /*
  1359.  * Directional filter, according to the algorithm described on p. 60 of:
  1360.  * _Algorithms_for_Graphics_and_Image_Processing_. Theo Pavlidis.
  1361.  * Computer Science Press, 1982.
  1362.  * For each pixel, detect the most prominent edge and apply a filter that
  1363.  * does not degrade that edge.
  1364.  */
  1365. Image *
  1366. ImageDirectionalFilter(Image * input)
  1367. {
  1368.     unsigned char *p00, *p10, *p_0, *p11, *p__, *p01, *p0_, *p_1, *p1_;
  1369.     int g00, g10, g_0, g11, g__, g01, g0_, g_1, g1_;
  1370.     int v0, v45, v90, v135, vmin, theta;
  1371.     int x, y;
  1372.     Image *output = ImageNew(input->width, input->height);
  1373.     unsigned char *op = output->data;
  1374.  
  1375.     /*
  1376.      * We don't process the border of the image tto avoid the hassle
  1377.      * of having do deal with boundary conditions. Hopefully no one
  1378.      * will notice.
  1379.      */
  1380.     /* copy first row unchanged */
  1381.     for (x = 0; x < input->width; x++) {
  1382.     p00 = ImagePixel(input, x, 0);
  1383.     *op++ = *p00++;
  1384.     *op++ = *p00++;
  1385.     *op++ = *p00++;
  1386.     }
  1387.  
  1388.     for (y = 1; y < input->height - 1; y++) {
  1389.     /* copy first column unchanged */
  1390.     p00 = ImagePixel(input, 0, y);
  1391.     *op++ = *p00++;
  1392.     *op++ = *p00++;
  1393.     *op++ = *p00++;
  1394.     for (x = 1; x < input->width - 1; x++) {
  1395.         /* find values of pixel and all neighbours */
  1396.         p00 = ImagePixel(input, x, y);
  1397.         p10 = ImagePixel(input, x + 1, y);
  1398.         p_0 = ImagePixel(input, x - 1, y);
  1399.         p11 = ImagePixel(input, x + 1, y + 1);
  1400.         p__ = ImagePixel(input, x - 1, y - 1);
  1401.         p01 = ImagePixel(input, x, y + 1);
  1402.         p0_ = ImagePixel(input, x, y - 1);
  1403.         p_1 = ImagePixel(input, x - 1, y + 1);
  1404.         p1_ = ImagePixel(input, x + 1, y - 1);
  1405.  
  1406.         /* get grayscale values */
  1407.         g00 = GRAY(p00);
  1408.         g01 = GRAY(p01);
  1409.         g10 = GRAY(p10);
  1410.         g11 = GRAY(p11);
  1411.         g0_ = GRAY(p0_);
  1412.         g_0 = GRAY(p_0);
  1413.         g__ = GRAY(p__);
  1414.         g_1 = GRAY(p_1);
  1415.         g1_ = GRAY(p1_);
  1416.  
  1417.         /* estimate direction of edge, if any */
  1418.         v0 = SQ(g00 - g10) + SQ(g00 - g_0);
  1419.         v45 = SQ(g00 - g11) + SQ(g00 - g__);
  1420.         v90 = SQ(g00 - g01) + SQ(g00 - g0_);
  1421.         v135 = SQ(g00 - g_1) + SQ(g00 - g1_);
  1422.  
  1423.         vmin = MIN(MIN(v0, v45), MIN(v90, v135));
  1424.         theta = 0;
  1425.         if (vmin == v45)
  1426.         theta = 1;
  1427.         else if (vmin == v90)
  1428.         theta = 2;
  1429.         else if (vmin == v135)
  1430.         theta = 3;
  1431.  
  1432.         /* apply filtering according to direction of edge */
  1433.         switch (theta) {
  1434.         case 0:        /* 0 degrees */
  1435.         *op++ = (*p_0++ + *p00++ + *p10++) / 3;
  1436.         *op++ = (*p_0++ + *p00++ + *p10++) / 3;
  1437.         *op++ = (*p_0++ + *p00++ + *p10++) / 3;
  1438.         break;
  1439.         case 1:        /* 45 degrees */
  1440.         *op++ = (*p__++ + *p00++ + *p11++) / 3;
  1441.         *op++ = (*p__++ + *p00++ + *p11++) / 3;
  1442.         *op++ = (*p__++ + *p00++ + *p11++) / 3;
  1443.         break;
  1444.         case 2:        /* 90 degrees */
  1445.         *op++ = (*p0_++ + *p00++ + *p01++) / 3;
  1446.         *op++ = (*p0_++ + *p00++ + *p01++) / 3;
  1447.         *op++ = (*p0_++ + *p00++ + *p01++) / 3;
  1448.         break;
  1449.         case 3:        /* 135 degrees */
  1450.         *op++ = (*p1_++ + *p00++ + *p_1++) / 3;
  1451.         *op++ = (*p1_++ + *p00++ + *p_1++) / 3;
  1452.         *op++ = (*p1_++ + *p00++ + *p_1++) / 3;
  1453.         break;
  1454.         }
  1455.     }
  1456.     /* copy last column unchanged */
  1457.     p00 = ImagePixel(input, x, y);
  1458.     *op++ = *p00++;
  1459.     *op++ = *p00++;
  1460.     *op++ = *p00++;
  1461.     }
  1462.  
  1463.     /* copy last row unchanged */
  1464.     for (x = 0; x < input->width; x++) {
  1465.     p00 = ImagePixel(input, x, y);
  1466.     *op++ = *p00++;
  1467.     *op++ = *p00++;
  1468.     *op++ = *p00++;
  1469.     }
  1470.     return output;
  1471. }
  1472.  
  1473. #if 0
  1474. #define ISFOREGND(p) ((p) ? \
  1475.               (deltaR-deltaRV <= p[0]) && (p[0] <= deltaR+deltaRV) && \
  1476.               (deltaG-deltaGV <= p[1]) && (p[1] <= deltaG+deltaGV) && \
  1477.               (deltaB-deltaBV <= p[2]) && (p[2] <= deltaB+deltaBV) : 1)
  1478.  
  1479. /*
  1480.  * Thicken an image.
  1481.  */
  1482. Image *
  1483. ImageThicken(Image * input)
  1484. {
  1485.     unsigned char *p00, *p10, *p_0, *p01, *p0_;
  1486.     int x, y, width, height, br, bg, bb;
  1487.     int deltaR, deltaRV, deltaG, deltaGV, deltaB, deltaBV;
  1488.     Image *output;
  1489.     unsigned char *op;
  1490.  
  1491.     width = input->width;
  1492.     height = input->height;
  1493.     output = ImageNew(width, height);
  1494.     op = output->data;
  1495.  
  1496.     /* Get RGB values of background colour */
  1497.     br = ImgProcessInfo.background->red / 256;
  1498.     bg = ImgProcessInfo.background->green / 256;
  1499.     bb = ImgProcessInfo.background->blue / 256;
  1500.  
  1501.     for (y = 0; y < height; y++)
  1502.     for (x = 0; x < width; x++) {
  1503.         /* find values of pixel and all d-neighbours */
  1504.         p00 = ImagePixel(input, x, y);
  1505.         p10 = x < width - 1 ? ImagePixel(input, x + 1, y) : NULL;
  1506.         p_0 = x > 0 ? ImagePixel(input, x - 1, y) : NULL;
  1507.         p01 = y < height - 1 ? ImagePixel(input, x, y + 1) : NULL;
  1508.         p0_ = y > 0 ? ImagePixel(input, x, y - 1) : NULL;
  1509.  
  1510.         if (ISFOREGND(p00)) {
  1511.         *op++ = *p00++;
  1512.         *op++ = *p00++;
  1513.         *op++ = *p00++;
  1514.         } else {
  1515.         *op++ = br;
  1516.         *op++ = bg;
  1517.         *op++ = bb;
  1518.         }
  1519.     }
  1520.  
  1521.     return output;
  1522. }
  1523. #endif
  1524.