home *** CD-ROM | disk | FTP | other *** search
/ rtsi.com / 2014.01.www.rtsi.com.tar / www.rtsi.com / OS9 / OSK / GRAPHICS / mgif.lzh / GIF / process.c < prev    next >
Text File  |  1992-01-13  |  23KB  |  1,429 lines
  1. #undef DBL_LOOP            /* def for for(x..., for(y... loops */
  2.  
  3. /*
  4.  *    process.c - all the image processing functions.
  5.  */
  6.  
  7. static char *sccsid = "@(#) process.c 1.1 91/6/3 rosenkra\0\0                 ";
  8.  
  9. #include <stdio.h>
  10. #include "mgif.h"
  11.  
  12. /*
  13.  *    local functions
  14.  */
  15. int        larger ();        /* frame processes... */
  16. int        smaller ();
  17. int        convolve ();
  18. int        blur ();        /* area processes... */
  19. int        median ();
  20. int        piksrt ();
  21. int        logscale ();        /* point processes... */
  22. int        Log2 ();
  23. int        contrast ();
  24. int        brighten ();
  25. int        invert ();
  26. int        threshold ();
  27. int        histeq ();
  28. int        redistribute ();
  29. int        copyrast ();        /* others... */
  30.  
  31. /*
  32.  *    external functions
  33.  */
  34. extern int    do_hist ();
  35. extern int    copyrast ();
  36.  
  37.  
  38. /*------------------------------*/
  39. /*    larger            */
  40. /*------------------------------*/
  41. int larger (pras, w, h, ptrans)
  42. uchar_t           *pras;
  43. int        w;
  44. int        h;
  45. uchar_t           *ptrans;
  46. {
  47.  
  48. /*
  49.  *    enlarge an image (so far 320x200 -> 640x400). return 1 if
  50.  *    successful, else 0.
  51.  */
  52.  
  53.     uchar_t               *pr;
  54.     uchar_t               *pt;
  55.     long            x;
  56.     long            y;
  57.     long            x2;
  58.     long            y2;
  59.     long            wold;
  60.     long            hold;
  61.     long            w2;
  62.     long            h2;
  63.  
  64.  
  65.     if ((w > 320) || (h > 200))
  66.         return (0);
  67.  
  68.  
  69.     if (pras == ptrans)
  70.     {
  71.         /*
  72.          *   do in place...
  73.          */
  74.         pr   = pras;
  75.         wold = w;
  76.         hold = h;
  77.         w2   = wold*2L;
  78.         h2   = hold*2L;
  79.  
  80.         /*
  81.          *   first shift existing data. this opens up pixels in
  82.          *   a row...
  83.          */
  84.         for (y = hold; y > 0; y--)
  85.         {
  86.             for (x = wold; x > 0; x--)
  87.             {
  88.                 x2 = (2*x)-1;
  89.                 y2 = (2*y)-1;
  90.                 pr[y2*w2 + x2] = pr[(y-1)*wold + (x-1)];
  91.             }
  92.         }
  93.  
  94.         /*
  95.          *   now make intermediate points within (new) row
  96.          */
  97.         for (y2 = 1; y2 < h2; y2 += 2)
  98.         {
  99.             /* first point is duplicate of second */
  100.             pr[y2*w2] = pr[y2*w2 + 1];
  101.  
  102.             /* others are average */
  103.             for (x2 = 2; x2 < w2; x2 += 2)
  104.             {
  105.                 pr[y2*w2 + x2] = (uchar_t) (((uint_t) pr[y2*w2 + (x2-1)]
  106.                            +             (uint_t) pr[y2*w2 + (x2+1)]) / 2);
  107.             }
  108.         }
  109.  
  110.         /*
  111.          *   first row is duplicate of second row
  112.          */
  113.         for (x2 = 0; x2 < w2; x2++)
  114.         {
  115.             pr[x2] = pr[w2 + x2];
  116.         }
  117.  
  118.         /*
  119.          *   finally, make new rows
  120.          */
  121.         for (y2 = 2; y2 < h2; y2 += 2)
  122.         {
  123.             for (x2 = 0; x2 < w2; x2++)
  124.             {
  125.                 pr[y2*w2 + x2] = (uchar_t) (((uint_t) pr[(y2-1)*w2 + x2]
  126.                            +             (uint_t) pr[(y2+1)*w2 + x2]) / 2);
  127.             }
  128.         }
  129.     }
  130.     else
  131.     {
  132.         /*
  133.          *   copy to new image
  134.          */
  135.         pr   = pras;
  136.         pt   = ptrans;
  137.         wold = w;
  138.         hold = h;
  139.         w2   = wold*2L;
  140.         h2   = hold*2L;
  141.  
  142.  
  143.         /*
  144.          *   first copy orig image to transform image
  145.          */
  146.         copyrast (pr, wold, hold, pt);
  147.  
  148.  
  149.         /*
  150.          *   now shift existing data...
  151.          */
  152.         for (y = hold; y > 0; y--)
  153.         {
  154.             for (x = wold; x > 0; x--)
  155.             {
  156.                 x2 = (2*x)-1;
  157.                 y2 = (2*y)-1;
  158.                 pt[y2*w2 + x2] = pt[(y-1)*wold + (x-1)];
  159.             }
  160.         }
  161.  
  162.         /*
  163.          *   now make intermediate points within (new) row
  164.          */
  165.         for (y2 = 1; y2 < h2; y2 += 2)
  166.         {
  167.             /* first point is duplicate of second */
  168.             pt[y2*w2] = pt[y2*w2 + 1];
  169.  
  170.             /* others are average */
  171.             for (x2 = 2; x2 < w2; x2 += 2)
  172.             {
  173.                 pt[y2*w2 + x2] = (uchar_t) (((uint_t) pt[y2*w2 + (x2-1)]
  174.                            +             (uint_t) pt[y2*w2 + (x2+1)]) / 2);
  175.             }
  176.         }
  177.  
  178.         /*
  179.          *   first row is duplicate of second row
  180.          */
  181.         for (x2 = 0; x2 < w2; x2++)
  182.         {
  183.             pt[x2] = pt[w2 + x2];
  184.         }
  185.  
  186.         /*
  187.          *   finally, make new rows
  188.          */
  189.         for (y2 = 2; y2 < h2; y2 += 2)
  190.         {
  191.             for (x2 = 0; x2 < w2; x2++)
  192.             {
  193.                 pt[y2*w2 + x2] = (uchar_t) (((uint_t) pt[(y2-1)*w2 + x2]
  194.                            +             (uint_t) pt[(y2+1)*w2 + x2]) / 2);
  195.             }
  196.         }
  197.     }
  198.  
  199.     return (1);
  200. }
  201.  
  202.  
  203.  
  204.  
  205. /*------------------------------*/
  206. /*    smaller            */
  207. /*------------------------------*/
  208. int smaller (pras, w, h, ptrans)
  209. uchar_t           *pras;
  210. int        w;
  211. int        h;
  212. uchar_t           *ptrans;
  213. {
  214.  
  215. /*
  216.  *    reduce an image by half. return 1 if successful, else 0.
  217.  */
  218.  
  219.     uchar_t               *pr;
  220.     uchar_t               *pt;
  221.     long            x;
  222.     long            y;
  223.     long            x2;
  224.     long            y2;
  225.     long            wold;
  226.     long            hold;
  227.     long            w2;
  228.     long            h2;
  229.  
  230.  
  231.     if ((w < 50) || (h < 50))
  232.         return (0);
  233.  
  234.  
  235.     if (pras == ptrans)
  236.     {
  237.         /*
  238.          *   do in place...
  239.          */
  240.         pr   = pras;
  241.         wold = w;
  242.         hold = h;
  243.         w2   = wold/2L;
  244.         h2   = hold/2L;
  245.  
  246.         /*
  247.          *   first shift existing data. this skips pixels in
  248.          *   a row...
  249.          */
  250.         for (y = 0; y < h2; y++)
  251.         {
  252.             for (x = 0; x < w2; x++)
  253.             {
  254.                 pr[y*w2 + x] = pr[2*y*wold + x*2];
  255.             }
  256.         }
  257.     }
  258.     else
  259.     {
  260.         /*
  261.          *   copy to new image
  262.          */
  263.         pr   = pras;
  264.         pt   = ptrans;
  265.         wold = w;
  266.         hold = h;
  267.         w2   = wold/2L;
  268.         h2   = hold/2L;
  269.  
  270.         /*
  271.          *   do it. skip odd lines and pixels
  272.          */
  273.         for (y = 0; y < h2; y++)
  274.         {
  275.             for (x = 0; x < w2; x++)
  276.             {
  277.                 pt[y*w2 + x] = pr[2*y*wold + x*2];
  278.             }
  279.         }
  280.     }
  281.  
  282.     return (1);
  283. }
  284.  
  285.  
  286.  
  287.  
  288. /*------------------------------*/
  289. /*    convolve        */
  290. /*------------------------------*/
  291.  
  292. /*
  293.  *    built-in convolution kernels
  294.  *
  295.  *    k[0..8] is kernel (always 3x3)
  296.  *    k[9] is scaling flag: -1=divide, 1=multiply, 0=no scaling
  297.  *    k[10] is scale
  298.  *
  299.  *    example:
  300.  *        for lp1 with flag -1 (divide) and scale 9, use
  301.  *
  302.  *            p = ((p0 * k0) +...+ (p8 * k8)) / 9
  303.  *
  304.  *        for hp2 with flag 0 (no scaling), use
  305.  *
  306.  *            p = (p0 * k0) +...+ (p8 * k8)
  307.  */
  308. int            lp1[11] = { 1, 1, 1,        /* low pass */
  309.                     1, 1, 1,
  310.                     1, 1, 1,    -1, 9};
  311. int            lp2[11] = { 1, 1, 1,
  312.                     1, 2, 1,
  313.                     1, 1, 1,    -1, 10};
  314. int            lp3[11] = { 1, 2, 1,
  315.                     2, 4, 2,
  316.                     1, 2, 1,    -1, 16};
  317.  
  318. int            hp1[11] = {-1,-1,-1,        /* hi pass */
  319.                    -1, 9,-1,
  320.                    -1,-1,-1,    0, 1};
  321. int            hp2[11] = { 0,-1, 0,
  322.                    -1, 5,-1,
  323.                     0,-1, 0,    0, 1};
  324. int            hp3[11] = { 1,-2, 1,
  325.                    -2, 5,-2,
  326.                     1,-2, 1,    0, 1};
  327.  
  328.                             /* shift edge */
  329. int            se1[11] = { 0, 0, 0,        /* vertical */
  330.                    -1, 1, 0,
  331.                     0, 0, 0,    0, 1};
  332. int            se2[11] = { 0,-1, 0,        /* horizontal */
  333.                     0, 1, 0,
  334.                     0, 0, 0,    0, 1};
  335. int            se3[11] = {-1, 0, 0,        /* hor and vert */
  336.                     0, 1, 0,
  337.                     0, 0, 0,    0, 1};
  338.  
  339.                             /* gradient edge */
  340. int            ge1[11] = { 1, 1, 1,        /* north */
  341.                     1,-2, 1,
  342.                    -1,-1,-1,    0, 1};
  343. int            ge2[11] = { 1, 1, 1,        /* northeast */
  344.                    -1,-2, 1,
  345.                    -1,-1, 1,    0, 1};
  346. int            ge3[11] = {-1, 1, 1,        /* east */
  347.                    -1,-2, 1,
  348.                    -1, 1, 1,    0, 1};
  349. int            ge4[11] = {-1,-1, 1,        /* southeast */
  350.                    -1,-2, 1,
  351.                     1, 1, 1,    0, 1};
  352. int            ge5[11] = {-1,-1,-1,        /* south */
  353.                     1,-2, 1,
  354.                     1, 1, 1,    0, 1};
  355. int            ge6[11] = { 1,-1,-1,        /* southwest */
  356.                     1,-2,-1,
  357.                     1, 1, 1,    0, 1};
  358. int            ge7[11] = { 1, 1,-1,        /* west */
  359.                     1,-2,-1,
  360.                     1, 1,-1,    0, 1};
  361. int            ge8[11] = { 1, 1, 1,        /* northwest */
  362.                     1,-2,-1,
  363.                     1,-1,-1,    0, 1};
  364.  
  365. int            le1[11] = { 0, 1, 0,        /* laplace edge */
  366.                     1,-4, 1,
  367.                     0, 1, 0,    0, 1};
  368. int            le2[11] = {-1,-1,-1,
  369.                    -1, 8,-1,
  370.                    -1,-1,-1,    0, 1};
  371. int            le3[11] = {-1,-1,-1,
  372.                    -1, 9,-1,
  373.                    -1,-1,-1,    0, 1};
  374. int            le4[11] = { 1,-2, 1,
  375.                    -2, 4,-2,
  376.                     1,-2, 1,    0, 1};
  377.  
  378. int convolve (pras, w, h, opt, kernel, ptrans)
  379. uchar_t           *pras;
  380. register int    w;
  381. int        h;
  382. int        opt;        /* 0=user, 1,2,...=built-in */
  383. int           *kernel;        /* user kernel, if any */
  384. uchar_t           *ptrans;
  385. {
  386.  
  387. /*
  388.  *    convolute an image (3x3). return 1 if successful, else 0.
  389.  *
  390.  *    the basic algorithm is:
  391.  *
  392.  *        new[i] = sum (k  * p ) for i pixels in 3x3 neighborhood
  393.  *                       i    i
  394.  */
  395.  
  396.     uchar_t               *pr;
  397.     uchar_t               *pt;
  398.     register int        val;
  399.     register long        x;
  400.     register long        y;
  401.     register uchar_t       *ps;
  402.     register int           *pk;
  403.  
  404.  
  405.     pr = pras;
  406.     pt = ptrans;
  407.  
  408.     if (pras != ptrans)
  409.     {
  410.         /*
  411.          *   first copy orig image to transform image
  412.          */
  413.         copyrast (pr, w, h, pt);
  414.     }
  415.  
  416.  
  417.     /*
  418.      *   set pointer to user or built-in kernels, depending on opt
  419.      */
  420.     switch (opt)
  421.     {
  422.     case USER_KERN:        pk = kernel;        break;
  423.  
  424.     case LP1_KERN:        pk = lp1;        break;
  425.     case LP2_KERN:        pk = lp2;        break;
  426.     case LP3_KERN:        pk = lp3;        break;
  427.  
  428.     case HP1_KERN:        pk = hp1;        break;
  429.     case HP2_KERN:        pk = hp2;        break;
  430.     case HP3_KERN:        pk = hp3;        break;
  431.  
  432.     case SE1_KERN:        pk = se1;        break;
  433.     case SE2_KERN:        pk = se2;        break;
  434.     case SE3_KERN:        pk = se3;        break;
  435.  
  436.     case GE1_KERN:        pk = ge1;        break;
  437.     case GE2_KERN:        pk = ge2;        break;
  438.     case GE3_KERN:        pk = ge3;        break;
  439.     case GE4_KERN:        pk = ge4;        break;
  440.     case GE5_KERN:        pk = ge5;        break;
  441.     case GE6_KERN:        pk = ge6;        break;
  442.     case GE7_KERN:        pk = ge7;        break;
  443.     case GE8_KERN:        pk = ge8;        break;
  444.  
  445.     case LE1_KERN:        pk = le1;        break;
  446.     case LE2_KERN:        pk = le2;        break;
  447.     case LE3_KERN:        pk = le3;        break;
  448.     case LE4_KERN:        pk = le4;        break;
  449.  
  450.     default:        return (0);
  451.     }
  452.  
  453.  
  454.     /*
  455.      *   do it...
  456.      */
  457.     ps = pt;
  458.     for (y = 1; y < h-1; y++)
  459.     {
  460.         for (x = 1; x < w-1; x++)
  461.         {
  462.             val  = ((uint_t) ps[(y-1)*w+(x-1)]) * pk[0];
  463.             val += ((uint_t) ps[(y-1)*w+(x  )]) * pk[1];
  464.             val += ((uint_t) ps[(y-1)*w+(x+1)]) * pk[2];
  465.             val += ((uint_t) ps[(y  )*w+(x-1)]) * pk[3];
  466.             val += ((uint_t) ps[(y  )*w+(x  )]) * pk[4];
  467.             val += ((uint_t) ps[(y  )*w+(x+1)]) * pk[5];
  468.             val += ((uint_t) ps[(y+1)*w+(x-1)]) * pk[6];
  469.             val += ((uint_t) ps[(y+1)*w+(x  )]) * pk[7];
  470.             val += ((uint_t) ps[(y+1)*w+(x+1)]) * pk[8];
  471.  
  472.             if ((pk[9] < 0) && pk[10])
  473.                 val /= pk[10];
  474.             else if (pk[9] > 0)
  475.                 val *= pk[10];
  476.  
  477.             if (val < 0)
  478.                 val = 0;
  479.             else if (val > 255)
  480.                 val = 255;
  481.  
  482.             ps[(y*w) + x] = (uchar_t) val;
  483.         }
  484.     }
  485.  
  486.     return (1);
  487. }
  488.  
  489.  
  490.  
  491.  
  492. /*------------------------------*/
  493. /*    blur            */
  494. /*------------------------------*/
  495. int blur (pras, w, h, ptrans)
  496. uchar_t           *pras;
  497. register int    w;
  498. int        h;
  499. uchar_t           *ptrans;
  500. {
  501.  
  502. /*
  503.  *    blur an image. return 1 if successful, else 0.
  504.  *
  505.  *    the basic algorithm is:
  506.  *
  507.  *        new[i] = average of neighboring cells
  508.  */
  509.  
  510.     uchar_t               *pr;
  511.     uchar_t               *pt;
  512.     register long        x;
  513.     register long        y;
  514.     register uchar_t       *ps;
  515.  
  516.  
  517.  
  518.     pr = pras;
  519.     pt = ptrans;
  520.  
  521.     if (pras != ptrans)
  522.     {
  523.         /*
  524.          *   first copy orig image to transform image
  525.          */
  526.         copyrast (pr, w, h, pt);
  527.     }
  528.  
  529.  
  530.     /*
  531.      *   do it...
  532.      */
  533.     ps = pt;
  534.     for (y = 1; y < h-1; y++)
  535.     {
  536.         for (x = 1; x < w-1; x++)
  537.         {
  538.             ps[(y*w) + x] = (uchar_t) ((
  539.                     (uint_t) ps[(y-1)*w+(x-1)] +
  540.                     (uint_t) ps[(y-1)*w+(x  )] +
  541.                     (uint_t) ps[(y-1)*w+(x+1)] +
  542.                     (uint_t) ps[(y  )*w+(x-1)] +
  543.                     (uint_t) ps[(y  )*w+(x  )] +
  544.                     (uint_t) ps[(y  )*w+(x+1)] +
  545.                     (uint_t) ps[(y+1)*w+(x-1)] +
  546.                     (uint_t) ps[(y+1)*w+(x  )] +
  547.                     (uint_t) ps[(y+1)*w+(x+1)]) / 9);
  548.         }
  549.     }
  550.  
  551.     return (1);
  552. }
  553.  
  554.  
  555.  
  556.  
  557. /*------------------------------*/
  558. /*    median            */
  559. /*------------------------------*/
  560. int median (pras, w, h, ptrans)
  561. uchar_t           *pras;
  562. int        w;
  563. int        h;
  564. uchar_t           *ptrans;
  565. {
  566.  
  567. /*
  568.  *    median filter an image. return 1 if successful, else 0.
  569.  *
  570.  *    the basic algorithm is:
  571.  *
  572.  *        new[i] = median of neighboring cells
  573.  */
  574.  
  575.     uchar_t               *pr;
  576.     uchar_t               *pt;
  577.     long            x;
  578.     long            y;
  579.     uint_t            l[9]; 
  580.  
  581.  
  582.     pr = pras;
  583.     pt = ptrans;
  584.  
  585.     if (pras != ptrans)
  586.     {
  587.         /*
  588.          *   first copy orig image to transform image
  589.          */
  590.         copyrast (pr, w, h, pt);
  591.     }
  592.  
  593.  
  594.     /*
  595.      *   do it...
  596.      */
  597.     for (y = 1; y < h-1; y++)
  598.     {
  599.         for (x = 1; x < w-1; x++)
  600.         {
  601.  
  602.             l[0] = (uint_t) pt[(y-1)*w+(x-1)];
  603.             l[1] = (uint_t) pt[(y-1)*w+(x  )];
  604.             l[2] = (uint_t) pt[(y-1)*w+(x+1)];
  605.             l[3] = (uint_t) pt[(y  )*w+(x-1)];
  606.             l[4] = (uint_t) pt[(y  )*w+(x  )];
  607.             l[5] = (uint_t) pt[(y  )*w+(x+1)];
  608.             l[6] = (uint_t) pt[(y+1)*w+(x-1)];
  609.             l[7] = (uint_t) pt[(y+1)*w+(x  )];
  610.             l[8] = (uint_t) pt[(y+1)*w+(x+1)];
  611.  
  612.             piksrt (9, l);
  613.  
  614.             pt[(y*w) + x] = l[4];
  615.         }
  616.     }
  617.  
  618.     return (1);
  619. }
  620.  
  621.  
  622.  
  623.  
  624. /*------------------------------*/
  625. /*    piksrt            */
  626. /*------------------------------*/
  627. piksrt (n, arr)
  628. register int        n;
  629. register uint_t           *arr;
  630. {
  631.  
  632. /*
  633.  *    straight insertion sort (num. rec. in C, pp 243). is N^2 sort,
  634.  *    probably ok for N=9. shell sort is 3x faster for N=9.
  635.  */
  636.  
  637.     register int    i;
  638.     register int    j;
  639.     register uint_t    a;
  640.  
  641.     for (j = 1; j < n; j++)
  642.     {
  643.         for (a = arr[j], i = j - 1; i && arr[i] > a; i--)
  644.             arr[i+1] = arr[i];
  645.  
  646.         arr[i+1] = a;
  647.     }
  648. }
  649.  
  650.  
  651.  
  652.  
  653. /*------------------------------*/
  654. /*    logscale        */
  655. /*------------------------------*/
  656. int logscale (pras, w, h, ptrans)
  657. uchar_t           *pras;
  658. int        w;
  659. int        h;
  660. uchar_t           *ptrans;
  661. {
  662.  
  663. /*
  664.  *    apply log scaling to an image. return 1 if successful, else 0.
  665.  *
  666.  *    the basic algorithm is:
  667.  *
  668.  *        new[i] = maxval * log(old[i]) / log(maxval)
  669.  *        new[i] = maxval * log2(old[i]) / log2(maxval)
  670.  */
  671.  
  672.     uchar_t               *pr;
  673.     uchar_t               *pt;
  674.     ulong_t            maxval;
  675.     ulong_t            l2maxval;
  676.     ulong_t            l2pr;
  677.     uchar_t            uctmp;
  678.     long            x;
  679.     long            y;
  680.     long            val;
  681.     register long        ii;
  682.     register long        lim;
  683.     register uchar_t       *ps;
  684.  
  685.  
  686.  
  687.     pr = pras;
  688.     pt = ptrans;
  689.  
  690.     if (pras != ptrans)
  691.     {
  692.         /*
  693.          *   first copy orig image to transform image
  694.          */
  695.         copyrast (pr, w, h, pt);
  696.     }
  697.  
  698.  
  699.     /*
  700.      *   find maxval
  701.      */
  702. #if 0
  703.     maxval = 0L;
  704.     for (y = 0; y < h; y++)
  705.     {
  706.         for (x = 0; x < w; x++)
  707.         {
  708.             uctmp = pt[(y*w) + x];
  709.             if (uctmp > maxval)
  710.                 maxval = (ulong_t) uctmp;
  711.         }
  712.     }
  713. #else
  714.     maxval = 255L;
  715. #endif
  716.  
  717.     /*
  718.      *   do it...
  719.      */
  720. #if 0
  721.     l2maxval = (ulong_t) Log2 ((uint_t) maxval);
  722.     if (l2maxval == 0)
  723.         return (0);
  724. #else
  725.     l2maxval = 16;
  726. #endif
  727.  
  728. #ifdef DBL_LOOP
  729.     for (y = 0; y < h; y++)
  730.     {
  731.         for (x = 0; x < w; x++)
  732.         {
  733.             l2pr          = (ulong_t) Log2 ((uint_t) pt[(y*w)+x]);
  734.             pt[(y*w) + x] = (uchar_t) (maxval * l2pr / l2maxval);
  735.         }
  736.     }
  737. #else
  738.     lim = (long) w * (long) h;
  739.     for (ps = pt, ii = 0L; ii < lim; ii++)
  740.     {
  741.         l2pr  = (ulong_t) Log2 ((uint_t) *ps);
  742. /*        *ps++ = (uchar_t) (maxval * l2pr / l2maxval);*/
  743.         *ps++ = (uchar_t) ((maxval * l2pr) >> 4);
  744.     }
  745. #endif
  746.  
  747.     return (1);
  748. }
  749.  
  750.  
  751.  
  752.  
  753.  
  754. /*------------------------------*/
  755. /*    Log2            */
  756. /*------------------------------*/
  757. int Log2(x)
  758. uint_t    x;
  759. {
  760. /*
  761.  *    VERY approximate log base 2 times 10. basically finds MS bit...
  762.  */
  763.     if (x & 0xff00 == 0)    goto lobyte;
  764.     if (x & 0x8000)        return (15);
  765.     if (x & 0x4000)        return (14);
  766.     if (x & 0x2000)        return (13);
  767.     if (x & 0x1000)        return (12);
  768.     if (x & 0x0800)        return (11);
  769.     if (x & 0x0400)        return (10);
  770.     if (x & 0x0200)        return (9);
  771.     if (x & 0x0100)        return (8);
  772. lobyte:
  773.     if (x == 0x00ff)
  774.         return (80);
  775.     if (x == 0x0000)
  776.         return (1);
  777.     if (x & 0x0080)
  778.     {
  779.         if (x & 0x0040)
  780.             return (78);
  781.         else if (x & 0x0020)
  782.             return (75);
  783.         else
  784.             return (71);
  785.     }
  786.     if (x & 0x0040)
  787.     {
  788.         if (x & 0x0020)
  789.             return (68);
  790.         else if (x & 0x0010)
  791.             return (65);
  792.         else
  793.             return (61);
  794.     }
  795.     if (x & 0x0020)
  796.     {
  797.         if (x & 0x0010)
  798.             return (58);
  799.         else if (x & 0x0008)
  800.             return (55);
  801.         else
  802.             return (51);
  803.     }
  804.     if (x & 0x0010)
  805.     {
  806.         if (x & 0x0008)
  807.             return (48);
  808.         else if (x & 0x0004)
  809.             return (45);
  810.         else
  811.             return (41);
  812.     }
  813.     if (x & 0x0008)
  814.     {
  815.         if (x & 0x0004)
  816.             return (38);
  817.         else if (x & 0x0002)
  818.             return (35);
  819.         else
  820.             return (31);
  821.     }
  822.     if (x & 0x0004)
  823.     {
  824.         if (x & 0x0002)
  825.             return (28);
  826.         else if (x & 0x0001)
  827.             return (25);
  828.         else
  829.             return (21);
  830.     }
  831.     if (x & 0x0002)
  832.     {
  833.         if (x & 0x0001)
  834.             return (18);
  835.         else
  836.             return (13);
  837.     }
  838.     if (x & 0x0001)        return (10);
  839.     return (0);
  840. }
  841.  
  842.  
  843.  
  844.  
  845. /*------------------------------*/
  846. /*    contrast        */
  847. /*------------------------------*/
  848. int contrast (pras, w, h, thresh, hist, ptrans)
  849. uchar_t           *pras;
  850. int        w;
  851. int        h;
  852. long        thresh;
  853. long           *hist;
  854. uchar_t           *ptrans;
  855. {
  856.  
  857. /*
  858.  *    apply contrast expansion to an image. return 1 if successful, else 0.
  859.  *
  860.  *    the basic algorithm is:
  861.  *
  862.  *        new[i] = maxval * (old[i] - lo) / (hi - lo)
  863.  *
  864.  *        where lo to hi is new contrast range
  865.  */
  866.  
  867.     uchar_t               *pr;
  868.     uchar_t               *pt;
  869.     ulong_t            maxval;
  870.     uchar_t            uctmp;
  871.     long            x;
  872.     long            y;
  873.     long            hi_lo;
  874.     long            newlo;
  875.     long            newhi;
  876.     int            i;
  877.     register long        val;
  878.     register long        rng;
  879.     register long        lim;
  880.     register long        ii;
  881.     register uchar_t       *ps;
  882.  
  883.  
  884.     pr = pras;
  885.     pt = ptrans;
  886.  
  887.     if (pras != ptrans)
  888.     {
  889.         /*
  890.          *   first copy orig image to transform image
  891.          */
  892.         copyrast (pr, w, h, pt);
  893.     }
  894.  
  895.  
  896.     /*
  897.      *   get current histogram...
  898.      */
  899.     if (do_hist (pt, w, h, hist) == 0)
  900.         return (0);
  901.  
  902.  
  903.     /*
  904.      *    find intensities on either side of hist where pixel count
  905.      *    exceeds thresh
  906.      */
  907.     for (i = 0; i < HISTSIZ; i++)
  908.     {
  909.         if (hist[i] > thresh)
  910.             break;
  911.     }
  912.     newlo = i;
  913.     for (i = HISTSIZ-1; i > newlo; i--)
  914.     {
  915.         if (hist[i] > thresh)
  916.             break;
  917.     }
  918.     newhi = i;
  919.  
  920.  
  921.     /*
  922.      *   do it. set pixels below lower threshold 0 and those above to 255.
  923.      *   contrast expand pixels between...
  924.      */
  925. #ifdef DBL_LOOP
  926.     rng = newhi - newlo + 1;
  927.     for (y = 0; y < h; y++)
  928.     {
  929.         for (x = 0; x < w; x++)
  930.         {
  931.             val = (ulong_t) pt[(y*w) + x];
  932.             if (val < newlo)
  933.                 pt[(y*w) + x] = 0;
  934.             else if (val > newhi)
  935.                 pt[(y*w) + x] = 255;
  936.             else
  937.             {
  938.                 pt[(y*w) + x] = (uchar_t) ((255L * (val - newlo)) / rng);
  939.             }
  940.         }
  941.     }
  942. #else
  943.     rng = newhi - newlo + 1;
  944.     lim = (long) w * (long) h;
  945.     for (ps = pt, ii = 0L; ii < lim; ii++, ps++)
  946.     {
  947.         val = (ulong_t) *ps;
  948.         if (val < newlo)
  949.             *ps = 0;
  950.         else if (val > newhi)
  951.             *ps = 255;
  952.         else
  953.             *ps = (uchar_t) ((255L * (val - newlo)) / rng);
  954.     }
  955. #endif
  956.  
  957.     return (1);
  958. }
  959.  
  960.  
  961.  
  962.  
  963. /*------------------------------*/
  964. /*    brighten        */
  965. /*------------------------------*/
  966. int brighten (pras, w, h, brite, ptrans)
  967. uchar_t           *pras;
  968. int        w;
  969. int        h;
  970. register int    brite;
  971. uchar_t           *ptrans;
  972. {
  973.  
  974. /*
  975.  *    brighten an image. return 1 if successful, else 0.
  976.  *
  977.  *    the basic algorithm is:
  978.  *
  979.  *        new[i] = old[i] + brite
  980.  */
  981.  
  982.     uchar_t               *pr;
  983.     uchar_t               *pt;
  984.     long            x;
  985.     long            y;
  986.     int            sign;
  987.     register uint_t        val;
  988.     register long        ii;
  989.     register long        lim;
  990.     register uchar_t       *ps;
  991.  
  992.  
  993.  
  994.     pr = pras;
  995.     pt = ptrans;
  996.  
  997.     if (pras != ptrans)
  998.     {
  999.         /*
  1000.          *   first copy orig image to transform image
  1001.          */
  1002.         copyrast (pr, w, h, pt);
  1003.     }
  1004.  
  1005.  
  1006.     sign = 0;
  1007.     if (brite < 0)
  1008.     {
  1009.         sign = 1;
  1010.         brite = -brite;
  1011.     }
  1012.  
  1013.  
  1014.     /*
  1015.      *   do it...
  1016.      */
  1017.     if (sign)
  1018.     {
  1019. #ifdef DBL_LOOP
  1020.         for (y = 0; y < h; y++)
  1021.         {
  1022.             for (x = 0; x < w; x++)
  1023.             {
  1024.                 val = (uint_t) pt[y*w + x];
  1025.                 if (val < brite)
  1026.                     val  = 0;
  1027.                 else
  1028.                     val -= brite;
  1029.                 pt[(y*w) + x] = (uchar_t) val;
  1030.             }
  1031.         }
  1032. #else
  1033.         lim = (long) w * (long) h;
  1034.         for (ps = pt, ii = 0L; ii < lim; ii++, ps++)
  1035.         {
  1036.             if ((uint_t) *ps < brite)
  1037.                 *ps  = 0;
  1038.             else
  1039.                 *ps -= brite;
  1040.         }
  1041. #endif
  1042.     }
  1043.     else
  1044.     {
  1045. #ifdef DBL_LOOP
  1046.         for (y = 0; y < h; y++)
  1047.         {
  1048.             for (x = 0; x < w; x++)
  1049.             {
  1050.                 val = (uint_t) ((uint_t) pt[y*w + x] + brite);
  1051.                 if (val > 255)
  1052.                     val = 255;
  1053.                 pt[(y*w) + x] = (uchar_t) val;
  1054.             }
  1055.         }
  1056. #else
  1057.         lim = (long) w * (long) h;
  1058.         for (ps = pt, ii = 0L; ii < lim; ii++, ps++)
  1059.         {
  1060.             val = (uint_t) ((uint_t) *ps + brite);
  1061.             if (val > 255)
  1062.                 val = 255;
  1063.             *ps = (uchar_t) val;
  1064.         }
  1065. #endif
  1066.     }
  1067.  
  1068.     return (1);
  1069. }
  1070.  
  1071.  
  1072.  
  1073.  
  1074. /*------------------------------*/
  1075. /*    invert            */
  1076. /*------------------------------*/
  1077. int invert (pras, w, h, thresh, ptrans)
  1078. uchar_t           *pras;
  1079. int        w;
  1080. int        h;
  1081. int        thresh;
  1082. uchar_t           *ptrans;
  1083. {
  1084.  
  1085. /*
  1086.  *    invert (negate) an image. return 1 if successful, else 0.
  1087.  *
  1088.  *    the basic algorithm is:
  1089.  *
  1090.  *        new[i] = 255 - old[i]
  1091.  */
  1092.  
  1093.     uchar_t               *pr;
  1094.     uchar_t               *pt;
  1095.     long            x;
  1096.     long            y;
  1097.     int            sign;
  1098.     register long        lim;
  1099.     register long        ii;
  1100.     register int        thr;
  1101.     register uchar_t       *ps;
  1102.  
  1103.  
  1104.     pr = pras;
  1105.     pt = ptrans;
  1106.  
  1107.     if (pras != ptrans)
  1108.     {
  1109.         /*
  1110.          *   first copy orig image to transform image
  1111.          */
  1112.         copyrast (pr, w, h, pt);
  1113.     }
  1114.  
  1115.  
  1116.     sign = 0;
  1117.     if (thresh < 0)
  1118.     {
  1119.         sign   = 1;
  1120.         thresh = -thresh;
  1121.     }
  1122.  
  1123.  
  1124.     /*
  1125.      *   do it...
  1126.      */
  1127. #ifdef DBL_LOOP
  1128.     if (thresh == 0)
  1129.     {
  1130.         for (y = 0; y < h; y++)
  1131.         {
  1132.             for (x = 0; x < w; x++)
  1133.             {
  1134.                 pt[(y*w) + x] = 255 - pt[(y*w) + x];
  1135.             }
  1136.         }
  1137.     }
  1138.     else if (sign)
  1139.     {
  1140.         /*
  1141.          *   invert only below thresh
  1142.          */
  1143.         for (y = 0; y < h; y++)
  1144.         {
  1145.             for (x = 0; x < w; x++)
  1146.             {
  1147.                 if ((uint_t) pt[(y*w) + x] < thresh)
  1148.                     pt[(y*w) + x] = 255 - pt[(y*w) + x];
  1149.             }
  1150.         }
  1151.     }
  1152.     else
  1153.     {
  1154.         /*
  1155.          *   invert only above thresh
  1156.          */
  1157.         for (y = 0; y < h; y++)
  1158.         {
  1159.             for (x = 0; x < w; x++)
  1160.             {
  1161.                 if ((uint_t) pt[(y*w) + x] > thresh)
  1162.                     pt[(y*w) + x] = 255 - pt[(y*w) + x];
  1163.             }
  1164.         }
  1165.     }
  1166. #else
  1167.     lim = (long) h * (long) w;
  1168.     thr = thresh;
  1169.     if (thr == 0)
  1170.     {
  1171.         for (ps = pt, ii = 0L; ii < lim; ii++, ps++)
  1172.         {
  1173.             *ps = 255 - *ps;
  1174.         }
  1175.     }
  1176.     else if (sign)
  1177.     {
  1178.         /*
  1179.          *   invert only below thresh
  1180.          */
  1181.         for (ps = pt, ii = 0L; ii < lim; ii++, ps++)
  1182.         {
  1183.             if ((uint_t) *ps < thr)
  1184.                 *ps = 255 - *ps;
  1185.         }
  1186.     }
  1187.     else
  1188.     {
  1189.         /*
  1190.          *   invert above below thresh
  1191.          */
  1192.         for (ps = pt, ii = 0L; ii < lim; ii++, ps++)
  1193.         {
  1194.             if ((uint_t) *ps > thr)
  1195.                 *ps = 255 - *ps;
  1196.         }
  1197.     }
  1198. #endif
  1199.  
  1200.     return (1);
  1201. }
  1202.  
  1203.  
  1204.  
  1205.  
  1206. /*------------------------------*/
  1207. /*    threshold        */
  1208. /*------------------------------*/
  1209. int threshold (pras, w, h, thresh, ptrans)
  1210. uchar_t           *pras;
  1211. int        w;
  1212. int        h;
  1213. int        thresh;
  1214. uchar_t           *ptrans;
  1215. {
  1216.  
  1217. /*
  1218.  *    threshold an image. return 1 if successful, else 0.
  1219.  *
  1220.  *    the basic algorithm is:
  1221.  *
  1222.  *        new[i] =   0 if old[i] <= thresh
  1223.  *        new[i] = 255 if old[i] >  thresh
  1224.  */
  1225.  
  1226.     uchar_t               *pr;
  1227.     uchar_t               *pt;
  1228.     long            x;
  1229.     long            y;
  1230.     uint_t            val;
  1231.     register long        ii;
  1232.     register long        lim;
  1233.     register uchar_t       *ps;
  1234.  
  1235.  
  1236.     pr = pras;
  1237.     pt = ptrans;
  1238.  
  1239.     if (pras != ptrans)
  1240.     {
  1241.         /*
  1242.          *   first copy orig image to transform image
  1243.          */
  1244.         copyrast (pr, w, h, pt);
  1245.     }
  1246.  
  1247.  
  1248.     /*
  1249.      *   do it...
  1250.      */
  1251. #ifdef DBL_LOOP
  1252.     for (y = 0; y < h; y++)
  1253.     {
  1254.         for (x = 0; x < w; x++)
  1255.         {
  1256.             val = (uint_t) pt[(y*w) + x];
  1257.             if (val > thresh)
  1258.                 pt[(y*w) + x] = 255;
  1259.             else
  1260.                 pt[(y*w) + x] = 0;
  1261.         }
  1262.     }
  1263. #else
  1264.     lim = (long) h * (long) w;
  1265.     for (ps = pt, ii = 0L; ii < lim; ii++, ps++)
  1266.     {
  1267.         if ((uint_t) *ps > thresh)
  1268.             *ps = 255;
  1269.         else
  1270.             *ps = 0;
  1271.     }
  1272. #endif
  1273.  
  1274.     return (1);
  1275. }
  1276.  
  1277.  
  1278.  
  1279.  
  1280. /*------------------------------*/
  1281. /*    histeq            */
  1282. /*------------------------------*/
  1283. int histeq (pras, w, h, hist, ptrans)
  1284. uchar_t           *pras;
  1285. int        w;
  1286. int        h;
  1287. long           *hist;
  1288. uchar_t           *ptrans;
  1289. {
  1290.  
  1291. /*
  1292.  *    histogram equalization of image. return 1 if successful, else 0.
  1293.  */
  1294.  
  1295.     long            p_r[HISTSIZ];
  1296.     long            s_k[HISTSIZ];
  1297.     long            s_kapprx[HISTSIZ];
  1298.     uchar_t               *pr;
  1299.     uchar_t               *pt;
  1300.     int            i;
  1301.     long            count;
  1302.  
  1303.  
  1304.  
  1305.     pr = pras;
  1306.     pt = ptrans;
  1307.  
  1308.     if (pras != ptrans)
  1309.     {
  1310.         /*
  1311.          *   first copy orig image to transform image
  1312.          */
  1313.         copyrast (pr, w, h, pt);
  1314.     }
  1315.  
  1316.  
  1317.     /*
  1318.      *   get current histogram...
  1319.      */
  1320.     if (do_hist (pt, w, h, hist) == 0)
  1321.         return (0);
  1322.  
  1323.  
  1324.     /*
  1325.      *   find prob density function (p_r), transformation function (s_k),
  1326.      *   and new distribution. note the 1000 to avoid floats...
  1327.      */
  1328.     count = (long) w * (long) h;
  1329.     for (i = 0; i < HISTSIZ; i++)
  1330.     {
  1331.         p_r[i] = (hist[i] * 1000L) / count;
  1332.     }
  1333.     for (s_k[0] = p_r[0], i = 1; i < HISTSIZ; i++)
  1334.     {
  1335.         s_k[i] = s_k[i-1] + p_r[i];
  1336.     }
  1337.     for (i = 0; i < HISTSIZ; i++)
  1338.     {
  1339.         /* the 500 is for rounding to nearest int... */
  1340.         s_kapprx[i] = ((s_k[i] * 255L) + 500L) / 1000L;
  1341.     }
  1342.  
  1343. #if 0
  1344.     printf ("\n    i      hist       p_r       s_k  s_kapprx\n");
  1345.     for (count = 0, i = 0; i < HISTSIZ; i++)
  1346.     {
  1347.         printf ("%5d%10ld%10ld%10ld%10ld\n",
  1348.             i,hist[i],p_r[i],s_k[i],s_kapprx[i]);
  1349.     }
  1350. #endif
  1351.  
  1352.     /*
  1353.      *   now redistribute. s_kapprx will contain the new value for
  1354.      *   a given pixel intensity
  1355.      */
  1356.     redistribute (pt, w, h, s_kapprx);
  1357.  
  1358.     return (1);
  1359. }
  1360.  
  1361.  
  1362.  
  1363.  
  1364. /*------------------------------*/
  1365. /*    redistribute        */
  1366. /*------------------------------*/
  1367. int redistribute (pras, w, h, s_kapprx)
  1368. uchar_t           *pras;
  1369. int        w;
  1370. int        h;
  1371. register long  *s_kapprx;
  1372. {
  1373.  
  1374. /*
  1375.  *    redistribute pixels. each pixel set to value given by s_kapprx
  1376.  */
  1377.  
  1378.     register uchar_t       *pr;
  1379.     register long        ii;
  1380.     register long        lim;
  1381.  
  1382.  
  1383.     lim = (long) w * (long) h;
  1384.     for (pr = pras, ii = 0L; ii < lim; ii++, pr++)
  1385.     {
  1386.         *pr = (uchar_t) s_kapprx[(uint_t) *pr];
  1387.     }
  1388. }
  1389.  
  1390.  
  1391.  
  1392.  
  1393. /*------------------------------*/
  1394. /*    copyrast        */
  1395. /*------------------------------*/
  1396. int copyrast (ps, w, h, pd)
  1397. uchar_t           *ps;        /* source */
  1398. int        w;
  1399. int        h;
  1400. uchar_t           *pd;        /* dest */
  1401. {
  1402.  
  1403. /*
  1404.  *    copy an image. return 1 if successful, else 0.
  1405.  */
  1406.  
  1407. #ifdef DBL_LOOP
  1408.     long        x;
  1409.     long        y;
  1410.  
  1411.     for (y = 0; y < h; y++)
  1412.     {
  1413.         for (x = 0; x < w; x++)
  1414.         {
  1415.             pd[(y*w) + x] = ps[(y*w) + x];
  1416.         }
  1417.     }
  1418. #else
  1419.     register long    ii;
  1420.     register long    lim;
  1421.  
  1422.     lim = (long) w * (long) h;
  1423.     for (ii = 0L; ii < lim; ii++)
  1424.         *pd++ = *ps++;
  1425. #endif
  1426.     return (1);
  1427. }
  1428.  
  1429.