home *** CD-ROM | disk | FTP | other *** search
/ rtsi.com / 2014.01.www.rtsi.com.tar / www.rtsi.com / OS9 / OSK / NETWORK / netpbm_src.lzh / NETPBM / PNM / pnmhisteq.c < prev    next >
C/C++ Source or Header  |  1996-11-18  |  9KB  |  372 lines
  1. /*
  2.                  pnmhisteq.c
  3.  
  4.            Equalise histogram for a PGM or PPM file
  5.  
  6.    Options:    -gray:      modify gray pixels only; leave colours unchanged
  7.         -rmap fn: read luminosity map from PGM file fn
  8.         -wmap wn: write luminosity map to PGM file fn
  9.         -verbose: print histogram and luminosity map
  10.  
  11.     Accepts PGM and PPM as input.  PBM input is allowed, but histogram
  12.     equalisation does not modify a PBM file.
  13.  
  14.       by John Walker (kelvin@fourmilab.ch) -- March MVM.
  15.            WWW home page: http://www.fourmilab.ch/
  16.  
  17.           Copyright (C) 1995 by John Walker
  18.  
  19.     Permission    to use, copy, modify, and distribute this software and
  20.     its documentation for  any    purpose  and  without  fee  is    hereby
  21.     granted, without any conditions or restrictions.  This software is
  22.     provided "as is" without express or implied warranty.
  23.  
  24. */
  25.  
  26. #include "pnm.h"
  27.  
  28. /* Prototypes */
  29.  
  30. static void hsv_rgb ARGS((double h, double s, double v,
  31.               double *r, double *g, double *b));
  32. static void rgb_hsv ARGS((double r, double g, double b,
  33.               double *h, double *s, double *v));
  34.  
  35. /*  HSV_RGB  --  Convert HSV colour specification to RGB  intensities.
  36.          Hue is specified as a    real  value  from  0  to  360,
  37.          Saturation  and  Intensity as reals from 0 to 1.  The
  38.          RGB components are returned as reals from 0 to 1. */
  39.  
  40. static void hsv_rgb(h, s, v, r, g, b)
  41.   double h, s, v;
  42.   double *r, *g, *b;
  43. {
  44.     int i;
  45.     double f, p, q, t;
  46.  
  47.     if (s == 0) {
  48.     *r = *g = *b = v;
  49.     } else {
  50.     if (h == 360.0) {
  51.         h = 0;
  52.     }
  53.     h /= 60.0;
  54.  
  55.     i = h;
  56.     f = h - i;
  57.     p = v * (1.0 - s);
  58.     q = v * (1.0 - (s * f));
  59.     t = v * (1.0 - (s * (1.0 - f)));
  60.     switch (i) {
  61.  
  62.         case 0:
  63.         *r = v;
  64.         *g = t;
  65.         *b = p;
  66.         break;
  67.  
  68.         case 1:
  69.         *r = q;
  70.         *g = v;
  71.         *b = p;
  72.         break;
  73.  
  74.         case 2:
  75.         *r = p;
  76.         *g = v;
  77.         *b = t;
  78.         break;
  79.  
  80.         case 3:
  81.         *r = p;
  82.         *g = q;
  83.         *b = v;
  84.         break;
  85.  
  86.         case 4:
  87.         *r = t;
  88.         *g = p;
  89.         *b = v;
  90.         break;
  91.  
  92.         case 5:
  93.         *r = v;
  94.         *g = p;
  95.         *b = q;
  96.         break;
  97.      }
  98.     }
  99. }
  100.  
  101. /*  RGB_HSV  --  Map R, G, B intensities in the range from 0 to 1 into
  102.          Hue, Saturation,  and    Value:    Hue  from  0  to  360,
  103.          Saturation  from  0  to  1,  and  Value  from 0 to 1.
  104.                  Special case: if Saturation is 0 (it's a  grey  scale
  105.          tone), Hue is undefined and is returned as -1.
  106.  
  107.          This follows Foley & van Dam, section 17.4.4. */
  108.  
  109. static void rgb_hsv(r, g, b, h, s, v)
  110.   double r, g, b;
  111.   double *h, *s, *v;
  112. {
  113.     double imax = max(r, max(g, b)),
  114.        imin = min(r, min(g, b)),
  115.        rc, gc, bc;
  116.  
  117.     *v = imax;
  118.     if (imax != 0) {
  119.     *s = (imax - imin) / imax;
  120.     } else {
  121.     *s = 0;
  122.     }
  123.  
  124.     if (*s == 0) {
  125.     *h = -1;
  126.     } else {
  127.     rc = (imax - r) / (imax - imin);
  128.     gc = (imax - g) / (imax - imin);
  129.     bc = (imax - b) / (imax - imin);
  130.     if (r == imax) {
  131.         *h = bc - gc;
  132.     } else if (g == imax) {
  133.         *h = 2.0 + rc - bc;
  134.     } else {
  135.         *h = 4.0 + gc - rc;
  136.     }
  137.     *h *= 60.0;
  138.     if (*h < 0.0) {
  139.         *h += 360.0;
  140.     }
  141.     }
  142. }
  143.  
  144. int main(argc, argv)
  145.   int argc;
  146.   char *argv[];
  147. {
  148.     FILE *ifp;
  149.     int argn = 1, i, j, verbose = 0, mono_only = 0;
  150.     gray lmin, lmax;
  151.     gray **lumamap;              /* Luminosity map */
  152.     long *lumahist;              /* Histogram of luminosity values */
  153.     int rows, hist_cols;          /* Rows, columns of input image */
  154.     xelval maxval;              /* Maxval of input image */
  155.     int format;               /* Format indicator (PBM/PGM/PPM) */
  156.     xel** xels;               /* Pixel array */
  157.     unsigned long pixels = 0, pixsum = 0, maxluma = 0;
  158.     double lscale;
  159.     xel *grayrow;
  160.     pixel *pixrow;
  161.     FILE *rmap = NULL, *wmap = NULL;
  162.     char *usage = "[-gray] [-verbose] [-rmap pgmfile] [-wmap pgmfile] [pnmfile]";
  163.  
  164.     pnm_init(&argc, argv);
  165.  
  166.     /* Check for flags. */
  167.  
  168.     while (argn < argc && argv[argn][0] == '-' && argv[argn][1] != '\0') {
  169.         if (pm_keymatch(argv[argn], "-gray", 1)) {
  170.         mono_only = 1;
  171.         } else if (pm_keymatch(argv[argn], "-verbose", 1)) {
  172.         verbose = 1;
  173.         } else if (pm_keymatch(argv[argn], "-rmap", 1)) {
  174.         if (rmap != NULL) {
  175.                 pm_error("already specified an input map");
  176.         }
  177.         argn++;
  178.             if (argn == argc || strcmp(argv[argn], "-") == 0) {
  179.         pm_usage(usage);
  180.         }
  181.         rmap = pm_openr(argv[argn]);
  182.         } else if (pm_keymatch(argv[argn], "-wmap", 1)) {
  183.         if (wmap != NULL) {
  184.                 pm_error("already specified an output map");
  185.         }
  186.         argn++;
  187.         if (argn == argc) {
  188.         pm_usage(usage);
  189.         }
  190.         wmap = pm_openw(argv[argn]);
  191.     } else {
  192.         pm_usage(usage);
  193.     }
  194.     argn++;
  195.     }
  196.  
  197.     if (--argc > argn) {
  198.     pm_usage(usage);
  199.     } else if (argc == argn) {
  200.     ifp = pm_openr(argv[argn]);
  201.     } else {
  202.     ifp = stdin;
  203.     }
  204.  
  205.     xels = pnm_readpnm(ifp, &hist_cols, &rows, &maxval, &format);
  206.     pm_close(ifp);
  207.  
  208.     /* Allocate histogram and luminosity map arrays.  If the
  209.        user has specified an input map file, read it in at
  210.        this point. */
  211.  
  212.     lumahist = (long *) pm_allocrow(maxval + 1, sizeof(long));
  213.     bzero((char *) lumahist, (maxval + 1) * sizeof(long));
  214.  
  215.     if (rmap == NULL) {
  216.     lumamap = pgm_allocarray(maxval + 1, 1);
  217.     } else {
  218.     int rmcols, rmrows; 
  219.     gray rmmaxv;
  220.  
  221.     lumamap = pgm_readpgm(rmap, &rmcols, &rmrows, &rmmaxv);
  222.     if (rmmaxv != maxval) {
  223.             pm_error("maxval in map file (%d) different from input (%d)",
  224.         rmmaxv, maxval);
  225.     }
  226.     if (rmrows != 1 || rmcols != rmmaxv) {
  227.             pm_error("map size (%d by %d) wrong; must be (%d by 1)",
  228.         rmcols, rmrows, maxval);
  229.     }
  230.     }
  231.  
  232.     /* Scan the image and build the luminosity histogram.  If
  233.        the input is a PPM, we calculate the luminosity of each
  234.        pixel from its RGB components. */
  235.  
  236.     lmin = maxval;
  237.     lmax = 0;
  238.     if (PNM_FORMAT_TYPE(format) == PGM_TYPE ||
  239.     PNM_FORMAT_TYPE(format) == PBM_TYPE) {
  240.  
  241.     /* Compute intensity histogram */
  242.  
  243.     pixels = ((unsigned long) rows) * ((unsigned long) hist_cols);
  244.     for (i = 0; i < rows; i++) {
  245.         xel *grayrow = xels[i];
  246.         for (j = 0; j < hist_cols; j++) {
  247.         gray l = PNM_GET1(grayrow[j]);
  248.         lmin = min(lmin, l);
  249.         lmax = max(lmax, l);
  250.         lumahist[l]++;
  251.         }
  252.     }
  253.     } else if (PNM_FORMAT_TYPE(format) == PPM_TYPE) {
  254.     for (i = 0; i < rows; i++) {
  255.         pixel *pixrow = (pixel *) xels[i];
  256.  
  257.         for (j = 0; j < hist_cols; j++) {
  258.         if (!mono_only ||
  259.             ((PPM_GETR(pixrow[j]) == PPM_GETG(pixrow[j])) &&
  260.              (PPM_GETR(pixrow[j]) == PPM_GETB(pixrow[j])))) {
  261.             gray l = (gray) PPM_LUMIN(pixrow[j]);
  262.             lmin = min(lmin, l);
  263.             lmax = max(lmax, l);
  264.             lumahist[l]++;
  265.             pixels++;
  266.         }
  267.         }
  268.     }
  269.     } else {
  270.         pm_error("unknown input format");
  271.     }
  272.  
  273.     /* The PGM and PPM branches rejoin here to calculate the
  274.        luminosity mapping table which gives the histogram-equalised
  275.        luminosity for each original luminosity. */
  276.  
  277.     /* Calculate initial histogram equalisation curve. */
  278.  
  279.     for (i = 0; i <= (int) maxval; i++) {
  280.  
  281.     /* Yick.  If PGM_BIGGRAYS is defined (I thought they were little
  282.        guys, about four foot, with funny eyes...) the following
  283.        calculation can overflow a 32 bit long.  So, we do it in
  284.        floating point.  Since this happens only maxval times, the
  285.        inefficiency is trivial compared to the every-pixel code above
  286.        and below. */
  287.  
  288.     lumamap[0][i] = (gray) (((((double) pixsum * maxval)) / pixels) + 0.5);
  289.     if (lumahist[i] > 0) {
  290.         maxluma = i;
  291.     }
  292.     pixsum += lumahist[i];
  293.     }
  294.  
  295.     /* Normalise so that the brightest pixels are set to
  296.        maxval. */
  297.  
  298.     lscale = ((double) maxval) / ((lumahist[maxluma] > 0) ?
  299.          ((double) lumamap[0][maxluma]) : ((double) maxval));
  300.     for (i = 0; i <= (int) maxval; i++) {
  301.     lumamap[0][i] = (gray)
  302.         min(((long) maxval), ((long) (lumamap[0][i] * lscale + 0.5)));
  303.     }
  304.  
  305.     /* If requested, print the luminosity map and original histogram. */
  306.  
  307.     if (verbose) {
  308.     fprintf(stderr,
  309.             "  Luminosity map    Number of\n Original    New     Pixels\n");
  310.     for (i = 0; i <= (int) maxval; i++) {
  311.         if (lumahist[i] > 0) {
  312.                 fprintf(stderr,"%6d -> %6d  %8d\n", i,
  313.             lumamap[0][i], lumahist[i]);
  314.         }
  315.     }
  316.     }
  317.  
  318.     switch (PNM_FORMAT_TYPE(format)) {
  319.     case PBM_TYPE:
  320.     case PPM_TYPE:
  321.         for (i = 0; i < rows; i++) {
  322.         pixrow = (pixel *) xels[i];
  323.         for (j = 0; j < hist_cols; j++) {
  324.             if (!mono_only ||
  325.             ((PPM_GETR(pixrow[j]) == PPM_GETG(pixrow[j])) &&
  326.              (PPM_GETR(pixrow[j]) == PPM_GETB(pixrow[j])))) {
  327.             double r, g, b, h, s, v;
  328.             int iv;
  329.  
  330.             r = (double) PPM_GETR(pixrow[j]) / ((double) maxval);
  331.             g = (double) PPM_GETG(pixrow[j]) / ((double) maxval);
  332.             b = (double) PPM_GETB(pixrow[j]) / ((double) maxval);
  333.             rgb_hsv(r, g, b, &h, &s, &v);
  334.             iv = (int) ((v * maxval) + 0.5);
  335.  
  336.             if (iv > ((int) maxval)) {
  337.                 iv = maxval;
  338.             }
  339.             v = ((double) lumamap[0][iv]) / ((double) maxval);
  340.             if (v > 1.0) {
  341.                 v = 1.0;
  342.             }
  343.             hsv_rgb(h, s, v, &r, &g, &b);
  344.             PPM_ASSIGN(pixrow[j], (int) (r * maxval),
  345.                 (int) (g * maxval), (int) (b * maxval));
  346.             }
  347.         }
  348.         }
  349.         break;
  350.  
  351.     case PGM_TYPE:
  352.         for (i = 0; i < rows; i++) {
  353.         grayrow = xels[i];
  354.         for (j = 0; j < hist_cols; j++) {
  355.             PNM_ASSIGN1(grayrow[j], lumamap[0][PNM_GET1(grayrow[j])]);
  356.         }
  357.         }
  358.         break;
  359.     }
  360.  
  361.     pnm_writepnm(stdout, xels, hist_cols, rows, maxval, format, 0);
  362.  
  363.     /* If requested, save the map as a PGM file. */
  364.  
  365.     if (wmap != NULL) {
  366.     pgm_writepgm(wmap, lumamap, maxval, 1, maxval, 0);
  367.     fclose(wmap);
  368.     }
  369.  
  370.     return 0;
  371. }
  372.